RでSEM（共分散構造分析／構造方程式モデリング）
RでSEM（共分散構造分析／構造方程式モデリング）をやってみる。
SEMのツールといえばメジャーなのはなんといってもAmosでしょう。他にもCALIS(SAS)、EQS、 LISRELなどがよく使われているようです。
そしてRにもでもsemパッケージというのがあります。（実はOpenMxという別のパッケージも存在します。）
ただ、Amosにはまだ及ばないかな、というのが率直な印象です。
それでは、心理データ解析Ａに載っている三つの事例をsemパッケージで再現してみます。
なお、以下のページを参考にしました。
- Rで共分散構造分析・構造方程式モデル - RjpWiki
- RでSEM - 人の嫌がることを進んでする
- 共分散構造分析
- Structural Equation Modeling With the sem Package in R ※PDFなので注意
まずは単純なパス解析モデル。
# 分析例１（測定変数を用いたパス解析） # http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/10_folder/da10_02.html library(sem) 肯定感 <- c(34, 31, 30, 17, 13, 18, 17, 28, 23, 29, 33, 25, 37, 24, 33, 38, 23, 20, 26, 37) 積極性 <- c( 7, 6, 4, 3, 2, 4, 4, 4, 5, 3, 7, 3, 5, 6, 6, 8, 5, 3, 6, 3) 満足度 <- c( 8, 4, 3, 5, 2, 5, 2, 7, 5, 4, 8, 2, 6, 5, 7, 7, 3, 2, 8, 7) dat1.cor <- cor(cbind(肯定感, 積極性, 満足度)) model1 <- specify.model() 肯定感 -> 満足度, b01, NA 肯定感 -> 積極性, b02, NA 積極性 -> 満足度, b03, NA 肯定感 <-> 肯定感, e01, NA 積極性 <-> 積極性, e02, NA 満足度 <-> 満足度, e03, NA dat1.sem <- sem(model1, dat1.cor, N=20) # SEMの実行 std.coef(dat1.sem) # 標準化推定値 # Std. Estimate # b01 b01 0.4083722 満足度 <--- 肯定感 # b02 b02 0.5776581 積極性 <--- 肯定感 # b03 b03 0.4060430 満足度 <--- 積極性 # e01 e01 1.0000000 肯定感 <--> 肯定感 # e02 e02 0.6663111 積極性 <--> 積極性 # e03 e03 0.4767904 満足度 <--> 満足度 summary(dat1.sem) # GFI、RMSEA、BICなど # # Model Chisquare = -4.2188e-15 Df = 0 Pr(>Chisq) = NA # Chisquare (null model) = 21.787 Df = 3 # Goodness-of-fit index = 1 # BIC = -4.2188e-15 # # Normalized Residuals # Min. 1st Qu. Median Mean 3rd Qu. Max. # -1.37e-15 -1.22e-15 -1.22e-15 -6.57e-16 0.00e+00 3.42e-16 # # Parameter Estimates # Estimate Std Error z value Pr(>|z|) # b01 0.40837 0.19407 2.1043 0.0353522 満足度 <--- 肯定感 # b02 0.57766 0.18727 3.0847 0.0020378 積極性 <--- 肯定感 # b03 0.40604 0.19407 2.0923 0.0364118 満足度 <--- 積極性 # e01 1.00000 0.32451 3.0816 0.0020590 肯定感 <--> 肯定感 # e02 0.66631 0.21624 3.0813 0.0020611 積極性 <--> 積極性 # e03 0.47679 0.15476 3.0809 0.0020637 満足度 <--> 満足度 # # Iterations = 1
標準化推定値はAmosと一致しています。
R^2や非標準化推定値の出し方についてはまだ把握できていません。
パス図はpath.diagram関数とオープンソースのグラフ描画ツールGraphvizを使うことで描けるようなのですが、私の環境では上手くGraphvizが動いてくれなかったため割愛します。
続けて、本命の潜在変数を含んだモデル。
# 分析例２（潜在変数間の因果関係） # http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/10_folder/da10_03.html 勉強量a <- c( 5, 4, 4, 5, 4, 5, 5, 6, 4, 4, 3, 6, 5, 8, 5, 5, 6, 4, 3, 4, 6, 4, 3, 3, 6, 3, 5, 5, 3, 3) 勉強量b <- c( 6, 4, 7, 5, 5, 7, 6, 7, 5, 5, 3, 5, 6, 8, 6, 6, 5, 6, 4, 3, 5, 5, 3, 4, 7, 4, 6, 6, 4, 4) 期待a <- c( 2, 5, 6, 5, 4, 3, 3, 5, 6, 5, 5, 5, 4, 6, 7, 7, 6, 5, 5, 4, 4, 6, 5, 6, 4, 3, 7, 3, 6, 7) 期待b <- c( 3, 6, 5, 4, 5, 2, 5, 5, 6, 4, 4, 6, 5, 5, 7, 6, 6, 6, 4, 4, 5, 4, 4, 5, 5, 2, 5, 4, 7, 6) 自信a <- c(36, 51, 62, 50, 60, 50, 45, 62, 48, 44, 59, 55, 57, 58, 67, 58, 48, 47, 32, 25, 44, 45, 28, 36, 45, 35, 51, 54, 38, 60) 自信b <- c(31, 45, 41, 28, 38, 34, 31, 56, 45, 35, 42, 51, 40, 54, 60, 53, 45, 31, 23, 24, 38, 40, 33, 41, 39, 36, 43, 48, 26, 55) dat2.cor <- cor(cbind(勉強量a, 勉強量b, 期待a, 期待b, 自信a, 自信b)) model2 <- specify.model() 勉強量 -> 勉強量a, b01, 1 勉強量 -> 勉強量b, b02, NA 期待 -> 期待a, b03, 1 期待 -> 期待b, b04, NA 自信 -> 自信a, b05, 1 自信 -> 自信b, b06, NA 勉強量 -> 自信, b07, NA 勉強量 -> 期待, b08, NA 期待 -> 自信, b09, NA 勉強量 <-> 勉強量, NA, 1 勉強量b <-> 勉強量b, e01, NA 勉強量a <-> 勉強量a, e02, NA 期待b <-> 期待b, e03, NA 期待a <-> 期待a, e04, NA 自信a <-> 自信a, e05, NA 自信b <-> 自信b, e06, NA 自信 <-> 自信, NA, 1 期待 <-> 期待, NA, 1 dat2.sem <- sem(model2, dat2.cor, N=30) std.coef(dat2.sem) # Std. Estimate # 1 b01 0.80603730 勉強量a <--- 勉強量 # 2 b02 0.91774434 勉強量b <--- 勉強量 # 3 b03 0.87370849 期待a <--- 期待 # 4 b04 0.75329851 期待b <--- 期待 # 5 b05 0.87505495 自信a <--- 自信 # 6 b06 0.86500049 自信b <--- 自信 # 7 b07 0.57349021 自信 <--- 勉強量 # 8 b08 -0.02412674 期待 <--- 勉強量 # 9 b09 0.58473569 自信 <--- 期待 # 10 1.00000000 勉強量 <--> 勉強量 # 11 e01 0.15774533 勉強量b <--> 勉強量b # 12 e02 0.35030388 勉強量a <--> 勉強量a # 13 e03 0.43254136 期待b <--> 期待b # 14 e04 0.23663348 期待a <--> 期待a # 15 e05 0.23427884 自信a <--> 自信a # 16 e06 0.25177416 自信b <--> 自信b # 17 0.34537449 自信 <--> 自信 # 18 0.99941790 期待 <--> 期待 summary(dat2.sem) # # Model Chisquare = 8.8562 Df = 6 Pr(>Chisq) = 0.18183 # Chisquare (null model) = 91.764 Df = 15 # Goodness-of-fit index = 0.91965 # Adjusted goodness-of-fit index = 0.71878 # RMSEA index = 0.12812 90% CI: (NA, 0.29380) # Bentler-Bonnett NFI = 0.90349 # Tucker-Lewis NNFI = 0.90698 # Bentler CFI = 0.9628 # SRMR = 0.05172 # BIC = -11.551 # # Normalized Residuals # Min. 1st Qu. Median Mean 3rd Qu. Max. # -0.417000000 -0.051800000 0.000000461 0.050400000 0.239000000 0.775000000 # # Parameter Estimates # Estimate Std Error z value Pr(>|z|) # b01 0.806037 0.18904 4.26378 0.0000200993 勉強量a <--- 勉強量 # b02 0.917744 0.19000 4.83030 0.0000013633 勉強量b <--- 勉強量 # b03 0.873455 0.19455 4.48956 0.0000071372 期待a <--- 期待 # b04 0.753079 0.19351 3.89165 0.0000995645 期待b <--- 期待 # b05 0.514257 0.15138 3.39723 0.0006807200 自信a <--- 自信 # b06 0.508348 0.13990 3.63354 0.0002795575 自信b <--- 自信 # b07 0.975845 0.41237 2.36641 0.0179617202 自信 <--- 勉強量 # b08 -0.024134 0.22642 -0.10659 0.9151138679 期待 <--- 勉強量 # b09 0.994691 0.44183 2.25131 0.0243662171 自信 <--- 期待 # e01 0.157745 0.23689 0.66590 0.5054760827 勉強量b <--> 勉強量b # e02 0.350304 0.20208 1.73351 0.0830049283 勉強量a <--> 勉強量a # e03 0.432541 0.20336 2.12695 0.0334241808 期待b <--> 期待b # e04 0.236634 0.23525 1.00586 0.3144817342 期待a <--> 期待a # e05 0.234279 0.16160 1.44975 0.1471279356 自信a <--> 自信a # e06 0.251774 0.16029 1.57077 0.1162371372 自信b <--> 自信b # # Iterations = 48
標準化推定値はAmosと一致しています。
とりあえず、勉強量 -> 期待 のパスが有意になっていないので、これを外したモデルで再挑戦してみます。
model21 <- specify.model() 勉強量 -> 勉強量a, b01, 1 勉強量 -> 勉強量b, b02, NA 期待 -> 期待a, b03, 1 期待 -> 期待b, b04, NA 自信 -> 自信a, b05, 1 自信 -> 自信b, b06, NA 勉強量 -> 自信, b07, NA 期待 -> 自信, b09, NA 勉強量 <-> 勉強量, NA, 1 自信 <-> 自信, NA, 1 期待 <-> 期待, NA, 1 勉強量a <-> 勉強量a, e03, NA 勉強量b <-> 勉強量b, e04, NA 期待a <-> 期待a, e05, NA 期待b <-> 期待b, e06, NA 自信a <-> 自信a, e07, NA 自信b <-> 自信b, e08, NA dat21.sem <- sem(model21, dat2.cor, N=30) std.coef(dat21.sem) # Std. Estimate # 1 b01 0.8052671 勉強量a <--- 勉強量 # 2 b02 0.9186220 勉強量b <--- 勉強量 # 3 b03 0.8702293 期待a <--- 期待 # 4 b04 0.7563099 期待b <--- 期待 # 5 b05 0.8768200 自信a <--- 自信 # 6 b06 0.8655016 自信b <--- 自信 # 7 b07 0.5671534 自信 <--- 勉強量 # 8 b09 0.5804504 自信 <--- 期待 # 9 1.0000000 勉強量 <--> 勉強量 # 10 0.3414143 自信 <--> 自信 # 11 1.0000000 期待 <--> 期待 # 12 e03 0.3515449 勉強量a <--> 勉強量a # 13 e04 0.1561337 勉強量b <--> 勉強量b # 14 e05 0.2427009 期待a <--> 期待a # 15 e06 0.4279954 期待b <--> 期待b # 16 e07 0.2311867 自信a <--> 自信a # 17 e08 0.2509071 自信b <--> 自信b summary(dat21.sem) # # Model Chisquare = 8.8675 Df = 7 Pr(>Chisq) = 0.26231 # Chisquare (null model) = 91.764 Df = 15 # Goodness-of-fit index = 0.91927 # Adjusted goodness-of-fit index = 0.7578 # RMSEA index = 0.095915 90% CI: (NA, 0.26027) # Bentler-Bonnett NFI = 0.90337 # Tucker-Lewis NNFI = 0.94787 # Bentler CFI = 0.97567 # SRMR = 0.050651 # BIC = -14.941 # # Normalized Residuals # Min. 1st Qu. Median Mean 3rd Qu. # -0.4550000000 -0.0891000000 0.0000000109 0.0077500000 0.2010000000 # Max. # 0.6960000000 # # Parameter Estimates # Estimate Std Error z value Pr(>|z|) # b01 0.80527 0.18876 4.26604 0.0000198969 勉強量a <--- 勉強量 # b02 0.91862 0.18965 4.84379 0.0000012739 勉強量b <--- 勉強量 # b03 0.87023 0.19214 4.52922 0.0000059201 期待a <--- 期待 # b04 0.75631 0.19059 3.96825 0.0000724009 期待b <--- 期待 # b05 0.51444 0.15185 3.38776 0.0007046465 自信a <--- 自信 # b06 0.50775 0.13909 3.65064 0.0002615862 自信b <--- 自信 # b07 0.97064 0.40851 2.37605 0.0174991886 自信 <--- 勉強量 # b09 0.99340 0.44091 2.25305 0.0242562550 自信 <--- 期待 # e03 0.35154 0.20125 1.74679 0.0806740937 勉強量a <--> 勉強量a # e04 0.15613 0.23629 0.66076 0.5087645764 勉強量b <--> 勉強量b # e05 0.24270 0.22584 1.07464 0.2825360169 期待a <--> 期待a # e06 0.42800 0.19854 2.15567 0.0311091319 期待b <--> 期待b # e07 0.23310 0.16112 1.44676 0.1479638022 自信a <--> 自信a # e08 0.25293 0.15965 1.58422 0.1131446284 自信b <--> 自信b # # Iterations = 45
BIC、RMSEA、AGFIはこちらのモデルのほうが良好なようです。
最後に双方向因果を含むモデルを試してみます。
# 分析例３（双方向の因果関係） # http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/10_folder/da10_04.html 不信 <- c(22, 40, 43, 28, 29, 28, 33, 23, 23, 44, 38, 43, 24, 34, 32, 27, 37, 59, 26, 35, 40, 23, 49, 33, 36, 29, 32, 26, 31, 44, 10, 38, 33, 43, 21, 40, 34, 41, 26, 22, 39, 35, 49, 31, 47, 22, 31, 27, 28, 39, 26, 21, 33, 28, 34, 33, 36, 38, 42, 41, 37, 26, 37, 39, 31, 36, 38, 21, 29, 39, 25, 31, 41, 28, 29, 20, 25, 34, 39, 30, 32, 35, 36, 34, 26, 22, 35, 30, 27, 27, 24, 42, 30, 38, 44, 33, 32, 30, 30, 29) 自信頼 <- c(28, 27, 26, 24, 25, 28, 27, 26, 19, 25, 25, 29, 23, 30, 28, 27, 19, 25, 34, 22, 28, 24, 22, 27, 23, 26, 25, 34, 23, 21, 18, 25, 24, 14, 25, 19, 21, 30, 25, 23, 16, 12, 18, 25, 10, 22, 31, 30, 25, 24, 27, 25, 27, 26, 24, 25, 24, 27, 27, 28, 19, 27, 22, 26, 22, 32, 31, 22, 23, 26, 24, 25, 26, 27, 26, 24, 23, 27, 26, 20, 23, 21, 23, 20, 24, 25, 26, 22, 25, 24, 22, 18, 27, 22, 18, 20, 21, 28, 22, 27) 他信頼 <- c(37, 34, 27, 38, 32, 32, 33, 37, 36, 38, 32, 19, 37, 32, 43, 32, 23, 13, 39, 37, 34, 37, 29, 36, 29, 31, 34, 27, 28, 26, 36, 37, 32, 23, 36, 27, 28, 34, 31, 32, 32, 28, 36, 36, 26, 31, 35, 40, 32, 35, 40, 38, 33, 39, 36, 32, 31, 33, 35, 38, 35, 37, 28, 34, 29, 28, 42, 36, 36, 35, 35, 36, 28, 39, 34, 37, 34, 38, 33, 34, 38, 36, 36, 29, 37, 37, 38, 34, 35, 29, 39, 27, 39, 30, 33, 37, 28, 31, 28, 32) 肯定 <- c(22, 6, 23, 25, 17, 22, 31, 32, 32, 18, 28, 23, 11, 24, 33, 24, 20, 25, 26, 22, 20, 24, 20, 28, 20, 21, 21, 28, 14, 30, 28, 17, 25, 20, 19, 18, 23, 25, 22, 22, 19, 22, 27, 25, 25, 21, 21, 29, 23, 28, 24, 21, 27, 24, 24, 24, 25, 30, 21, 29, 26, 24, 25, 26, 25, 22, 23, 30, 22, 21, 21, 19, 22, 34, 20, 20, 21, 32, 23, 27, 22, 28, 23, 23, 27, 28, 24, 26, 22, 21, 28, 26, 28, 27, 36, 24, 13, 26, 25, 31) 否定 <- c(27, 35, 20, 26, 21, 32, 23, 15, 16, 30, 18, 22, 30, 21, 18, 25, 26, 23, 22, 23, 20, 22, 30, 20, 20, 29, 25, 17, 26, 7, 19, 25, 21, 24, 23, 25, 22, 24, 23, 24, 24, 17, 24, 24, 20, 17, 26, 17, 21, 27, 24, 26, 28, 22, 19, 20, 27, 24, 28, 19, 27, 14, 22, 22, 22, 22, 25, 18, 20, 27, 21, 25, 24, 17, 28, 25, 22, 19, 22, 18, 26, 22, 25, 26, 20, 20, 21, 19, 27, 27, 22, 21, 19, 17, 15, 26, 25, 21, 21, 20) 修復 <- c(15, 15, 16, 24, 18, 18, 26, 23, 28, 18, 25, 8, 19, 17, 26, 19, 18, 5, 25, 17, 18, 25, 16, 23, 14, 18, 20, 14, 21, 21, 30, 18, 19, 15, 19, 18, 18, 20, 16, 20, 17, 19, 14, 24, 20, 19, 20, 23, 22, 21, 21, 19, 18, 22, 22, 15, 20, 26, 17, 22, 18, 20, 22, 14, 20, 18, 22, 25, 27, 19, 22, 23, 17, 26, 21, 18, 21, 24, 21, 19, 24, 18, 18, 20, 22, 21, 20, 22, 25, 18, 22, 18, 23, 22, 30, 19, 20, 19, 18, 21) 崩壊 <- c(20, 24, 19, 15, 15, 20, 16, 6, 13, 22, 13, 15, 13, 13, 9, 15, 15, 20, 14, 14, 15, 14, 17, 16, 20, 13, 17, 19, 11, 9, 9, 15, 14, 21, 16, 15, 15, 14, 12, 16, 16, 12, 19, 9, 16, 18, 18, 10, 13, 12, 12, 20, 17, 11, 16, 20, 18, 15, 18, 16, 16, 8, 15, 15, 18, 17, 10, 10, 10, 18, 14, 18, 16, 11, 18, 20, 11, 5, 16, 16, 8, 13, 18, 13, 16, 13, 23, 13, 11, 15, 16, 16, 18, 13, 5, 19, 20, 13, 19, 22) dat3.cor <- cor(cbind(不信, 自信頼, 他信頼, 肯定, 否定, 修復, 崩壊)) model3 <- specify.model() 信頼感 -> 不信 , b01, NA 信頼感 -> 自信頼, b02, NA 信頼感 -> 他信頼, b03, 1 態度 -> 肯定, b04, 1 態度 -> 否定, b05, NA 志向性 -> 修復, b06, 1 志向性 -> 崩壊, b07, NA 信頼感 -> 志向性, b08, NA 志向性 -> 態度, b09, NA 態度 -> 志向性, b10, NA 不信 <-> 不信, e01, NA 自信頼 <-> 自信頼, e02, NA 他信頼 <-> 他信頼, e03, NA 肯定 <-> 肯定, e04, NA 否定 <-> 否定, e05, NA 修復 <-> 修復, e06, NA 崩壊 <-> 崩壊, e07, NA 態度 <-> 態度, e08, NA 志向性 <-> 志向性, e09, NA 信頼感 <-> 信頼感, NA, 1 dat3.sem <- sem(model3, dat3.cor, N=100) std.coef(dat3.sem) # Std. Estimate # b01 b01 -0.5930496 不信 <--- 信頼感 # b02 b02 0.2426302 自信頼 <--- 信頼感 # b03 b03 0.8365695 他信頼 <--- 信頼感 # b04 b04 0.8563447 肯定 <--- 態度 # b05 b05 -0.7298114 否定 <--- 態度 # b06 b06 0.9074742 修復 <--- 志向性 # b07 b07 -0.6123338 崩壊 <--- 志向性 # b08 b08 0.6451703 志向性 <--- 信頼感 # b09 b09 0.2353028 態度 <--- 志向性 # b10 b10 0.3882149 志向性 <--- 態度 # e01 e01 0.6482921 不信 <--> 不信 # e02 e02 0.9411306 自信頼 <--> 自信頼 # e03 e03 0.3001514 他信頼 <--> 他信頼 # e04 e04 0.2666738 肯定 <--> 肯定 # e05 e05 0.4673753 否定 <--> 否定 # e06 e06 0.1764905 修復 <--> 修復 # e07 e07 0.6250473 崩壊 <--> 崩壊 # e08 e08 0.7864972 態度 <--> 態度 # e09 e09 0.2908700 志向性 <--> 志向性 # 1.0000000 信頼感 <--> 信頼感 summary(dat3.sem) # # Model Chisquare = 33.771 Df = 9 Pr(>Chisq) = 0.000097978 # Chisquare (null model) = 217.11 Df = 21 # Goodness-of-fit index = 0.91241 # Adjusted goodness-of-fit index = 0.72751 # RMSEA index = 0.16674 90% CI: (0.10915, 0.22843) # Bentler-Bonnett NFI = 0.84446 # Tucker-Lewis NNFI = 0.70528 # Bentler CFI = 0.8737 # SRMR = 0.076056 # BIC = -7.6759 # # Normalized Residuals # Min. 1st Qu. Median Mean 3rd Qu. Max. # -1.8000000 -0.0000086 0.2170000 0.2440000 0.6730000 1.7300000 # # Parameter Estimates # Estimate Std Error z value Pr(>|z|) # b01 -0.59305 0.11604 -5.1109 0.000000320653195 不信 <--- 信頼感 # b02 0.24263 0.12168 1.9940 0.046152914923030 自信頼 <--- 信頼感 # b03 0.83657 0.12305 6.7986 0.000000000010563 他信頼 <--- 信頼感 # b04 0.86451 NaN NaN NaN 肯定 <--- 態度 # b05 -0.73677 NaN NaN NaN 否定 <--- 態度 # b06 0.77629 NaN NaN NaN 修復 <--- 志向性 # b07 -0.52381 NaN NaN NaN 崩壊 <--- 志向性 # b08 0.75420 NaN NaN NaN 志向性 <--- 信頼感 # b09 0.19939 NaN NaN NaN 態度 <--- 志向性 # b10 0.45814 NaN NaN NaN 志向性 <--- 態度 # e01 0.64829 0.12543 5.1687 0.000000235776981 不信 <--> 不信 # e02 0.94113 0.13813 6.8135 0.000000000009526 自信頼 <--> 自信頼 # e03 0.30015 0.16073 1.8674 0.061849798956312 他信頼 <--> 他信頼 # e04 0.26667 0.15161 1.7590 0.078579415632115 肯定 <--> 肯定 # e05 0.46738 0.12564 3.7200 0.000199218052814 否定 <--> 否定 # e06 0.17649 0.12833 1.3753 0.169042177848088 修復 <--> 修復 # e07 0.62505 0.10574 5.9110 0.000000003399387 崩壊 <--> 崩壊 # e08 0.77172 NaN NaN NaN 態度 <--> 態度 # e09 0.39748 NaN NaN NaN 志向性 <--> 志向性 # # Iterations = 25 # # Aliased parameters: b04 b05 b06 b07 b08 b09 b10 e08 e09 # 警告メッセージ： # In sqrt(diag(object$cov)) : 計算結果が NaN になりました
う〜ん。
標準化推定値は一致しているけれど、警告が出てるし、NaNもいっぱいある。
どこかモデルの書き方が間違っているのか、双方向因果にはパッケージが対応できていないのか……
今日はここまで。
いずれにせよ、対話的にモデルを構築するにはやっぱりAmosの方がべんりだな。
あとはOpenMXの今後に期待ということで。
