MPTinRを使ってMPTモデリングした
Multinomial Processing Tree model（MPTモデリング）をした時に，Rパッケージである「MPTinR」（Singmann & Kellen, 2013）を使った。日本語資料がネットの海に無かったので，使った話を備忘録的に残しておく。
参考文献：Singmann, H., & Kellen, D. (2013). MPTinR: Analysis of multinomial processing tree models in R. Behavioral Research, 45, 560-575.
想定する実験
学習段階で形容詞（暖かい，寒い，丸い，鋭いetc...）が，赤色または青色で呈示された後，再認段階では提示される形容詞が，学習段階で赤色で提示された（Red）のか，青色で呈示された（Blue）のか，はたまた呈示されなかったのか（非呈示：Not）を聞かれる実験を想定，再認段階で以下のモデル・パラメータを仮定する。再認対象項目の性質は，学習段階での形容詞の呈示方法（要は再認課題の正解）である。参加者の回答は，それぞれの再認対象項目の性質での，参加者の判断カテゴリを表す。
それぞれのパラメータについて説明すると，dRed/dBlueは再認対象の項目が赤色/青色で呈示されたと正しく再認する割合，dNotは呈示されていない（非呈示）と正しく再認する割合である。bは再認対象の項目が呈示されたと判断するバイアスの割合で，gRedは再認対象の項目が赤色で呈示されたと推測する割合である。
このモデルは，例えば，学習段階の際に赤色で呈示された項目について，再認段階で赤色で呈示されたと再認することは，「赤色で呈示されたと正しく再認する」ことに加えて，正しく再認できなくても「呈示されたとバイアスで判断し，赤色で呈示されたと推測する」ことでも達成されるということを表している。
この時，各判断カテゴリの割合は以下のように表すことができるので，あとはそれぞれについて，実験をして実際に観測された割合を与えると，MPTinRが方程式を解いて最も近い割合が得られるようなパラメータを推定してくれる。
--------------------------------------------------------------------------------------
赤色で呈示された項目を
赤と判断する割合＝dRed + (1 - dRed)*b*gRed
青と判断する割合＝(1 - dRed)*b*(1 - gRed)
呈示されていないと判断する割合＝(1 - dRed)*(1 -b)
青色で呈示された項目を
赤と判断する割合＝(1 - dBlue)*b*gRed
青と判断する割合＝dBlue + (1 - dBlue)*b*(1 - gRed)
呈示されていないと判断する割合＝(1 - dBlue)*(1 -b)
呈示されていない項目を
赤と判断する割合＝(1 - dNot)*b*gRed
青と判断する割合＝(1 - dNot)*b*(1 - gRed)
呈示されていないと判断する割合＝dNot + (1 - dNot)*(1 -b)
--------------------------------------------------------------------------------------
MPTinRを使ってみる
準備
（まず，MPTinRをインストールして呼び出しておく）
まず，モデルを指定する。モデルの指定というか，各判断カテゴリの割合を表す方程式の未知数（パラメータ）が入っている辺（???）を木毎に記述する。今回，赤色で呈示された項目を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合の順に，以下のように指定した。木毎に一行空ける必要があるので，木と木の間，加えて最終行に空行が入っている。最終行に空行がないとエラーを吐かれる。
txtファイルで保存しておけばいい。今回はmodel.txtとして保存した。
#赤色で呈示された項目を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
dRed + (1 - dRed)*b*gRed
(1 - dRed)*b*(1 - gRed)
(1 - dRed)*(1 -b)
#青色で呈示された項目を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dBlue)*b*gRed
dBlue + (1 - dBlue)*b*(1 - gRed)
(1 - dBlue)*(1 -b)
#非呈示項目を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dNot)*b*gRed
(1 - dNot)*b*(1 - gRed)
dNot + (1 - dNot)*(1 -b) 次に，パラメータの制限を指定する。MPTinRでのパラメータの制限は，「<」か「=」を使って行うことが出来る。不等式による制限は，等式による制限の前に指定する必要がある。また，制限式の中で変数は何度も指定できるが，不等式・等式の左辺に指定することが出来るのは一度だけである。そして，1行につき一つの制限式が指定され，指定毎に一行空ける必要があるので，ここでも最終行に空行が入っている。最終行に空行がないとエラーを吐かれる。
今回，dRedとdBlueの値は同じだと制限し，restrict.txtとして保存した。
#dRed = dBlue = dNotのように3つ以上の変数を同時に指定することも可能
dRed = dBlue
最後にデータセットを準備する。上記のようにデータセットの列がモデルで指定した判断カテゴリの順番に対応している必要がある。行はここでは参加者となっている。参加者IDは要らないので，枠線の中の部分だけで良い。なお，上記のデータは架空のものを適当に入力した。
実際にMPTinRで分析
check.mpt("モデル名")で，モデルの概要が確認できる。
> check.mpt("model.txt")今回は
--------------------------------------------------------------------------------------
$probabilities.eq.1
[1] TRUE
$n.trees
[1] 3
$n.model.categories
[1] 9
$n.independent.categories
[1] 6
$n.params
[1] 5
$parameters
[1] "b" "dBlue" "dNot" "dRed" "gRed"
--------------------------------------------------------------------------------------
と返って来た。
$probabilities.eq.1は，それぞれの木の割合の合計が1になっているかをTRUE or Falseで教えてくれる。普通は1になるはずである。$n.treesは木の数で，今回は3個である。$n.model.categoriesは，モデルが持つ判断のカテゴリ数である。今回，再認項目の種類（3：赤，青，非呈示）*参加者の判断（3：赤，青，非呈示）なので9個である。$n.independent.categoriesは，モデルが提供する独立した判断のカテゴリ数で，今回は6個である。$n.paramsは，パラメータの数で，今回は5個である。$parametersは，パラメータの内容で，アルファベット順に表示されている。
fit.mpt(データ，"モデルファイル"，"制限ファイル")で，推定が行われる。今回は，パラメータを制限しないモデル，
results <- fit.mpt(dat, "model.txt")制限したモデル，
Restrictedresults <- fit.mpt(dat, "model.txt", "restrict.txt")それぞれについて推定してみる。
デフォルトでは5回推定した結果の内，一番良かった結果を報告してくれる。n.optimで推定回数を指定できる。後FIAはデフォルトでは算出されないが，fiaでサンプル数を指定すると，モンテカルロ法によって算出される。
集約された結果は，
results[["goodness.of.fit"]][["aggregated"]]
results[["parameters"]][["aggregated"]]
Restrictedresults[["goodness.of.fit"]][["aggregated"]]
Restrictedresults[["parameters"]][["aggregated"]]で確認出来る。今回は，
と返って来た。
G^2は，観察された割合と，推定したモデルから予測された割合がどれくらいずれているかを表している（Singmann & Kellen, 2013）ので，小さければ小さいほどよい。また，G^2は漸近的にχ二乗分布に従うため，χ二乗検定を行うことでG^2の大きさが0と有意に異なるか確認できる。今回，制限の有無に関わらずG^2は0と有意に異ならないため，適合度は良好だと言える。パラメータの推定値を見ると，制限の有無に関わらず値はほぼ同じで，赤色/青色/非呈示を正しく再認する割合は6割程度，バイアスによって呈示されたと判断する割合は7割程度，バイアスによって呈示されたと判断した後に赤色だったと推測する割合は5割程度だったと解釈できる。
パラメータに制限をかけることの正当性を評価するためには，select.mpt関数を使用する。
select.mpt(list(r = results, rr = Restrictedresults), output = "full")
#output = "standard"で比較結果のみが返ってくる。今回は
と返って来た。
～.bestの列には，それぞれのデータセットが比較された時に，最良の説明を提供する頻度が示されている。今回は，パラメータに制限がかかったモデルにおいて，AICでもBICでも10人中9人（＝9 - 0）で最良の説明が提供された。
なお，FIAを算出していればこの表に表示される。FIAだとパラメータ値に制限がかかり，取り得る値の範囲が狭い場合など，単純なモデルの方が高く評価される（Singmann & Kellen, 2013）らしい。また，FIAはサンプル数が少ないと，最悪の場合競合するモデルの複雑さの項の順位が逆になることが指摘されている（Heck, Moshagen & Erdfelder, 2014）。Heck, Moshagen & Erdfelder（2014）は現実的な対処法として，FIAを用いる場合に推奨される下限のサンプル数の算出を行っている。
パラメータ推定値を参加者内で比較
実際の研究では，パラメータ推定値を算出して終わりではなく，パラメータ推定値を条件間で比較したい場合が殆どであると思われるので，まだ少し続く。パラメータ推定値の比較，つまりパラメータ推定値が有意に異なるかどうかは，パラメータが異なるとしたモデルと，パラメータが同じであると制限をかけたモデルとの間で，モデルの適合度を比較する，あるいはG2統計量の差が有意に異なるかどうかによって解釈できる。つまり，条件毎に対応する，異なるパラメータの制約をモデルにかけた時の，それぞれモデルの当てはまりの良さを比較し，判断する。
今回，想定実験では「学習段階で形容詞（暖かい，寒い，にぎやかな，さびしいetc...）が赤色または青色で呈示された後，再認段階では提示される形容詞が，学習段階で赤色で提示された（Red）のか，青色で呈示された（Blue）のか，はたまた呈示されなかったのか（非呈示：Not）を聞かれる」としている。この時，「暖かい」「にぎやかな」といったポジティブな形容詞と，「寒い」「さびしい」といったネガティブな形容詞とで，パラメータの値が異なるかどうか検定する。
仮説として，「gRedの推定値は，ネガティブな形容詞とポジティブな形容詞を再認する場合とで異なり，ネガティブな形容詞よりも，ポジティブな形容詞に置いて大きい」について検定する。なお，根拠は何もないことを述べておく。
準備
検定するにあたり，まず，想定するモデルとパラメータを若干改変する。具体的には，最初に出した図のモデル（木3つ）をポジティブな形容詞，ネガティブな形容詞毎に作る。つまり，木が6つあるモデルが出来る。この時，パラメータがそのままでは意味がないので，検定を行いたいパラメータ，「gRed」についてそれぞれポジティブな形容詞の場合は末尾に「p」，ネガティブな形容詞の場合は末尾に[n」でもつけておく。今回はModel2.txtとして保存した。
#赤色で呈示されたポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
dRed + (1 - dRed)*b*gRedp
(1 - dRed)*b*(1 - gRedp)
(1 - dRed)*(1 -b)
#青色で呈示されたポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dBlue)*b*gRedp
dBlue + (1 - dBlue)*b*(1 - gRedp)
(1 - dBlue)*(1 -b)
#非呈示のポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dNot)*b*gRedp
(1 - dNot)*b*(1 - gRedp)
dNot + (1 - dNot)*(1 -b)
#赤色で呈示されたネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
dRed + (1 - dRed)*b*gRedn
(1 - dRed)*b*(1 - gRedn)
(1 - dRed)*(1 -b)
#青色で呈示されたネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dBlue)*b*gRedn
dBlue + (1 - dBlue)*b*(1 - gRedn)
(1 - dBlue)*(1 -b)
#非呈示のネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
(1 - dNot)*b*gRedn
(1 - dNot)*b*(1 - gRedn)
dNot + (1 - dNot)*(1 -b)次に，上記のモデルについて，制限を2つ作成する。今回はリストで指定する。まず作成するのは，dRed = dBlueという制限のみのものである。次に作成するのはdRed = dBlueに加えて，gRedn = gRedpであるとする制限である。前者は仮説を支持する制限，後者は仮説を支持しない制限である。
r1 <- list("dRed = dBlue")
r2 <- list("dRed = dBlue", "gRedn = gRedp") 最後に，モデルが変わったので，データセットの列がモデルで指定した判断カテゴリの順番に対応するように変更する。この場合，
赤色で呈示されたポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
→青色で呈示されたポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
→非呈示のポジティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
→赤色で呈示されたネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
→青色で呈示されたネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
→非呈示のネガティブな形容詞を赤と判断する割合，青と判断する割合，呈示されていないと判断する割合
の順に18列ならぶことになる。
データは適当に入力したものを用いた，繰り返しになるが，仮説の根拠は何もない。
実際に検定
まず，MPTinRでパラメータを推定してもらう。
library("MPTinR")
dat2 <- read.csv("Imagine data2.csv")
r1 <- list("dRed = dBlue")
r2 <- list("dRed = dBlue", "gRedn = gRedp")
results1 <- fit.mpt(dat2, "model2.txt", r1)
results2 <- fit.mpt(dat2, "model2.txt", r2)とりあえず，モデル適合度と，パラメータ推定値を見てみる。
results1[["goodness.of.fit"]][["aggregated"]]
results1[["parameters"]][["aggregated"]]
results2[["goodness.of.fit"]][["aggregated"]]
results2[["parameters"]][["aggregated"]]今回は，
> results1[["goodness.of.fit"]][["aggregated"]]
Log.Likelihood G.Squared df p.value
1 -40.4139 0.431527 7 0.9996606
> results2[["goodness.of.fit"]][["aggregated"]]
Log.Likelihood G.Squared df p.value
1 -40.48077 0.5652571 8 0.9997877
> results1[["goodness.of.fit"]][["aggregated"]]
Log.Likelihood G.Squared df p.value
1 -40.4139 0.431527 7 0.9996606
> results1[["parameters"]][["aggregated"]]
estimates lower.conf upper.conf restricted.parameter
b 0.7186289 0.45713745 0.9801204
dBlue 0.5952360 0.35670172 0.8337703
dNot 0.6366538 0.32486910 0.9484384
dRed 0.5952360 0.35670172 0.8337703 dBlue
gRedn 0.3920774 -0.01764863 0.8018034
gRedp 0.5017784 0.08469464 0.9188622
> results2[["goodness.of.fit"]][["aggregated"]]
Log.Likelihood G.Squared df p.value
1 -40.48077 0.5652571 8 0.9997877
> results2[["parameters"]][["aggregated"]]
estimates lower.conf upper.conf restricted.parameter
b 0.7178431 0.4557775 0.9799087
dBlue 0.5963632 0.3583412 0.8343853
dNot 0.6362560 0.3239450 0.9485670
dRed 0.5963632 0.3583412 0.8343853 dBlue
gRedn 0.4472290 0.1525777 0.7418802 gRedp
gRedp 0.4472290 0.1525777 0.7418802 と返って来た。仮説に関わる場所だけ抜き出してくると，gRedpがgRednと異なるとした場合，その値はそれぞれgRedp = 0.502，gRedn = 0.392となっており，値だけ見れば仮説通りの結果となっている。繰り返しになるが，仮説の根拠は何もない。
検定を行うには，モデル適合度の比較を行うと述べた。まず，前述の select.mpt関数で，どちらの制限がモデルの適合が良いのか確認してみる。
> select.mpt(list(r = results1, rr = results2), output = "standard")
model n.parameters G.Squared.sum df.sum p.sum p.smaller.05 FIA.penalty.sum delta.FIA.sum FIA.best
1 r 5 6.641723 63 1 0 8.189363 0.000000 9
2 rr 4 7.589235 72 1 0 9.267269 1.551662 0
delta.AIC.sum wAIC.sum AIC.best delta.BIC.sum wBIC.sum BIC.best
1 17.05249 0.000198 0 34.95334 0 0
2 0.00000 0.999802 9 0.00000 1 9この結果を見ると，rrモデルの適合がより良い。つまりgRedpはgRednと等しいという帰無仮説が採択されたと言える。また，前述したAIC・BICとFIAの結果の逆転が起きている。
パラメータ推定値を参加者間で比較
パラメータ推定値の比較1ではパラメータ推定値の比較を参加者内行ったが，参加者間でも行うことが可能である。例えば，想定実験の結果を大学生と高齢者，幼児と成人で比較するといった研究を想定する。
わざわざ項目立てしたが，なんということはなく参加者内と同じように扱える。MPTinRで用いるデータセットは列がモデルで指定した判断カテゴリの順番に対応しているが，行は特に何かに対応しているわけではない（ここまででは行が参加者に対応していた）。データセットは方程式の解を得るための値を与えているのに過ぎないので，それぞれの参加者条件に対応するパラメータだけ違うモデルを複数つくり，列にはモデルで指定した判断カテゴリの順番に対応するようデータを並べればよい。同じ行に複数人の参加者のデータが入るが問題ない。
この記事が気に入ったらサポートをしてみませんか？