こがねいろ

いろいろ備忘録

反復測定分散分析と変量効果

はじめに

反復測定分散分析について頭が混乱してきたのでまとめる。


参考文献

千野 直仁（2003）．反復測定データの分析とその周辺　教育心理学年報，42，107-118．
https://doi.org/10.5926/arepj1962.42.0_107

その他ウェブ上の情報など
http://www.cardio.med.tohoku.ac.jp/2005/newmember/pdf/ms/15_6S.pdf

www.slideshare.net


ノート

反復測定分散分析の概要

  • 同一個体（被験者など）から複数の変数を単一時点または複数時点で得たデータに対して利用できるモデル（ほかにGMANOVAなどもある）。
  • 伝統的な分散分析におけるF分布や（対比検定を行う際の）t分布は，反復測定デザインの場合，データにより歪む可能性がある。
  • 球面性仮定（球状性仮定，球形仮定）が満たされている場合には，検定量は歪まないとされる。
    • 球面性：反復測定変数の水準間の差の分散が等しいこと
    • 経時的な反復測定変数の場合，時点が離れるほど相関は低下するので，球面性が満たされないことも多い？

反復測定分散分析の弱点

  • 原則的には母数モデルが仮定されている。
    • 母数モデル：誤差項を除く未知パラメータがすべて定数［母数］であるモデル
    • 個体による差を固定効果で表現している（パラメータの数が増えてしまう）ということ？
    • 使用した刺激の効果など，変量効果で表現するべき効果を表現しきれない場合がある？
  • 反復測定変数の中に欠損値があると，欠損値を含む個体ごと除外する必要がある。
  • 経時的な反復測定変数は連続量ではなく，カテゴリカルな水準として扱われる？

G and R side effect

はじめに

線形混合モデルにおいて経時的な反復測定データを扱う場合，例えば「個体ごとのランダム効果」と「時点間の相関構造」を同時に組み込みたい場合がある。

SAS（SAS PROC Glimmix）では，前者をG side effect，後者をR side effectとし，両者をモデルに組み込むことができるようである。

SPSS（線型混合モデル）でも，前者を変数選択により，後者を「反復測定共分散」により指定することができる。

一方で，Rではパッケージにより指定可能かどうか，指定の方法が異なる。


混合モデルを扱うRの関数

{nlme}
誤差にガウス分布しか指定できない。
共分散構造を組み込むことができる。
（反復測定間の相関を指定できる）

{lme4}
リンク関数が豊富である。
共分散構造を組み込むことができない（難しい？）。

{brms}
共分散構造を組み込むことができる。
ベイジアンアプローチを採用している。
Rにおけるモデルの構文は {lme4} と同様のもの。


関連する記事・サイトのURL

{nlme} と {lme4} の選択
stackexchange How to choose nlme or lme4 R library for mixed effects models? - Cross Validated

空間的に自己相関のある残差への対処
stackoverflow r - How to deal with spatially autocorrelated residuals in GLMM - Stack Overflow

空間的自己相関を考量したベイズモデル
久保氏 簡単な例題で理解する空間統計モデル : HUSCAP

SPSSにおける共分散構造の指定
線型混合モデル

{nlme} における共分散構造の指定
橋本氏のウェブページ 橋本洸哉のサイト - ランダム効果とcorrelation structure

{brms} における共分散構造の指定
Specifying R-side and G-side Random Effects in brms - brms - The Stan Forums

縦断的データの解析手法に関する講義
Marie Davidian氏 ST 732 - Spring 2022

エラーなく実験が終了してしまう（ように見える）現象 [PsychoPy]

PsychoPy

概要と原因

PsychoPyで実験を作成していて，エラー表示がなく実験が終了してしまう「ように見える」ことがあります。

実際は正常に実験が終了しているのですが，ユーザー側には実験がエラーを吐かずに終了してしまうように感じられる場合を紹介します。

例えば，Routine-AとBが連続しているとき，それぞれのRoutine内に置かれたJoystickコンポーネントのボタン入力に応じてRoutineを終了させる場合を考えます。

# In Routine_A
# Press a button to quit Routine_A
if js01.getAllButtons()[0]:
    continueRoutine = False

# In Routine_B
# Press a button to quit Routine_B
if js02.getAllButtons()[0]:
    continueRoutine = False

このとき，Routine-Aで指定したボタンの入力を検出すると，Routine-Bへと移行します。

