Rによる因子分析法、いくつかあるので比較してみた。遠く離れた学生に対する補講の意味と、自分に対する備忘録をかねて。
比較した方法は、
- psychパッケージによる主因子解
- factanal関数による最尤解
- MCMCパッケージによるMCMC的（ベイズ的）因子分析法
である。
取りあえず、ソースは以下の通り。これをコピペすると、再現できるはず。
ちなみに、こちらの実行環境はR2.9.0である。
library<span class="synSpecial">(</span>psych<span class="synSpecial">)</span> library<span class="synSpecial">(</span>MCMCpack<span class="synSpecial">)</span> <span class="synComment">#サンプルデータ/psychパックから。</span> data<span class="synSpecial">(</span>bfi<span class="synSpecial">)</span> <span class="synComment">#平行分析による因子数の決定</span> fa.parallel<span class="synSpecial">(</span>bfi<span class="synSpecial">)</span> <span class="synComment">#主因子解</span> fpa.out <span class="synStatement"><-</span> factor.pa<span class="synSpecial">(</span>bfi<span class="synSpecial">,</span>nfactors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotate=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> print<span class="synSpecial">(</span>fpa.out<span class="synSpecial">,</span>cut=<span class="synConstant">0</span><span class="synSpecial">)</span> <span class="synComment">#最尤解</span> fml.out <span class="synStatement"><-</span> factanal<span class="synSpecial">(</span>~.<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotation=<span class="synConstant">"none"</span><span class="synSpecial">)</span> pro <span class="synStatement"><-</span> promax<span class="synSpecial">(</span>loadings<span class="synSpecial">(</span>fml.out<span class="synSpecial">),</span>m=<span class="synConstant">3</span><span class="synSpecial">)</span> solve<span class="synSpecial">(</span>t<span class="synSpecial">(</span>pro$rotmat<span class="synSpecial">)</span>%*%pro$rotmat<span class="synSpecial">)</span> <span class="synComment">#MCMCpack</span> mcmc.out<span class="synStatement"><-</span> MCMCfactanal<span class="synSpecial">(</span>~.<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>burnin=<span class="synConstant">10000</span><span class="synSpecial">,</span>mcmc=<span class="synConstant">50000</span><span class="synSpecial">)</span> summary<span class="synSpecial">(</span>mcmc.out<span class="synSpecial">)</span> plot<span class="synSpecial">(</span>mcmc.out<span class="synSpecial">)</span>
解説。まずはライブラリの読み込み。
library<span class="synSpecial">(</span>psych<span class="synSpecial">)</span> library<span class="synSpecial">(</span>MCMCpack<span class="synSpecial">)</span>
psychパッケージは心理系のデータ解析に必要なものが、MCMCパッケージにはマルコフ連鎖モンテカルロ法による統計に必要なものが、一揃えはいってます。
持っていない人は、次のコードを先に実行して、ライブラリをインストールしておいてね。
install.packages<span class="synSpecial">(</span><span class="synConstant">"psych"</span><span class="synSpecial">)</span> install.packages<span class="synSpecial">(</span><span class="synConstant">"MCMCpack"</span><span class="synSpecial">)</span>
psychパッケージに含まれるサンプルデータを使います。
data<span class="synSpecial">(</span>bfi<span class="synSpecial">)</span>
続いて因子数を決定する作業。psychパッケージにある、平行分析関数が便利。
<span class="synComment">#平行分析による因子数の決定</span> fa.parallel<span class="synSpecial">(</span>bfi<span class="synSpecial">)</span>
これを実行すると、こうなる。
△の印が因子分析のスクリープロット。これを見て因子数を決めることができる。ちなみに×は主成分解、つまり共通性の推定をしないときのスクリープロット。
点線(二種)が平行分析、要するに同じデータセットの乱数表から作られる相関行列や、データセットの数値を適当にリサンプリングして作られる相関行列のスクリープロットである。乱数やバラバラになったデータセットから得られる固有値の並びに「意味」はなく、当然そのラインより上に来る因子が有意味、と判断される。今回は3因子を仮定しました。
で、本番。
psychパッケージに含まれる主因子解の関数はfactor.pa()。引数として、データ、因子の数、回転方法、因子得点の算出等々が選べる。回転は”none”, “varimax”, “promax” “oblimin”から選択。今回は因子数とプロマックス回転だけ指定して実行。詳しくはヘルプを参照。
<span class="synComment">#主因子解</span> fpa.out <span class="synStatement"><-</span> factor.pa<span class="synSpecial">(</span>bfi<span class="synSpecial">,</span>nfactors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotate=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> print<span class="synSpecial">(</span>fpa.out<span class="synSpecial">,</span>cut=<span class="synConstant">0</span><span class="synSpecial">)</span>
次はデフォルトで入っているfactanal関数。次の箇所がそれ。
<span class="synComment">#最尤解</span> fml.out <span class="synStatement"><-</span> factanal<span class="synSpecial">(</span>~.<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotation=<span class="synConstant">"none"</span><span class="synSpecial">)</span>
factanal、と書いた次の引数だけど、「~.」は与えられたデータセットの中に含まれる変数を全部分析対象にして、という意味。データセットは次の「data=bfi」で指定している。例えばこの中で、変数A1〜A5だけを使いたいのであれば、
factanal<span class="synSpecial">(</span>~A1+A2+A3+A4+A5<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotation=<span class="synConstant">"none"</span><span class="synSpecial">)</span>
のようにする。これはRの定まった書き方ですな。
さて、続きの引数で因子数(今度はfactorsと書く。factor.paではnfactorsでしたが、nが要らないことに注意)や、回転ナシを指定している。なぜ回転しないかというと、回転用の関数promaxを使うからです。なぜわざわざそうするかというと、factanalでpromax回転をしてくれるのだけど、因子間相関が出ないからです。因子間相関を出すための式が次。
pro <span class="synStatement"><-</span> promax<span class="synSpecial">(</span>loadings<span class="synSpecial">(</span>fml.out<span class="synSpecial">),</span>m=<span class="synConstant">3</span><span class="synSpecial">)</span> solve<span class="synSpecial">(</span>t<span class="synSpecial">(</span>pro$rotmat<span class="synSpecial">)</span>%*%pro$rotmat<span class="synSpecial">)</span>
ちなみに、出力が綺麗なfactanal2関数というのが青木先生のサイトからダウンロードできる。このサイトはとても便利です。この関数を借りると、こんな感じになる。
source<span class="synSpecial">(</span><span class="synConstant">"http://aoki2.si.gunma-u.ac.jp/R/src/factanal2.R"</span><span class="synSpecial">,</span> encoding=<span class="synConstant">"<a class="keyword" href="http://d.hatena.ne.jp/keyword/euc">euc</a>-jp"</span><span class="synSpecial">)</span> factanal2<span class="synSpecial">(</span>bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>rotation=<span class="synConstant">"promax"</span><span class="synSpecial">)</span> source<span class="synSpecial">(</span><span class="synConstant">"http://aoki2.si.gunma-u.ac.jp/R/src/factor.correlation.R"</span><span class="synSpecial">,</span> encoding=<span class="synConstant">"<a class="keyword" href="http://d.hatena.ne.jp/keyword/euc">euc</a>-jp"</span><span class="synSpecial">)</span> factor.correlation<span class="synSpecial">(</span>~.<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">)</span>
source〜とあるのがネットから関数を取ってきなさいよ、という命令文。ということで、ネット環境があるときにこれを使う。なければソースを丸ごとぺたりと貼ればよい。これの二行目、四行目で上と同じことをしているが、出力が良いのでこちらを好まれる向きもあるかと。
さて、最後、今回の目玉（？）マルコフ連鎖モンテカルロ法(以下MCMC)による因子分析モデルの実行。MCMCの実行については、BRugsやR2WinBugsなど、OpenBugsを援用するのが一つのやり方。もう一つのやり方が、MCMCpackの利用。とても簡単にできます。因子分析はfactanalの前にMCMCがつくだけである。次の箇所がそう。
<span class="synComment">#MCMCpack</span> mcmc.out<span class="synStatement"><-</span> MCMCfactanal<span class="synSpecial">(</span>~.<span class="synSpecial">,</span>data=bfi<span class="synSpecial">,</span>factors=<span class="synConstant">3</span><span class="synSpecial">,</span>burnin=<span class="synConstant">10000</span><span class="synSpecial">,</span>mcmc=<span class="synConstant">50000</span><span class="synSpecial">)</span> summary<span class="synSpecial">(</span>mcmc.out<span class="synSpecial">)</span> plot<span class="synSpecial">(</span>mcmc.out<span class="synSpecial">)</span>
factanal関数と違うのは、mcmc=50000とburnin=10000という引数。前者はシミュレーション回数で、後者は最初の10000回は数値が落ち着かないのでデータとして見なしませんよ、という意味。バーン・イン期間と呼ばれます。回数は何回でもいいんだけど、数千回はするべき。ちなみに今回、五万回回して前の一万回を捨てる、という設定だが、これで要した時間は約12分でした(環境：Intel Core2 Duo E6550@2.33GHz,2GのRAM,Windows VISTA)。
さて、結果の比較。因子負荷量を表にまとめてみました。
| pfa1 | pfa2 | pfa3 | ml1 | ml2 | ml3 | mcmc1 | mcmc2 | mcmc3 | |
| A1 | -0.09 | 0.13 | -0.13 | 0.244 | -0.005 | -0.017 | 0.02 | 0.09 | 0.01 |
| A2 | 0.46 | 0.00 | 0.09 | -0.366 | 0.304 | -0.093 | 0.08 | -0.09 | -0.01 |
| A3 | 0.48 | -0.02 | 0.08 | -0.377 | 0.326 | -0.083 | 0.09 | -0.09 | -0.01 |
| A4 | 0.18 | -0.16 | 0.12 | -0.299 | 0.074 | 0.027 | 0.00 | -0.09 | -0.01 |
| A5 | 0.56 | -0.10 | 0.02 | -0.507 | 0.314 | -0.147 | 0.08 | -0.14 | -0.02 |
| C1 | -0.11 | 0.04 | 0.69 | -0.197 | 0.237 | 0.600 | 0.03 | 0.00 | 0.05 |
| C2 | -0.11 | 0.09 | 0.64 | -0.138 | 0.229 | 0.543 | 0.04 | 0.02 | 0.05 |
| C3 | 0.05 | 0.03 | 0.57 | -0.266 | 0.302 | 0.476 | 0.06 | -0.02 | 0.04 |
| C4 | 0.05 | 0.26 | -0.65 | 0.446 | -0.077 | -0.507 | 0.04 | 0.12 | -0.02 |
| C5 | 0.04 | 0.31 | -0.47 | 0.404 | 0.002 | -0.384 | 0.06 | 0.13 | -0.01 |
| E1 | -0.70 | -0.05 | 0.18 | 0.365 | -0.427 | 0.332 | -0.14 | 0.08 | 0.02 |
| E2 | -0.65 | 0.15 | 0.12 | 0.529 | -0.332 | 0.247 | -0.09 | 0.15 | 0.03 |
| E3 | 0.69 | 0.11 | 0.00 | -0.394 | 0.512 | -0.157 | 0.16 | -0.06 | 0.00 |
| E4 | 0.71 | -0.11 | -0.18 | -0.507 | 0.353 | -0.339 | 0.11 | -0.14 | -0.03 |
| E5 | 0.52 | 0.12 | 0.15 | -0.336 | 0.511 | 0.032 | 0.15 | -0.03 | 0.02 |
| N1 | 0.08 | 0.79 | -0.09 | 0.647 | 0.540 | -0.016 | 0.24 | 0.32 | 0.06 |
| N2 | 0.05 | 0.72 | -0.04 | 0.605 | 0.506 | 0.031 | 0.22 | 0.31 | 0.06 |
| N3 | 0.01 | 0.71 | -0.04 | 0.553 | 0.404 | 0.001 | 0.18 | 0.27 | 0.05 |
| N4 | -0.22 | 0.59 | -0.01 | 0.598 | 0.161 | 0.062 | 0.10 | 0.25 | 0.05 |
| N5 | -0.14 | 0.48 | 0.10 | 0.431 | 0.202 | 0.120 | 0.10 | 0.20 | 0.04 |
| O1 | 0.31 | 0.17 | 0.19 | -0.171 | 0.353 | 0.088 | 0.11 | 0.00 | 0.02 |
| O2 | -0.22 | -0.10 | 0.23 | -0.024 | -0.147 | 0.213 | -0.07 | -0.02 | 0.01 |
| O3 | 0.45 | 0.20 | 0.15 | -0.242 | 0.431 | 0.032 | 0.13 | -0.01 | 0.02 |
| O4 | 0.07 | 0.29 | 0.23 | 0.065 | 0.250 | 0.155 | 0.08 | 0.07 | 0.03 |
| O5 | -0.14 | -0.04 | -0.13 | 0.124 | -0.126 | -0.076 | -0.03 | 0.02 | -0.01 |
| 因子間相関 | 1.00 | -0.18 | 0.40 | 1.000 | -0.184 | -0.079 | |||
| 1.00 | -0.07 | 1.000 | 0.333 | ||||||
| 1.00 | 1.000 |
ちょっとわかりにくいかなー。
取りあえず見て思うのは、主因子解と最尤解が結構違うと言うこと。あと、MCMCの答えがどちらかというと主因子解に近いということ。実際、各セルの差の二乗をズレと考えて指標化してみると、主因子解と最尤解はセル平均0.324のズレ。主因子解とMCMCは0.206、最尤解とMCMCは0.276のズレとなり、直感を裏付けるわけです。
三つの手法、それぞれ背景の哲学が違うから、どれに真実を見るかは人それぞれ。ユーザー視点の適当な比較ですから、ご参考まで。
ところで、二つほど注意点を。
一つは、MCMC法では回転についてふれなかった。マニュアルを見ても、特に記述がないようなので、おそらく因子間相関の仮定はあると思われる。因子得点を算出させて、その相関を見ればいいのだけど、今回そこまでしませんでした。これで因子得点まで算出させると、マシンパワーが結構いるので要注意、とマニュアルにありました。あるいは、BRugsを使ってきちんとモデル化して、因子間共分散を算出させればいいのだけど・・・色々やってみたけど、うまくモデル化できなかったんじゃよー！（T_T)自分の馬鹿さ加減にあきれます。豊田(2008)*1に因子分析モデルの例があるのだけど、確認的で単純因子構造のモデルで、うまく応用できませんでした。探索的FAのモデルの書き方、調べて教えてください。＞W
二つめ。MCMC法のメリットは、母数の分布がわかること。最後の行、
plot<span class="synSpecial">(</span>mcmc.out<span class="synSpecial">)</span>
を実行すると、次のような画面が出る。
今回表の中にのせたのはあくまでも平均であって、母数(ここでは因子負荷量)の平均、分散、パーセンタイルや、図のような分布が見える。因子負荷量と因子得点には正規分布を、独自性にはガンマ分布を事前分布として与えるのだけど、事後分布の形はまぁ・・・すごいもんです。例としては悪いサンプルでしたかね。でもま、現実のデータはもっと誤差まみれなわけで。ともかくこういう、今まではなかった特徴が見られるのがMCMCの売りなのです。
今回、学会発表のデータ分析に当たって、ちょっとMCMCで遊んでみるかと思ってやってみたら、MCMCfactanalというひろいものもあったし、結構勉強になりました。面白かったです。このデータが読者の何かの足しになれば幸いです。
何か間違っているところがあったら、どなたでも遠慮なくご指摘ください。よろしくお願いします。
*1:豊田秀樹(編著）2008 マルコフ連鎖モンテカルロ法,朝倉書店