R Advent Calendarに参加しちゃったので,なんか考えないと行かんなあ,ということで今までやってきたことをぐるりと振り返ると,やっぱり心理学的調査の業界で生きてきたので,調査データの因子分析というのが(俺の中の)王道なわけです。なので,思うところを書いてみます。
枕が長くなりますので,Rを使ったお話だけ分かればいいやという人は,こちらのサイト(Rpubsコード)に飛んでください。
さて。
私が学生だった15年ほど前は,まだ大学でもメインフレームにTSSで接続するような時代で,調査実習もSPSS/PCをつかうという頃。一回の因子分析に20分かかったので,因子数の決定にも慎重に,ドキドキしながらやったものです。
そういう「ちょっと古い」時代をしっている人間は,データは5〜7件法でとって,それを因子分析するときは主因子解で,Varimax回転するというのが王道だと習いました。斜交回転や最尤法といった技は,97-98年頃,SPSSのバージョンが8?9?にあがった頃にオプションとしてついたのであって,それまでは理論はあるけど計算機がついてこないという時代だったのです。
大学を卒業する頃には最尤法,promax回転できまり!という時代が来たわけだけど,「データは5件法で」というあたりは特に問題視されることもなくスルーされていました。統計学者に言わせると,7件法でギリギリ,9か11件法ぐらいでないと間隔尺度水準と見なすことができないとなることは知られていたんだけど,「だってしょうがないじゃない」という感じで普通に因子分析をしておったわけです。間隔尺度水準がアヤシイとされる3,4件法になるんなら,もういっそ数量化三類(=双対尺度法。名義尺度水準の分析法ですな)にしなさいと言われたり。
ところが10年ほど前から,この「似非間隔尺度水準,実質は順序尺度水準じゃないか」というデータに対して,統計業界から項目反応理論(IRT) を使えばいいじゃない,という話が出てきた。IRTの中でも段階反応モデル(GRM)は,反応が段階でとれるものに対するモデルであり,順序尺度水準のままで分析できるときたもんだ。しかもGRMは数学的にはカテゴリカル因子分析と同じ,すなわち因子分析の姿を変えただけのものであり,ポイントは普通の因子分析を始める時につかう積率相関係数のかわりに,ポリコリック相関係数を使えばいいだけ,という。まあ理解もしやすいわけですね。
じゃあそのポリコリック相関係数が算出できるソフトは,といえば,やはりR。IRTを扱う専門ソフトはBILOG-MGとかあったし,GRM(PCM)を扱うソフトはParscaleというのがあったんだけど,海外ソフトを輸入して買うという,お金がないとやってられないじだいだったわけです。
Rもまだまだヨチヨチ歩きのころで,周りに十分なテキストや解説サイトがなかったので,私も敷居が高く感じていました。
お金がないから仕方がない,とRに手を出したのが,私とRとのそもそもの出会いでもあったわけで。
さて,Rのパッケージltmにpolycor関数が入っている!これで勝つる!と思われたけど,ひとつ問題が。それは,IRTはテスト理論を背景にしていることから,基本的に一因子モデルなんですね。これは困る。心理屋さんは基本的に多因子モデルが好きなんです。もちろん,単純構造を目指すという原則があるから,因子数がわかれば因子ごとにIRT(GRM)をやって因子得点の算出に向かえばいいんだけど,やはり多因子でないとねえ・・・。
ところがどっこい。最近,mirtパッケージがこの問題を解決してくれたのです。なんと因子間相関までみとめて頂ける!モデルも結構自由に書ける!
いやー,Rの発展,展開は最近目覚ましいものがありますね。これで今のところ,変に新しいソフトを買ってその使い方を習熟して,という苦労することなく(Parscaleをdisっているわけではない。パー助はパー助で可愛かったんです。),いつものRスクリプトでやりたいことが全部できるようになっちゃった。
ということで,話が長くなりましたが,これからの因子分析は
1.5件法ぐらいであれば,ポリコリック相関係数をもとに並行分析等で因子数を決定し,
2.多次元段階反応モデルwith最尤推定&プロマックス回転で尺度水準,抽出法,回転法もばっちり!
という方向性に進んでいくと思われます。
「5件法をまんま主因子法varimax回転」が許されるのは小学生までだよねー!
という時代が来るかどうか分かりませんが・・・ウヒヒ。
さあでは,実際にRでやってみましょう。以下ではコードと出力の一部を書いていきますが,結果を伴うRpubsの方にもリンクを貼っておきます。
それでは4件法の順序尺度水準で得られたデータを因子分析する例についてお話します。
従来通りの因子分析の方法,段階反応モデルを使った方法,多次元IRTを使う方法,の三段階にわけて比較検討しながら進めてまいります。
今回は次のパッケージをご用意ください。
> library<span class="synSpecial">(</span>psych<span class="synSpecial">)</span> > library<span class="synSpecial">(</span>ltm<span class="synSpecial">)</span> > library<span class="synSpecial">(</span>mirt<span class="synSpecial">)</span>
サンプルデータはltmパッケージにあるScienceデータを使います。4件法で7つの項目があります。
> data<span class="synSpecial">(</span>Science<span class="synSpecial">)</span> > summary<span class="synSpecial">(</span>Science<span class="synSpecial">)</span>
従来はこうした4件法であっても,(無理矢理)間隔尺度水準とみなして分析していたわけです。
すなわち,相関係数の出し方がピアソンの積率相関係数で,それに基づく因子分析だった。
ちょっとやってみます。まず間隔尺度水準に置き換えます。
Science.n <span class="synStatement"><-</span> <span class="synType">data.frame</span><span class="synSpecial">(</span>data.matrix<span class="synSpecial">(</span>Science<span class="synSpecial">))</span>
積率相関係数とポリコリック相関係数を比較します。ポリコリック相関係数は,ltmパッケージが読み込むpolycorパッケージにあるhetcor関数で,変数の型にあった適切な相関係数を出してくれる関数。
> cor<span class="synSpecial">(</span>Science.n<span class="synSpecial">)</span> > result.het <span class="synStatement"><-</span> hetcor<span class="synSpecial">(</span>Science<span class="synSpecial">)</span> > result.het$correlations
これが結果。ピアソンの積率相関係数。
Com Env Wor Fut Tec Ind Ben Com <span class="synConstant">1.000</span> <span class="synConstant">0.077</span> <span class="synConstant">0.151</span> <span class="synConstant">0.284</span> <span class="synConstant">0.071</span> <span class="synConstant">0.135</span> <span class="synConstant">0.334</span> Env <span class="synConstant">0.077</span> <span class="synConstant">1.000</span> -<span class="synConstant">0.071</span> -<span class="synConstant">0.030</span> <span class="synConstant">0.387</span> <span class="synConstant">0.328</span> -<span class="synConstant">0.032</span> Wor <span class="synConstant">0.151</span> -<span class="synConstant">0.071</span> <span class="synConstant">1.000</span> <span class="synConstant">0.401</span> -<span class="synConstant">0.089</span> -<span class="synConstant">0.019</span> <span class="synConstant">0.167</span> Fut <span class="synConstant">0.284</span> -<span class="synConstant">0.030</span> <span class="synConstant">0.401</span> <span class="synConstant">1.000</span> -<span class="synConstant">0.028</span> <span class="synConstant">0.062</span> <span class="synConstant">0.313</span> Tec <span class="synConstant">0.071</span> <span class="synConstant">0.387</span> -<span class="synConstant">0.089</span> -<span class="synConstant">0.028</span> <span class="synConstant">1.000</span> <span class="synConstant">0.345</span> -<span class="synConstant">0.012</span> Ind <span class="synConstant">0.135</span> <span class="synConstant">0.328</span> -<span class="synConstant">0.019</span> <span class="synConstant">0.062</span> <span class="synConstant">0.345</span> <span class="synConstant">1.000</span> <span class="synConstant">0.086</span> Ben <span class="synConstant">0.334</span> -<span class="synConstant">0.032</span> <span class="synConstant">0.167</span> <span class="synConstant">0.313</span> -<span class="synConstant">0.012</span> <span class="synConstant">0.086</span> <span class="synConstant">1.000</span>
ポリコリック相関係数はこちら。
Com Env Wor Fut Tec Ind Ben Com <span class="synConstant">1.000</span> <span class="synConstant">0.099</span> <span class="synConstant">0.2012</span> <span class="synConstant">0.346</span> <span class="synConstant">0.090</span> <span class="synConstant">0.1857</span> <span class="synConstant">0.408</span> Env <span class="synConstant">0.099</span> <span class="synConstant">1.000</span> -<span class="synConstant">0.0828</span> -<span class="synConstant">0.028</span> <span class="synConstant">0.464</span> <span class="synConstant">0.4108</span> -<span class="synConstant">0.037</span> Wor <span class="synConstant">0.201</span> -<span class="synConstant">0.083</span> <span class="synConstant">1.0000</span> <span class="synConstant">0.479</span> -<span class="synConstant">0.104</span> -<span class="synConstant">0.0076</span> <span class="synConstant">0.209</span> Fut <span class="synConstant">0.346</span> -<span class="synConstant">0.028</span> <span class="synConstant">0.4786</span> <span class="synConstant">1.000</span> -<span class="synConstant">0.036</span> <span class="synConstant">0.1027</span> <span class="synConstant">0.377</span> Tec <span class="synConstant">0.090</span> <span class="synConstant">0.464</span> -<span class="synConstant">0.1039</span> -<span class="synConstant">0.036</span> <span class="synConstant">1.000</span> <span class="synConstant">0.4348</span> -<span class="synConstant">0.014</span> Ind <span class="synConstant">0.186</span> <span class="synConstant">0.411</span> -<span class="synConstant">0.0076</span> <span class="synConstant">0.103</span> <span class="synConstant">0.435</span> <span class="synConstant">1.0000</span> <span class="synConstant">0.118</span> Ben <span class="synConstant">0.408</span> -<span class="synConstant">0.037</span> <span class="synConstant">0.2086</span> <span class="synConstant">0.377</span> -<span class="synConstant">0.014</span> <span class="synConstant">0.1180</span> <span class="synConstant">1.000</span>
結果をみると,ポリコリック相関係数のほうが数値が大きい。逆に言うと,順序尺度水準の変数を無理矢理間隔尺度水準とみなして積率相関係数を出すことは,不適切なモデル適用によって値が過小評価されちゃうことでもあるわけです。
で,普通の因子分析というのは,連続変数とみなしたデータに対して固有値分解し,因子数を決めたりするわけですね。
因子数の決定には,psychパッケージのfa.parallel関数による並行分析がいいかも。
> fa.parallel<span class="synSpecial">(</span>Science.n<span class="synSpecial">)</span> Parallel analysis suggests that the number of factors = <span class="synConstant">3</span> and the number of components = <span class="synConstant">2</span> > fa.parallel<span class="synSpecial">(</span>result.het$correlations<span class="synSpecial">,</span>n.obs=<span class="synConstant">392</span><span class="synSpecial">)</span> Parallel analysis suggests that the number of factors = <span class="synConstant">3</span> and the number of components = <span class="synConstant">2</span>
まあ結果は変わらないんですけどね。
因子分析も連続変数と見なした古典的方法と,ポリコリック相関係数をつかった方法,両方で見てみましょう。
まずは古典的方法から。
> fa.result.n <span class="synStatement"><-</span> fa<span class="synSpecial">(</span>Science.n<span class="synSpecial">,</span>nfactors=<span class="synConstant">3</span><span class="synSpecial">,</span>fm=<span class="synConstant">"ml"</span><span class="synSpecial">,</span>rotate=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> > print<span class="synSpecial">(</span>fa.result.n<span class="synSpecial">,</span>sort=T<span class="synSpecial">,</span>digit=<span class="synConstant">3</span><span class="synSpecial">)</span> Factor Analysis using method = ml Call<span class="synSpecial">:</span> fa<span class="synSpecial">(</span>r = Science.n<span class="synSpecial">,</span> nfactors = <span class="synConstant">3</span><span class="synSpecial">,</span> rotate = <span class="synConstant">"promax"</span><span class="synSpecial">,</span> fm = <span class="synConstant">"ml"</span><span class="synSpecial">)</span> Standardized loadings <span class="synSpecial">(</span>pattern <span class="synType">matrix</span><span class="synSpecial">)</span> based upon correlation <span class="synType">matrix</span> item ML2 ML1 ML3 h2 u2 Technology <span class="synConstant">5</span> <span class="synConstant">0.639</span> -<span class="synConstant">0.021</span> -<span class="synConstant">0.046</span> <span class="synConstant">0.403</span> <span class="synConstant">0.597</span> Environment <span class="synConstant">2</span> <span class="synConstant">0.617</span> <span class="synConstant">0.002</span> -<span class="synConstant">0.080</span> <span class="synConstant">0.372</span> <span class="synConstant">0.628</span> Industry <span class="synConstant">6</span> <span class="synConstant">0.545</span> <span class="synConstant">0.040</span> <span class="synConstant">0.058</span> <span class="synConstant">0.314</span> <span class="synConstant">0.686</span> Future <span class="synConstant">4</span> <span class="synConstant">0.018</span> <span class="synConstant">0.793</span> -<span class="synConstant">0.010</span> <span class="synConstant">0.620</span> <span class="synConstant">0.380</span> Work <span class="synConstant">3</span> -<span class="synConstant">0.082</span> <span class="synConstant">0.548</span> -<span class="synConstant">0.067</span> <span class="synConstant">0.273</span> <span class="synConstant">0.727</span> Benefit <span class="synConstant">7</span> -<span class="synConstant">0.076</span> -<span class="synConstant">0.044</span> <span class="synConstant">0.799</span> <span class="synConstant">0.587</span> <span class="synConstant">0.413</span> Comfort <span class="synConstant">1</span> <span class="synConstant">0.119</span> <span class="synConstant">0.170</span> <span class="synConstant">0.341</span> <span class="synConstant">0.236</span> <span class="synConstant">0.764</span> ML2 ML1 ML3 SS loadings <span class="synConstant">1.103</span> <span class="synConstant">0.952</span> <span class="synConstant">0.751</span> Proportion Var <span class="synConstant">0.158</span> <span class="synConstant">0.136</span> <span class="synConstant">0.107</span> Cumulative Var <span class="synConstant">0.158</span> <span class="synConstant">0.294</span> <span class="synConstant">0.401</span> Proportion Explained <span class="synConstant">0.393</span> <span class="synConstant">0.339</span> <span class="synConstant">0.268</span> Cumulative Proportion <span class="synConstant">0.393</span> <span class="synConstant">0.732</span> <span class="synConstant">1.000</span> With factor correlations of ML2 ML1 ML3 ML2 <span class="synConstant">1.000</span> -<span class="synConstant">0.008</span> <span class="synConstant">0.154</span> ML1 -<span class="synConstant">0.008</span> <span class="synConstant">1.000</span> <span class="synConstant">0.560</span> ML3 <span class="synConstant">0.154</span> <span class="synConstant">0.560</span> <span class="synConstant">1.000</span> <span class="synSpecial">(</span>以下略<span class="synSpecial">)</span>
ポリコリック相関係数を使った方。
> fa.result <span class="synStatement"><-</span> fa<span class="synSpecial">(</span>result.het$correlations<span class="synSpecial">,</span>nfactors=<span class="synConstant">3</span><span class="synSpecial">,</span>n.obs=<span class="synConstant">392</span><span class="synSpecial">,</span>rotate=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> > print<span class="synSpecial">(</span>fa.result<span class="synSpecial">,</span>sort=T<span class="synSpecial">,</span>digit=<span class="synConstant">3</span><span class="synSpecial">)</span> Factor Analysis using method = ml Call<span class="synSpecial">:</span> fa<span class="synSpecial">(</span>r = result.het$correlations<span class="synSpecial">,</span> nfactors = <span class="synConstant">3</span><span class="synSpecial">,</span> n.obs = <span class="synConstant">392</span><span class="synSpecial">,</span> rotate = <span class="synConstant">"promax"</span><span class="synSpecial">,</span> fm = <span class="synConstant">"ml"</span><span class="synSpecial">)</span> Standardized loadings <span class="synSpecial">(</span>pattern <span class="synType">matrix</span><span class="synSpecial">)</span> based upon correlation <span class="synType">matrix</span> item ML2 ML3 ML1 h2 u2 Technology <span class="synConstant">5</span> <span class="synConstant">0.700</span> -<span class="synConstant">0.049</span> -<span class="synConstant">0.047</span> <span class="synConstant">0.485</span> <span class="synConstant">0.515</span> Environment <span class="synConstant">2</span> <span class="synConstant">0.674</span> -<span class="synConstant">0.009</span> -<span class="synConstant">0.089</span> <span class="synConstant">0.443</span> <span class="synConstant">0.557</span> Industry <span class="synConstant">6</span> <span class="synConstant">0.630</span> <span class="synConstant">0.065</span> <span class="synConstant">0.055</span> <span class="synConstant">0.421</span> <span class="synConstant">0.579</span> Future <span class="synConstant">4</span> <span class="synConstant">0.016</span> <span class="synConstant">0.832</span> -<span class="synConstant">0.009</span> <span class="synConstant">0.686</span> <span class="synConstant">0.314</span> Work <span class="synConstant">3</span> -<span class="synConstant">0.090</span> <span class="synConstant">0.633</span> -<span class="synConstant">0.090</span> <span class="synConstant">0.351</span> <span class="synConstant">0.649</span> Benefit <span class="synConstant">7</span> -<span class="synConstant">0.075</span> -<span class="synConstant">0.050</span> <span class="synConstant">0.895</span> <span class="synConstant">0.736</span> <span class="synConstant">0.264</span> Comfort <span class="synConstant">1</span> <span class="synConstant">0.143</span> <span class="synConstant">0.213</span> <span class="synConstant">0.352</span> <span class="synConstant">0.293</span> <span class="synConstant">0.707</span> ML2 ML3 ML1 SS loadings <span class="synConstant">1.364</span> <span class="synConstant">1.132</span> <span class="synConstant">0.919</span> Proportion Var <span class="synConstant">0.195</span> <span class="synConstant">0.162</span> <span class="synConstant">0.131</span> Cumulative Var <span class="synConstant">0.195</span> <span class="synConstant">0.357</span> <span class="synConstant">0.488</span> Proportion Explained <span class="synConstant">0.399</span> <span class="synConstant">0.331</span> <span class="synConstant">0.269</span> Cumulative Proportion <span class="synConstant">0.399</span> <span class="synConstant">0.731</span> <span class="synConstant">1.000</span> With factor correlations of ML2 ML3 ML1 ML2 <span class="synConstant">1.000</span> <span class="synConstant">0.025</span> <span class="synConstant">0.160</span> ML3 <span class="synConstant">0.025</span> <span class="synConstant">1.000</span> <span class="synConstant">0.573</span> ML1 <span class="synConstant">0.160</span> <span class="synConstant">0.573</span> <span class="synConstant">1.000</span> (以下略)
これまた結果は大きく変わらないのですが,後者の方がやや大きい負荷量が算出されている。前者の方がやや「目減り」していたわけです。
後者がカテゴリカル因子分析のやり方ですが,いったんhetcor関数でポリコリック相関係数を算出させるあたりがちょっといやですね。
ということで,ルートとしては,この因子分析の結果を受けてIRT(GRM)ということが考えられます。
最近,psychパッケージにirt.faというポリコリック相関係数から因子分析できる関数が追加されました。
変数は数値型(Factor型ではない)で渡す必要がありますが,次のようにします。
> result.irtFa <span class="synStatement"><-</span> irt.fa<span class="synSpecial">(</span>Science.n<span class="synSpecial">,</span>nfactors=<span class="synConstant">3</span><span class="synSpecial">)</span> > result.irtFa$fa Factor Analysis using method = ml Call<span class="synSpecial">:</span> fa<span class="synSpecial">(</span>r = r<span class="synSpecial">,</span> nfactors = nfactors<span class="synSpecial">,</span> n.obs = n.obs<span class="synSpecial">,</span> fm = <span class="synConstant">"ml"</span><span class="synSpecial">)</span> Standardized loadings <span class="synSpecial">(</span>pattern <span class="synType">matrix</span><span class="synSpecial">)</span> based upon correlation <span class="synType">matrix</span> ML2 ML3 ML1 h2 u2 Comfort <span class="synConstant">0.17</span> <span class="synConstant">0.24</span> <span class="synConstant">0.33</span> <span class="synConstant">0.28</span> <span class="synConstant">0.72</span> Environment <span class="synConstant">0.67</span> -<span class="synConstant">0.02</span> -<span class="synConstant">0.06</span> <span class="synConstant">0.44</span> <span class="synConstant">0.56</span> Work -<span class="synConstant">0.09</span> <span class="synConstant">0.61</span> -<span class="synConstant">0.06</span> <span class="synConstant">0.35</span> <span class="synConstant">0.65</span> Future <span class="synConstant">0.02</span> <span class="synConstant">0.82</span> <span class="synConstant">0.02</span> <span class="synConstant">0.70</span> <span class="synConstant">0.30</span> Technology <span class="synConstant">0.69</span> -<span class="synConstant">0.05</span> -<span class="synConstant">0.02</span> <span class="synConstant">0.49</span> <span class="synConstant">0.51</span> Industry <span class="synConstant">0.63</span> <span class="synConstant">0.07</span> <span class="synConstant">0.07</span> <span class="synConstant">0.42</span> <span class="synConstant">0.58</span> Benefit -<span class="synConstant">0.02</span> -<span class="synConstant">0.01</span> <span class="synConstant">0.90</span> <span class="synConstant">0.80</span> <span class="synConstant">0.20</span> ML2 ML3 ML1 SS loadings <span class="synConstant">1.37</span> <span class="synConstant">1.14</span> <span class="synConstant">0.96</span> Proportion Var <span class="synConstant">0.20</span> <span class="synConstant">0.16</span> <span class="synConstant">0.14</span> Cumulative Var <span class="synConstant">0.20</span> <span class="synConstant">0.36</span> <span class="synConstant">0.50</span> Proportion Explained <span class="synConstant">0.39</span> <span class="synConstant">0.33</span> <span class="synConstant">0.28</span> Cumulative Proportion <span class="synConstant">0.39</span> <span class="synConstant">0.72</span> <span class="synConstant">1.00</span> With factor correlations of ML2 ML3 ML1 ML2 <span class="synConstant">1.00</span> <span class="synConstant">0.00</span> <span class="synConstant">0.06</span> ML3 <span class="synConstant">0.00</span> <span class="synConstant">1.00</span> <span class="synConstant">0.49</span> ML1 <span class="synConstant">0.06</span> <span class="synConstant">0.49</span> <span class="synConstant">1.00</span> <span class="synSpecial">(</span>以下略)
結果の全体像をみるには,print関数でshort=Fオプションをつけましょう。
> print<span class="synSpecial">(</span>result.irtFa<span class="synSpecial">,</span>shor=F<span class="synSpecial">)</span> (長いので略)
まあ他にも色々できるようですが,まだちょっと関数の挙動が怪しかったりします。
IRT(GRM)をやるには,以前からある専門のltmパッケージ使った方が確実かも。
今回第一因子はTechnology,Environment,Industryからなり,第二因子はFuture,Work,第三因子はBenefit,Comfortからなるわけですから,ltmパッケージのgrm関数を使って,
> grm<span class="synSpecial">(</span>Science<span class="synSpecial">[</span>c<span class="synSpecial">(</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">)])</span> <span class="synComment">#Factor 1</span> > grm<span class="synSpecial">(</span>Science<span class="synSpecial">[</span>c<span class="synSpecial">(</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">)])</span> <span class="synComment">#Factor 2</span> > grm<span class="synSpecial">(</span>Science<span class="synSpecial">[</span>c<span class="synSpecial">(</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">)])</span> <span class="synComment">#Factor 3</span> Call<span class="synSpecial">:</span> grm<span class="synSpecial">(</span>data = Science<span class="synSpecial">[</span>c<span class="synSpecial">(</span><span class="synConstant">5</span><span class="synSpecial">,</span> <span class="synConstant">2</span><span class="synSpecial">,</span> <span class="synConstant">6</span><span class="synSpecial">)])</span> Coefficients<span class="synSpecial">:</span> Extrmt1 Extrmt2 Extrmt3 Dscrmn Technology -<span class="synConstant">2.386</span> -<span class="synConstant">0.855</span> <span class="synConstant">0.624</span> <span class="synConstant">1.751</span> Environment -<span class="synConstant">2.112</span> -<span class="synConstant">0.772</span> <span class="synConstant">0.618</span> <span class="synConstant">1.626</span> Industry -<span class="synConstant">3.005</span> -<span class="synConstant">1.589</span> <span class="synConstant">0.301</span> <span class="synConstant">1.522</span> Log.Lik<span class="synSpecial">:</span> -<span class="synConstant">1312.057</span> <span class="synSpecial">(</span>以下略)
とすることができます。
もっとも,これだと因子間相関がでないですね。IRT(GRM)は基本的に単因子モデルですから。
これを多次元モデルにするのがmirtパッケージのmirt関数です。数値型の変数を渡して,
> result.mirt <span class="synStatement"><-</span> mirt<span class="synSpecial">(</span>Science.n<span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">,</span>itemtype=<span class="synConstant">"graded"</span><span class="synSpecial">,</span>rotate=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> > summary<span class="synSpecial">(</span>result.mirt<span class="synSpecial">)</span> Rotation<span class="synSpecial">:</span> promax Rotated factor loadings<span class="synSpecial">:</span> F_1 F_2 F_3 h2 Comfort <span class="synConstant">0.293</span> -<span class="synConstant">0.218</span> <span class="synConstant">0.347</span> <span class="synConstant">0.421</span> Environment -<span class="synConstant">0.113</span> -<span class="synConstant">0.665</span> -<span class="synConstant">0.039</span> <span class="synConstant">0.429</span> Work -<span class="synConstant">0.148</span> <span class="synConstant">0.048</span> <span class="synConstant">0.705</span> <span class="synConstant">0.390</span> Future <span class="synConstant">0.049</span> <span class="synConstant">0.038</span> <span class="synConstant">0.747</span> <span class="synConstant">0.591</span> Technology <span class="synConstant">0.077</span> -<span class="synConstant">0.679</span> -<span class="synConstant">0.141</span> <span class="synConstant">0.449</span> Industry -<span class="synConstant">0.069</span> -<span class="synConstant">0.680</span> <span class="synConstant">0.107</span> <span class="synConstant">0.487</span> Benefit <span class="synConstant">1.000</span> <span class="synConstant">0.068</span> <span class="synConstant">0.011</span> <span class="synConstant">0.996</span> Rotated SS loadings<span class="synSpecial">:</span> <span class="synConstant">1.133</span> <span class="synConstant">1.422</span> <span class="synConstant">1.207</span> Factor correlations<span class="synSpecial">:</span> F_1 F_2 F_3 F_1 <span class="synConstant">1.000</span> -<span class="synConstant">0.145</span> <span class="synConstant">0.569</span> F_2 -<span class="synConstant">0.145</span> <span class="synConstant">1.000</span> -<span class="synConstant">0.212</span> F_3 <span class="synConstant">0.569</span> -<span class="synConstant">0.212</span> <span class="synConstant">1.000</span>
となります(警告が出るのでちょっと推定がうまくいってないかもです)
ちなみにmirt関数は確認的にモデルを書くこともできます。
confmirt.model()関数で3つの因子と対応する項目番号から,次のようにモデルを描くことができます。
> cmodel <span class="synStatement"><-</span> confmirt.model<span class="synSpecial">()</span> > F1=<span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">6</span> > F2=<span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">3</span> > F3=<span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">1</span> > > result.cmirt <span class="synStatement"><-</span> mirt<span class="synSpecial">(</span>Science<span class="synSpecial">,</span>cmodel<span class="synSpecial">,</span>itemtype=<span class="synConstant">"graded"</span><span class="synSpecial">)</span> > summary<span class="synSpecial">(</span>result.cmirt<span class="synSpecial">)</span> > Factor loadings metric<span class="synSpecial">:</span> > F1 F2 F3 h2 > Comfort <span class="synConstant">0.000</span> <span class="synConstant">0.000</span> <span class="synConstant">0.311</span> <span class="synConstant">0.0969</span> > Environment <span class="synConstant">0.609</span> <span class="synConstant">0.000</span> <span class="synConstant">0.000</span> <span class="synConstant">0.3713</span> > Work <span class="synConstant">0.000</span> <span class="synConstant">0.433</span> <span class="synConstant">0.000</span> <span class="synConstant">0.1874</span> > Future <span class="synConstant">0.000</span> <span class="synConstant">0.445</span> <span class="synConstant">0.000</span> <span class="synConstant">0.1980</span> > Technology <span class="synConstant">0.587</span> <span class="synConstant">0.000</span> <span class="synConstant">0.000</span> <span class="synConstant">0.3447</span> > Industry <span class="synConstant">0.655</span> <span class="synConstant">0.000</span> <span class="synConstant">0.000</span> <span class="synConstant">0.4293</span> > Benefit <span class="synConstant">0.000</span> <span class="synConstant">0.000</span> <span class="synConstant">0.357</span> <span class="synConstant">0.1276</span> > > SS loadings<span class="synSpecial">:</span> <span class="synConstant">1.145</span> <span class="synConstant">0.385</span> <span class="synConstant">0.224</span> > > Factor correlations<span class="synSpecial">:</span> > F1 F2 F3 > F1 <span class="synConstant">1</span> <span class="synConstant">0</span> <span class="synConstant">0</span> > F2 <span class="synConstant">0</span> <span class="synConstant">1</span> <span class="synConstant">0</span> > F3 <span class="synConstant">0</span> <span class="synConstant">0</span> <span class="synConstant">1</span>
とまあ,このようにかなりカテゴリカルな因子分析ができる環境が整ってきているわけです。
ltmは安定した答えを出してくれますが,psychパッケージやmirtパッケージはまだちょっと不十分かな,とおもうところがなくはありません。
その辺は今後に期待ということで。
余談ですが,もっとちゃんとしたいという人は,M-plusというソフトがありますよ。
補遺)時代背景については,あくまでも私の個人史からのお話ですので,ご容赦ください。ちなみに私は1994年にKUに入学,1998年〜2003年はKGにおりました。