しかし，人間側のボタン押しは，多くの場合Routine-Bに移行した後も継続してしまいます。

（※ 人間が「一瞬の間」ボタンを押して離すまでの間に，コンピュータはRoutine-AからBへと移行してしまうため，Routine-Bでもキーを押していることになります。）

すると，Routine-Bにおけるif文の処理により，Routine-Bも即座に終了します。

このとき，人間側はRoutine-Bが「スキップされたような印象」を受けます。

これがエラーなく実験が中断されてしまうように見える現象の一因になります。


対策

簡単な対策としては，「Routine-Bでキーを検出しはじめるまでに少しだけタイムラグを持たせる」とか，「Routine-AとBでRoutineを終了させるキーを変える」などが挙げられるでしょう。

あるいは，「キーの”押し”と”離し”を分離して評価できるようなコードに変える」こともできるでしょう。

他に何か対策があればご教示ください。

参考までに，キーを検出しはじめるまでにタイムラグ（5秒）を持たせるコードを記しておきます。

# In Routine_B
# Press a button to quit Routine_B
if Routine_BClock.getTime() > 5:
    if js02.getAllButtons()[0]:
        continueRoutine = False


# The PsychoPy experiment ends without error (or so it seems)

日付の表記方法の変更 [PsychoPy]

PsychoPy Python

はじめに

PsychoPyのデフォルトの出力では，月の名前を省略した英名で出力します。

そのため，例えば3月と4月に実験を行うと，"Mar"と"Apr"という単語が混在することになります。

しかし，Windowsのエクスプローラー上では，多くの場合アルファベット順にファイルが整理されますので，4月のファイル"Apr"の方が3月"Mar"のファイルよりも上にくることになります。

この問題の解決策を考えます。

月名の省略
月 英語 省略系
1月 January Jan
2月 February Feb
3月 March Mar
4月 April Apr
5月 May May
6月 June Jun
7月 July Jul
8月 August Aug
9月 September Sep
10月 October Oct
11月 November Nov
12月 December Dec


2月にまたがる場合を考えると，

1-2月，3-4月，5-6月，6-7月，7-8月，9-10月，10-11月，11-12月

がアルファベット順の逆転現象が起こることになるでしょうか。

意外と多いですね。


ファイル名の指定

PsychoPyのBuilderから，プロパティ（歯車マーク）を開きます。

データタブにおいて「データファイル名」を指定できますので，expInfo['date']を次のように書き換えます。

# before
u'data/%s_%s_%s' % (expInfo['participant'], expName, expInfo['date'])

# after
u'data/%s_%s_%s' % (expInfo['participant'], expName, data.getDateStr(format="%Y%m%d-%H%M"))

参考記事

Changing the date format in the filename (%b to %m) - Builder - PsychoPy

時系列データの保存 [PsychoPy]

PsychoPy Python

はじめに

PsychoPyで，デフォルトの出力以外に実験結果を保存する方法を考えます。
Builderで基本的な骨組みを作り，一部をCoderで補完する方法をとります。

デフォルトの出力

PsychoPyでは，実験中に取得したデータを実験終了時に書き出すことができます。

例えば，「刺激」コンポーネント（例：Textコンポーネント）で次のように設定しておくと，刺激呈示の開始・終了時刻が最後に出力されます。

f:id:gh_nishiki:20210904171408p:plain
Textコンポーネント

「反応」コンポーネント（例：Keyboardコンポーネント）で次のような設定を行うと，押されたキーを保存することができます。

f:id:gh_nishiki:20210904171210p:plain
keyboardコンポーネント

「反応」コンポーネントのうち，一部のコンポーネント（例：Mouseコンポーネント）では，次のような設定でフレームごとの状態を保存することができます。

f:id:gh_nishiki:20210904172040p:plain
Mouseコンポーネント

概要

しかし，実験によっては，より応用的な処理が必要になることがあります。

例えば，ジョイスティックでスライダー（評定尺度のような横線）を操作して，その評定値をフレームごとに取得したいといった場合です。

その場合，ジョイスティックの状態（例：x軸方向の傾き）を取得するのではなく，スライダーの評定値を保存できれば便利です。

今回はその方法を考えていきます。

空の配列を用意する

Builderで，カスタムコンポーネントからCodeコンポーネントを追加しましょう。

Codeコンポーネントを開き，「フレーム毎」より前に処理されるタブ（「実験初期化中」または「実験開始時」）を選択します。

コード欄に次のように書き込み，空の配列を作成しましょう。

