RでSEM(共分散構造分析/構造方程式モデリング)
R を使った SEM(共分散構造分析/構造方程式モデリング)のやり方を勉強する機会があったので、メモしておきます。
環境は Windows Vista Home Premium SP2 × R 2.10.1 × Graphviz 2.26.3 です。
SEM を行うときのデータは『心理統計学の基礎:統合的理解のために』(南風原朝和)の第10章 第5節にある子どもの協調性の発達に関する例を使いました。上の画像のように、「母親が子どもの協調性に価値をおいている程度」と「幼稚園・保育園での友だちとの相互作用の経験の量」「子どもの協調性」の因果関係を考えるモデルです。それぞれの潜在変数を y01、y02、y03、y04 の4つ、y05、y06、y07、y08 の4つ、y09、y10、y11、y12 の4つで測定したデータです。
R で SEM をやる手順は、以下のとおりです。順に説明していきます。コードの書き方が R っぽくなくて慣れた方には見づらいかもしれませんが、とりあえず自分の見やすさを優先しております。悪しからず。
sem
ライブラリを読み込む
R で SEM をやるには、まずは sem
ライブラリを読み込みます。他にも SEM ができるライブラリはあるようですが、これが最も一般的なようです。
1 | library(sem) |
相関(共分散)行列を作る
続いて、相関(共分散)行列を作ります。R の sem
ライブラリでは、ローデータからではなく相関(共分散)行列を読み込んで解析を行うからです。『心理統計学の基礎』 p. 354 表10-8を以下のようにして読み込ませます。変数名はプログラマとしてどうかとは思いますが日本語で「相関行列」と付けました。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | 相関行列 <- read.moments( diag=TRUE, names=c( "y01", "y02", "y03", "y04", "y05", "y06", "y07", "y08", "y09", "y10", "y11", "y12" ) ) 1 0.160 1 0.302 0.341 1 0.461 0.400 0.372 1 0.299 0.404 0.552 0.302 1 0.152 0.320 0.476 0.225 0.708 1 0.134 0.403 0.467 0.256 0.623 0.324 1 0.182 0.374 0.572 0.255 0.776 0.769 0.724 1 0.251 0.285 0.316 0.164 0.361 0.295 0.260 0.284 1 0.372 0.100 0.408 0.236 0.294 0.206 0.071 0.142 0.295 1 0.157 0.291 0.393 0.229 0.472 0.351 0.204 0.320 0.290 0.468 1 0.206 -0.014 0.369 0.224 0.342 0.202 0.152 0.189 0.418 0.351 0.385 1 |
read.moments
は相関(共分散)行列を作る関数です。行列の対角にあたる「r = 1」の部分が含まれているなら diag
に TRUE
を、含まれていないなら FALSE
を指定します。names
部分では個々の変数名をベクトルで指定します。
もしローデータから相関行列を作るのであれば、以下のように read.csv
関数などを使ってまずデータフレームを作成し、それから相関行列を作ることになります。
1 2 | データフレーム <- read.csv("example.csv") 相関行列 <- cor(データフレーム) # 共分散行列なら cov(データフレーム) |
モデルを作る
次にモデルを作ります。R で SEM をやる際に最も引っかかりやすいところかもしれません。
『心理統計学の基礎』 p. 353 図10-6のパス図を以下のように文字列で表現します。「モデル」という変数名を付けました。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | モデル <- specify.model() 母親価値 -> y01, NA, 1 母親価値 -> y02, b12, NA 母親価値 -> y03, b13, NA 母親価値 -> y04, b14, NA 相互作用経験 -> y05, NA, 1 相互作用経験 -> y06, b22, NA 相互作用経験 -> y07, b23, NA 相互作用経験 -> y08, b24, NA 協調性 -> y09, NA, 1 協調性 -> y10, b32, NA 協調性 -> y11, b33, NA 協調性 -> y12, b34, NA 母親価値 -> 相互作用経験, g11, NA 母親価値 -> 協調性, g12, NA y01 <-> y01, e01, NA y02 <-> y02, e02, NA y03 <-> y03, e03, NA y04 <-> y04, e04, NA y05 <-> y05, e05, NA y06 <-> y06, e06, NA y07 <-> y07, e07, NA y08 <-> y08, e08, NA y09 <-> y09, e09, NA y10 <-> y10, e10, NA y11 <-> y11, e11, NA y12 <-> y12, e12, NA 母親価値 <-> 母親価値, delta01, NA 相互作用経験 <-> 相互作用経験, delta02, NA 協調性 <-> 協調性, delta03, NA |
変数 A から変数 B への因果が想定されているなら「A -> B
」、A と B の間の相関(分散)であれば「A <-> B
」、A の中の分散であれば「A <-> A
」のように表現します。
3行目の「母親価値 -> y02, b12, NA
」は、母親価値という潜在変数から y02 という観測変数への因果が想定されており、そのパスの名前が「b12」であり、係数の値が「NA
」である=値を推定するということを示しています。パス名は自分の好きな名前で OK です。
1行目はパスの名前が「NA
」で値が「1」となっていますが、これはモデルの識別のために係数の値を固定しておく必要があるのでこのようになっています。2行目、10行目も同様です。
ここのモデル記述でミスをするとあとの解析時にエラーが出ます。
解析を実行する
やっと解析実行まで来ました。このまでの道のりがけっこう長いですよねw
ここでは、「結果」という変数を作り解析結果を代入します。
1 | 結果 <- sem(モデル, 相関行列, N=50) |
sem
関数の第1引数に上で苦労して書いたモデルを、第2引数に相関(共分散)行列を、第3引数にサンプルサイズを指定します。相関(共分散)行列から解析するのでサンプルサイズをちゃんと指定してあげないといけないんですね。
結果は summary
関数で見ることができます。
1 | summary(結果) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 | Model Chisquare = 74.298 Df = 52 Pr(>Chisq) = 0.022864 Chisquare (null model) = 291.59 Df = 66 Goodness-of-fit index = 0.82725 Adjusted goodness-of-fit index = 0.74087 RMSEA index = 0.093548 90% CI: (0.036288, 0.13902) Bentler-Bonnett NFI = 0.7452 Tucker-Lewis NNFI = 0.87454 Bentler CFI = 0.90115 SRMR = 0.082692 BIC = -129.13 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -1.52000 -0.28700 0.00505 0.02800 0.32000 1.62000 Parameter Estimates Estimate Std Error z value Pr(>|z|) b11 0.42984 0.155145 2.7706 5.5959e-03 y01 <--- 母親価値 b12 0.48779 0.152358 3.2016 1.3668e-03 y02 <--- 母親価値 b13 0.79919 0.137259 5.8225 5.7977e-09 y03 <--- 母親価値 b14 0.52057 0.152250 3.4192 6.2814e-04 y04 <--- 母親価値 b22 0.93750 0.146553 6.3970 1.5848e-10 y06 <--- 相互作用経験 b23 0.85768 0.158411 5.4142 6.1548e-08 y07 <--- 相互作用経験 b24 1.13826 0.155903 7.3011 2.8555e-13 y08 <--- 相互作用経験 b32 1.15469 0.419885 2.7500 5.9591e-03 y10 <--- 協調性 b33 1.23854 0.443518 2.7925 5.2295e-03 y11 <--- 協調性 b34 1.09677 0.379612 2.8892 3.8623e-03 y12 <--- 協調性 g11 0.59725 0.147620 4.0459 5.2129e-05 相互作用経験 <--- 母親価値 g12 0.38974 0.137732 2.8297 4.6589e-03 協調性 <--- 母親価値 e01 0.81524 0.176746 4.6125 3.9789e-06 y01 <--> y01 e02 0.76207 0.169396 4.4987 6.8363e-06 y02 <--> y02 e03 0.36130 0.134066 2.6949 7.0404e-03 y03 <--> y03 e04 0.72901 0.166446 4.3799 1.1876e-05 y04 <--> y04 e05 0.29807 0.096051 3.1032 1.9142e-03 y05 <--> y05 e06 0.38307 0.087711 4.3674 1.2574e-05 y06 <--> y06 e07 0.48365 0.107430 4.5020 6.7309e-06 y07 <--> y07 e08 0.09056 0.073183 1.2374 2.1592e-01 y08 <--> y08 e09 0.70787 0.168704 4.1959 2.7176e-05 y09 <--> y09 e10 0.61050 0.157975 3.8645 1.1130e-04 y10 <--> y10 e11 0.55188 0.155441 3.5504 3.8464e-04 y11 <--> y11 e12 0.64859 0.162135 4.0003 6.3256e-05 y12 <--> y12 delta02 0.34522 0.119392 2.8915 3.8341e-03 相互作用経験 <--> 相互作用経験 delta03 0.14023 0.092364 1.5182 1.2896e-01 協調性 <--> 協調性 Iterations = 43 |
上のほうにχ2 や GFI、AGFI、RMSEA、BIC などといった適合度の指標が表示され、下のほうに推定されたパス係数の値がずらっと並びます。有意確率もここに出ています。
適合度は「モデルとデータがどれくらい一致しているか」という指標です。一般的な値の目安は以下のような感じです(ソースは授業など)。
- カイ二乗検定は有意になってはいけない
- 有意ということは「モデルとデータに統計上有意な差(違い/ズレ)がある」ということだから
- ただ、サンプルサイズが大きいと小さな差まで検出して有意になってしまうのであまりあてにならない(N が数百くらいまでなら OK)
- GFI、AGFI は1に近いほどよい
- 0.95以上が望ましい
- AGFI は自由度で調整した GFI
- RMSEA は0に近いほどよい
- 0.05以下が望ましい
- BIC はモデル間比較に使われる指標で絶対的な基準はないが小さいほうがより適合したモデル
- AIC、CAIC という指標もあるが R の
sem
ライブラリでは算出してくれない
- AIC、CAIC という指標もあるが R の
今回の例を見てみると、サンプルサイズは N = 50 しかないのにカイ二乗検定は有意になっていて、GFI は0.83、AGFI も0.74、 RMSEA も0.09と、あまり適合していません(『心理統計学の基礎』でも触れられているとおり)。適合度を上げる方法は後述します。
なお、ここで表示される係数の値は(共分散行列ではなく相関行列を使ったとしても)標準解ではありませんので、改めて求めなければいけません。標準解は std.coef
関数で求められます。
1 | std.coef(結果) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | Std. Estimate 1 b11 0.42984061 y01 <--- 母親価値 2 b12 0.48778544 y02 <--- 母親価値 3 b13 0.79918831 y03 <--- 母親価値 4 b14 0.52056707 y04 <--- 母親価値 5 0.83781311 y05 <--- 相互作用経験 6 b22 0.78544879 y06 <--- 相互作用経験 7 b23 0.71857315 y07 <--- 相互作用経験 8 b24 0.95364543 y08 <--- 相互作用経験 9 0.54048970 y09 <--- 協調性 10 b32 0.62409956 y10 <--- 協調性 11 b33 0.66942010 y11 <--- 協調性 12 b34 0.59279624 y12 <--- 協調性 13 g11 0.71286880 相互作用経験 <--- 母親価値 14 g12 0.72109403 協調性 <--- 母親価値 15 e01 0.81523705 y01 <--> y01 16 e02 0.76206537 y02 <--> y02 17 e03 0.36129805 y03 <--> y03 18 e04 0.72900993 y04 <--> y04 19 e05 0.29806919 y05 <--> y05 20 e06 0.38307019 y06 <--> y06 21 e07 0.48365263 y07 <--> y07 22 e08 0.09056039 y08 <--> y08 23 e09 0.70787089 y09 <--> y09 24 e10 0.61049974 y10 <--> y10 25 e11 0.55187673 y11 <--> y11 26 e12 0.64859262 y12 <--> y12 27 1.00000000 母親価値 <--> 母親価値 28 delta02 0.49181808 相互作用経験 <--> 相互作用経験 29 delta03 0.48002340 協調性 <--> 協調性 |
『心理統計学の基礎』 p. 353 図10-6のパス図とパス係数の値が一致してそうな雰囲気です。でもこれではやはりわかりづらいので、実際にパス図を描いてグラフィカルに見てみます。
……と、その前に、モデルの適合度を上げるための工夫もしてみましょう。mod.indices
関数を使うと、どの変数とどの変数の間にパスを引く(A matrix)あるいは相関を認める(P matrix)と適合度が上がるかを教えてくれます。
1 | mod.indices(結果) |
1 2 3 4 5 6 7 | 5 largest modification indices, A matrix: y06:y07 y07:y06 y05:協調性 y05:y11 y04:y01 8.386419 6.437148 4.924923 4.744598 4.120285 5 largest modification indices, P matrix: y07:y06 協調性:y05 y04:y01 協調性:y08 y08:y07 18.029447 5.965713 5.222873 4.332659 4.164359 |
どうやら、y07 と y06 の間に相関を認めると多少改善されるようです。さきほどのモデルの1番下に y06 <-> y07, c01, NA
という1行を加えて、またモデルの読み込み → 解析の実行とやってみます。
1 2 3 4 5 6 7 8 | モデル <- specify.model() 母親価値 -> y01, NA, 1 母親価値 -> y02, b12, NA (中略) 母親価値 <-> 母親価値, NA, 1 相互作用経験 <-> 相互作用経験, delta02, NA 協調性 <-> 協調性, delta03, NA y06 <-> y07, c01, NA |
1 2 3 4 5 6 7 8 9 10 11 12 | Model Chisquare = 45.974 Df = 51 Pr(>Chisq) = 0.67304 Chisquare (null model) = 291.59 Df = 66 Goodness-of-fit index = 0.87279 Adjusted goodness-of-fit index = 0.80544 RMSEA index = 0 90% CI: (NA, 0.075654) Bentler-Bonnett NFI = 0.84233 Tucker-Lewis NNFI = 1.0288 Bentler CFI = 1 SRMR = 0.077552 BIC = -153.54 (後略) |
見てみると、カイ二乗検定は有意 → n.s.、GFI は0.83 → 0.87、AGFI も0.74 → 0.81、 RMSEA も0.09 → 0.00、BIC も-129.13 → -153.54と適合度がよくなっています。y07 と y06 の間にパスを引くことが理論的にどうなのかはとりあえずおいておいて、このようにモデルを改善していくのも SEM の作業のうちの1つです(場合によっては1番時間がかかるところかも)。ただし、理論的に考えて「ここにパスを引く(相関を認める)のはおかしい」という場合は、いくら適合度のためとはいえパスを引いたり相関を認めたりしてはいけません。
パス図を描く
では、パス図を描いてみましょう。
図形を描く専用のソフトを使ったり Word のオートシェイプを使ったりしたほうがきれいなパス図になりやすいですが、R(とGraphviz)で描く方法も悪くないです。慣れればさくっと描けますし、パスの係数の値まで自動で書き込んでくれますから。
R で SEM のパス図を描くには、path.diagram
関数を使って書き出されたコードあるいは Dot ファイルを Graphviz という描画ソフトに読み込ませます。
path.diagram
関数は以下のように書きます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | path.diagram( 結果, file="D:/Desktop/path_diagram", output.type="dot", ignore.double=FALSE, edge.labels="values", digits=3, standardize=TRUE, node.font=c("Osaka", 20), edge.font=c("Osaka", 16), min.rank="y01, y02, y03, y04", max.rank="y05, y06, y07, y08, y09, y10, y11, y12", same.rank="相互作用経験, 協調性" ) |
path.diagram
関数はいろんな引数をとります。第1引数は SEM の解析結果です。生成された Dot ファイルまたはコードを Graphviz に(一工夫して)読み込ませればパス図のできあがりです。簡単っ!
上で使っている引数の説明は以下のとおり。help(path.diagram)
でヘルプを見れば他にも載っています。
file |
出力するファイル名を拡張子を除いて指定する。ファイルのパスの含めて指定すれば R の作業ディレクトリ以外の場所にも保存できる。 |
output.type |
出力するファイルの種類を指定する。デフォルトでは Dot ファイルと PDF ファイルの両方を出力するが、dot を指定すると Dot ファイルのみになる。PDF ファイルは基本的に日本語が化けるので Dot ファイルのみでいいと思われる。 |
ignore.double |
双方向矢印(分散あるいは相関)は表示しない。デフォルトは TRUE 。 |
edge.labels |
矢印の横に表示するものを指定する。names だとパスの名前を、values だと係数の値を、both だと両方を表示する。 |
digits |
係数の値で小数点以下何桁まで表示するか指定する。デフォルトは 2 。 |
standardize |
標準解を表示する指定。デフォルトはFALSE 。 |
node.font |
変数名を表示するフォント自体とフォントサイズをベクトルで指定する。 |
edge.font |
矢印横のパス名/係数の値を表示するフォント自体とフォントサイズをベクトルで指定する。 |
min.rank |
パス図の最左部(縦に流れるパス図なら最上部)に位置する変数を指定する。複数指定する場合はカンマで区切る。 |
max.rank |
パス図の最右部(縦に流れるパス図なら最下部)に位置する変数を指定する。複数指定する場合はカンマで区切る。 |
same.rank |
縦位置(縦に流れるパス図なら横位置)を揃えたい変数をカンマで区切って指定する。 |
上記の path.diagram
関数を実行すると、デスクトップ(デスクトップのパスは環境によって異なります)に path_diagram.dot
というファイルができていると思います。こいつを Graphviz に読み込ませれば OK! ……と言いたいところですが、どうやら文字コードの関係で(?)そのままではダメみたいですorz よって以下の2つの方法のどちらかがいいのではと思います。
- EmEditor の外部ツールから Graphviz に送る
- Graphviz の代わりに EasyGraphViz を使う
まず、「EmEditor の外部ツールから Graphviz に送る」方法の説明をします。EmEditorというのは僕がメインで使っているエディタです。EmEditor には「外部ツール」という機能があって、他のアプリケーションを呼び出すことができます。似たような機能を持つエディタであれば EmEditor でなくても OK です。
僕は以下のように設定しています。アイコンをツールバーに表示して、クリック1つでパス図が描けるようにしています。設定の際は、Graphviz – ノート – livedoor Wiki(ウィキ)を参考にさせていただきました。ありがとうございます。
- タイトル:
Graphviz (to PNG)
- コマンド:
C:\Program Files\Graphviz2.26.3\bin\dot.exe
- 引数:
-Tpng -oD:\Desktop\$(Filename).png $(FilenameEx)
-Tpdf -oD:\Desktop\$(Filename).pdf $(FilenameEx)
とすれば PDF になる(なぜか Adobe Reader では図が崩れる疑惑。Foxit Reader なら OK)
- 初期ディレクトリ:
$(Dir)
- ファイルを保存する:オン
- アウトプット バーを使用する:オン
- 終了時に閉じる:オン
- 入力:無し
- 出力:ツールチップとして表示
- エンコード:システム規定(932, Shift_JIS)
- 標準エラー:ツールチップとして表示
デスクトップのパスは環境により異なるので適宜変更する必要があります。また、文字コードは UTF-8(BOM なし)にしないとうまくいかないかもしれません。
下が、EmEditor の外部ツールを使って作成した PNG 画像です。
下は、EmEditor の外部ツールを使って作成した PDF です。
そこそこきれいなパス図じゃあないかなと思います。
次に、「Graphviz の代わりに EasyGraphViz を使う」方法の説明をします。Graphviz でうまくいかない場合、EasyGraphViz という Graphviz を簡単に扱えるようにしたソフトに読み込ませてみてください。たぶん文字化けせずにパス図が描けると思います。
下の画像が EasyGraphViz で描いたパス図です。Graphviz で描いたほうが少しだけきれいな気がします。環境のせいかもしれませんけどね。
以上です。多母集団モデルを扱えない点や GUI で直感的に操作できない点など、まだ Amos などの商用ソフトに敵わない部分がありますね。それでも sem
ライブラリと R はすばらしいと思います。とりあえず、Amos の Studen Version が Vista にも対応してくれると嬉しいと思う今日この頃です。
よろしければ以下の関連(してそうな)記事もどうぞ!
- ミニノートPCは危険? クラッキングやウイルスに注意!
- 英語学習に使えるサイト、ページ8選:リスニング、ライティングのトレーニングに
- JavaScript のお勉強でストップウォッチを作ってみた
- インターネット調査に変化の兆し:クロス・マーケティング他3社がモニターを相互利用
- リンク売買とSEO
- « 前の記事:FirefoxのEvernoteアドオンの設定画面が真っ白になって動かない不具合を解決する方法
- » 次の記事:Amazonへのリンクの横にブクログへのリンクを追加するGreasemonkeyスクリプト「Add To-Booklog Links」を書いた
きにきじNewPost:: RでSEM(共分散構造分析/構造方程式モデリング) http://www.kagitaku.com/diary…
鍵山琢実 @ライブレボリューション
» このコメントを引用してコメントする
きにきじNewPost:: RでSEM(共分散構造分析/構造方程式モデリング) http://j.mp/gXd94Z
takutakku
» このコメントを引用してコメントする
RでSEM(共分散構造分析/構造方程式モデリング) | きにきじ: SEM を行うときのデータは『心理統計学の基礎:統合的理解のために』(南風原朝和)の第10章 第5節にある子どもの協調性の発… http://bit.ly/gt23Cz http://www.seo-log.net
seo_log
» このコメントを引用してコメントする