きにきじ」:今日の気になる記事をきまぐれにご紹介

RでSEM(共分散構造分析/構造方程式モデリング)

Posted at 23:57 on January 31, 2011

Last updated at 23:25 on February 11, 2011

Category: Non-News, Note

Tags: , , , , ,


[画像] 母親価値 → {相互作用経験, 協調性}

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 っぽくなくて慣れた方には見づらいかもしれませんが、とりあえず自分の見やすさを優先しております。悪しからず。

  1. sem ライブラリを読み込む
  2. 相関(共分散)行列を作る
  3. モデルを作る
  4. 解析を実行する
  5. パス図を描く

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」の部分が含まれているなら diagTRUE を、含まれていないなら 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

上のほうにχ2GFIAGFIRMSEABIC などといった適合度の指標が表示され、下のほうに推定されたパス係数の値がずらっと並びます。有意確率もここに出ています。

適合度は「モデルとデータがどれくらい一致しているか」という指標です。一般的な値の目安は以下のような感じです(ソースは授業など)。

  • カイ二乗検定は有意になってはいけない
    • 有意ということは「モデルとデータに統計上有意な差(違い/ズレ)がある」ということだから
    • ただ、サンプルサイズが大きいと小さな差まで検出して有意になってしまうのであまりあてにならない(N が数百くらいまでなら OK)
  • GFIAGFI は1に近いほどよい
    • 0.95以上が望ましい
    • AGFI は自由度で調整した GFI
  • RMSEA は0に近いほどよい
    • 0.05以下が望ましい
  • BIC はモデル間比較に使われる指標で絶対的な基準はないが小さいほうがより適合したモデル
    • AICCAIC という指標もあるが R の sem ライブラリでは算出してくれない

今回の例を見てみると、サンプルサイズは 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つの方法のどちらかがいいのではと思います。

  1. EmEditor の外部ツールから Graphviz に送る
  2. 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 にも対応してくれると嬉しいと思う今日この頃です。

▲上に戻る▲


よろしければ以下の関連(してそうな)記事もどうぞ!


3 Responses to “RでSEM(共分散構造分析/構造方程式モデリング)”

  1. きにきじNewPost:: RでSEM(共分散構造分析/構造方程式モデリング) http://www.kagitaku.com/diary…

      

    » このコメントを引用してコメントする

  2. takutakku より:

    きにきじNewPost:: RでSEM(共分散構造分析/構造方程式モデリング) http://j.mp/gXd94Z

      

    » このコメントを引用してコメントする

  3. seo_log より:

    RでSEM(共分散構造分析/構造方程式モデリング) | きにきじ: SEM を行うときのデータは『心理統計学の基礎:統合的理解のために』(南風原朝和)の第10章 第5節にある子どもの協調性の発… http://bit.ly/gt23Cz http://www.seo-log.net

      

    » このコメントを引用してコメントする

Leave a Reply


Copyright © 2008-2012 鍵山琢実 (KAGIYAMA, Takumi). All rights reserved.

This site's design was checked by IE 6.0+, Firefox 3.5+, GChrome 2.0+, Safari 4.0+, Opera 10.0+, and Sleipnir 2.8+ (all for Windows).
And JavaScript is used for some details. I am so sorry if your browser is not supported.

正当なCSSです! 私はチーム・マイナス6%です

↓ Today's My Favorite Phrase ↓

「てめーはこの空条承太郎がじきじきにブチのめす」

From: 荒木飛呂彦 『ジョジョの奇妙な冒険』第19巻 p. 180