今回は，スライダーの評定値（slider_x）とデータを取得した時刻（slider_t）を入れておくために，2つの配列を作成しています。

slider_x = []
slider_t = []
フレームごとに状態を取得して配列に格納する

続いて，Codeコンポーネントの「フレーム毎」のタブを選択し，次のコードを書き込みます。

フレーム毎にスライダーの状態と，現在の時刻を取得して配列に格納しています。

slider_x.append(YourSliderName.markerPos)
slider_t.append(YourRoutine-ClockName.getTime())


ここで，”YourSliderName”と"YourRoutine-ClockName"という名前が出てきましたが，この箇所には各人が実験で使用している変数の名前が入ります。

評定値を取得したいSliderコンポーネントに，"my_sl01"という名前を付けている場合は，"my_sl01.markerPos"と書くことで"my_sl01"のマーカーの位置を取得することができます。

また，実験のRoutineに"my_R01"という名前を付けている場合は，"my_R01Clock.getTime()"でRoutine開始からの経過時間を取得することができます。

（※ "Clock"と"getTime"の大文字-小文字に注意しましょう。）


詳しくは，下記のサイトを参照してください。
スライダーのマーカー：https://psychopy.org/api/visual/slider.html
ルーチンの経過時間：Pythonで心理実験 - 例題18-5 — 十河研究室


データを書き出す

ここまでの処理で，Pythonの配列の中に①スライダーの評定値と②それに対応する時刻が入れられている状態になりました。

最後にそのデータを書き出します。

データの書き出しは，必ずしもデータを取得したCodeコンポーネント内で行う必要はありません。

データ取得以降であれば，別のRoutineで行っても差し支えありません。

今回は，同一Codeコンポーネント内で書き出しまで行ってしまいましょう。


Codeコンポーネントの「Routine修了時」のタブを選択し，次のコードを書き込みます。

# ライブラリのインポート
import csv

# データの書き込み
my_data = [list(x) for x in list(zip(slider_t, slider_x))]
path = r'data/{}_{}_{}_additional.csv'.format(data.getDateStr(format="%Y%m%d-%H%M"), expInfo['participant'], expName)
with open(path, 'w', newline='') as f:
    w = csv.writer(f)
    w.writerows(my_data)


はじめに，ライブラリ（csv）をインポートしています。

これはcsvファイルの書き出しに使用しています。

（※ インポートはいつでも構わないので，基本的に私は管理しやすいように実験の冒頭で行います。）


続いて，zipを用いて，1次元の配列2個を一つの2次元配列にまとめています。

「何秒目に，どの評定値か」をセットして縦に（列方向に）並べているイメージでしょうか。

この処理を挟むことで，末尾にある"writerows"を使用できます。


"path="の行では，csvファイルをコンピュータ上のどこに作成するかを指定しています。

デフォルト状態のPsychoPyは，実験プログラムが置かれているディレクトリ上にある"data"のフォルダにデータを出力するので，それに合わせて"data/"からパスを書き始めるとよいかもしれません。


なお，デフォルトのPsychoPyは，月の名前を省略した英名で出力します。

そのため，例えば3月と4月に実験を行うと，"Mar"と"Apr"という単語が混在することになります。

しかし，Windowsのエクスプローラー上では，多くの場合アルファベット順にファイルが整理されますので，4月のファイル"Apr"の方が3月"Mar"のファイルよりも上にくることになります。

この状態は直感に反するので，実験実行日時の出力を書き換えるために下記の書き方をしています。

data.getDateStr(format="%Y%m%d-%H%M")


また，"expInfo['participant']"ですが，実験開始時に入力する実験情報ダイアログに"participant"という項目を設定しない場合は，表現を適宜変更してください。


おわりに

以上の処理を行うことで，時系列データを書き出すことができました。

配列操作の備忘録 [Numpy + Built-in Functions]

Numpy Python

はじめに

様々なサイトでまとめられていますが，個人用に基本的な配列操作をまとめます。
pythonの組み込み型として使用できる配列と，numpyをインポートとして使用するnumpy arrayに関して記述します。

配列の作成

# ライブラリ
import numpy as np


# 配列の作成（1）基本
# Python組み込み型とnumpy array
l   = ['A', 'B', 'C', 'Z', 'Z']
lnp = np.array(['A', 'B', 'C', 'Z', 'Z'])

# 2次元配列（list of list）の場合
l_2d   = [['A', 'B'], ['C', 'D', 'E']]
lnp_2d = np.array([['A', 'B', 'X'], ['C', 'D', 'Y']])

# *** numpy arrayの場合はdtypeとshapeの指定が重要
lnp_2d = np.array([['A', 'B', 'X'], ['C', 'D']])  # エラー


# 配列の作成（2）繰り返しの配列
# 指定した個数（12個）だけ同一要素をくり返す配列
l_r1   = ['A'] * 12
lnp_r1 = np.full(12, 'A') # 指定した値をくり返す
lnp_r2 = np.ones(12)      # 1をくり返す
lnp_r3 = np.zeros(12)     # 0をくり返す


# 配列の作成（3）繰り返しの配列
# ある配列（l）の要素を指定した個数(12個)だけ繰り返す配列
l_r2 = [l[i % len(l)] for i in range(12)]


# 配列の作成（4）等差数列
# 等差数列（0以上，10以下の範囲：要素数12を指定）
lnp_s1 = np.linspace(0, 10, 12)
# 等差数列（0以上，10未満の範囲：要素数12を指定）
lnp_s2 = np.linspace(0, 10, 12, endpoint=False)
# 等差数列（0以上，10未満の範囲：公差0.2を指定）
lnp_s3 = np.arange(0, 10, 0.2)
# 等差数列（0以上，10未満の範囲：公差が整数の場合）
l_s3   = [i for i in range(0, 10, 2)]


# 配列の作成（5）空の配列
l_emp   = []
l_emp   = list()
lnp_emp = np.empty(12) # np.emptyでは要素数shapeを指定する

要素の指定

# 要素の指定・抽出（1）1次元
# インデックスを指定する（先頭が0）
l[0]
# インデックスを指定する（負の整数で末尾から数える）
l[-1]
# インデックスを指定する（1以上から3未満まで = インデックス1と2）
l[1:3]
# インデックスを指定する（3未満のすべて = インデックス0, 1, 2）
l[:3]
# インデックスを指定する（1以上のすべて = インデックス1, 2,...）
l[1:]
# インデックスを指定する（0以上から6未満まで2個ずつ数える＝1個飛ばしで）
l[0:6:2]
# 指定したインデックス以外の要素（numpy arrayを返す）
np.delete(l, 2)
np.delete(l, [2, 4])


# 要素の指定・抽出（2）2次元
# l_2dの0番目のリストの1番目の要素
l_2d[0][1]
# numpy arrayではインデックスによる指定が楽
# インデックスを指定する（行は1番目，列は2番目）
lnp_2d[1, 2]
# インデックスを指定する（行はすべて，列は2番目）
lnp_2d[:, 2]
# インデックスを指定する（行はすべて，列は1列飛ばし）
lnp_2d[:, ::2]

要素の検索

# 要素の検索
# 配列は指定した値を含むか否か（ブール値を返す）
'A' in l
# 配列は指定した値を何個含むか（個数を返す）
l.count('A')
# 配列は指定した値をどの位置に含むか（最初のインデックスだけを返す）
l.index('Z')
# 配列は指定した値をどの位置に含むか（すべてのインデックスを返す）
[ind for ind, item in enumerate(l) if item == 'Z']


# 配列の長さ
# 1次元配列の場合
len(l)
lnp.size
# 2次元のnumpy arrayに対して，shapeで行列数を取得できる
row, col = lnp_2d.shape

要素の追加と削除

# 要素の追加（1）末尾に追加する
l.append('X')
np.append(lnp, 'X')
# *** numpy配列の場合，appendはメソッドとして使えない
lnp.append('X')  # エラー


# 要素の追加（2）位置を指定して追加する
l.insert(1, 'Y')
np.insert(lnp, 1, 'Y')


# 要素の削除（1）全ての要素を削除する
l.clear()

# 要素の削除（2）指定した要素を削除する
# 指定した要素を削除する（リスト先頭から該当する1個のみ）
l.remove('Z')
# 指定した要素を削除する（リスト中からすべて）
[i for i in l if i != 'Z']

# *** numpy arrayではremoveメソッドは使えない
lnp.remove('Z')  # エラー
# 条件式を用いて表現できる（該当箇所すべてを取り除く）
lnp[~(lnp=='Z')]

# 行列のインデックスを指定して削除する
np.delete(lnp_2d, 1, axis=0)       # axis=0： 行を削除
np.delete(lnp_2d, [0, 2], axis=1)  # axis=1： 列を削除