Lesson6.CFA,SEM,Multiple Group SEM
Kosugitti
2015年11月27日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
この度ファイルをちょっと更新しました。2015がついていないbaseball,batter,pitcherデータを用意してください。
ピッチャーのデータ，pitcher.csvをb.dataというオブジェクトに入れているとします。 できてないひとは，次の通りやってみよう。
p.data <- read.csv("pitcher.csv")今日使うパッケージ
- psych
- lavaan
回帰分析と因子分析
前回「因子分析」の回で，目に見えないものが見えるようになりました。 例えば，身長と体重を足し合わせることで，「選手力」と呼ぶことができるのです。 そして，「選手力」が防御率ERAに影響しているかもしれない，ということを考えるのは自然なことですよね。
この因子分析と回帰分析を足し合わせたような方法として，パス解析の時に紹介したSEM(構造方程式モデリング，Structural Equation Modeling)を使うことが可能です。
SEMはlavaanパッケージを使うのでしたね。これを使って，選手力が年俸に影響しているかどうか，調べてみましょう。
library(lavaan)## Warning: package 'lavaan' was built under R version 3.2.2## This is lavaan 0.5-20
## lavaan is BETA software! Please report any bugs.model1 <- '
body =~ height + weight
ERA ~ body
'
result.fit <- sem(model1,data=p.data)
summary(result.fit,fit.measures=T,standardized=T)## lavaan (0.5-20) converged normally after 83 iterations
##
## Number of observations 66
##
## Estimator ML
## Minimum Function Test Statistic 0.000
## Degrees of freedom 0
##
## Model test baseline model:
##
## Minimum Function Test Statistic 52.614
## Degrees of freedom 3
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -545.536
## Loglikelihood unrestricted model (H1) -545.536
##
## Number of free parameters 6
## Akaike (AIC) 1103.073
## Bayesian (BIC) 1116.211
## Sample-size adjusted Bayesian (BIC) 1097.322
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.000
## P-value RMSEA <= 0.05 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## body =~
## height 1.000 7.107 0.949
## weight 1.162 0.443 2.620 0.009 8.255 0.752
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## ERA ~
## body -0.042 0.022 -1.929 0.054 -0.302 -0.303
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## height 5.631 18.488 0.305 0.761 5.631 0.100
## weight 52.467 26.533 1.977 0.048 52.467 0.435
## ERA 0.903 0.161 5.620 0.000 0.903 0.908
## body 50.505 20.865 2.421 0.015 1.000 1.000library(semPlot)## Warning: package 'semPlot' was built under R version 3.2.2semPaths(result.fit,what="stand",style="lisrel")
このように，潜在変数を作るときには「=~」で項を結びます。複数の潜在変数を置いてモデルを組んだり，潜在変数同士の関係を回帰的に表現したりすることだって可能です！
model2 <- '
strength =~ CG + games
weakness =~ DB + W.pitch + balk + ER
strength ~ weakness
'
result.fit2 <- sem(model2,data=p.data)
summary(result.fit2,fit.measures=T,standardized=T)## lavaan (0.5-20) converged normally after 43 iterations
##
## Number of observations 66
##
## Estimator ML
## Minimum Function Test Statistic 20.833
## Degrees of freedom 8
## P-value (Chi-square) 0.008
##
## Model test baseline model:
##
## Minimum Function Test Statistic 92.123
## Degrees of freedom 15
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.834
## Tucker-Lewis Index (TLI) 0.688
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1055.319
## Loglikelihood unrestricted model (H1) -1044.902
##
## Number of free parameters 13
## Akaike (AIC) 2136.637
## Bayesian (BIC) 2165.103
## Sample-size adjusted Bayesian (BIC) 2124.176
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.156
## 90 Percent Confidence Interval 0.075 0.240
## P-value RMSEA <= 0.05 0.021
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.088
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength =~
## CG 1.000 0.819 0.487
## games -14.223 4.509 -3.154 0.002 -11.654 -0.646
## weakness =~
## DB 1.000 1.613 0.570
## W.pitch 0.594 0.201 2.955 0.003 0.958 0.413
## balk 0.072 0.065 1.105 0.269 0.116 0.143
## ER 11.147 2.803 3.977 0.000 17.980 0.956
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength ~
## weakness 0.442 0.147 3.016 0.003 0.870 0.870
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 2.162 0.431 5.011 0.000 2.162 0.763
## games 189.355 53.906 3.513 0.000 189.355 0.582
## DB 5.406 1.058 5.108 0.000 5.406 0.675
## W.pitch 4.470 0.803 5.567 0.000 4.470 0.830
## balk 0.640 0.112 5.734 0.000 0.640 0.980
## ER 30.396 57.895 0.525 0.600 30.396 0.086
## strength 0.163 0.224 0.729 0.466 0.243 0.243
## weakness 2.601 1.137 2.289 0.022 1.000 1.000library(semPlot)
semPaths(result.fit2,what="stand",style="lisrel")
ここでlavaanの式の書き方をまとめておきましょう。
| 言いたいこと | 書き方 |
|---|---|
| 回帰(因果) | ~ |
| 因子(潜在変数) | =~ |
| 相関がある | ~~ |
| 切片 | ~1 |
| 直交 | ~0 |
修正指標 modification index
モデルが自由に組めるようになったわけですが，どうも適合度がうまく上がらない・・・ そんな時は，修正指数を参考にしてみましょう。
modificationindices(result.fit2)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 16 strength =~ DB 2.591 3.894 3.191 1.128 1.128
## 17 strength =~ W.pitch 5.436 -3.913 -3.207 -1.381 -1.381
## 18 strength =~ balk 0.410 -0.362 -0.297 -0.367 -0.367
## 19 strength =~ ER 0.501 23.159 18.977 1.009 1.009
## 23 CG ~~ DB 9.707 1.492 1.492 0.313 0.313
## 24 CG ~~ W.pitch 0.074 0.112 0.112 0.029 0.029
## 25 CG ~~ balk 0.962 -0.148 -0.148 -0.109 -0.109
## 26 CG ~~ ER 7.091 -12.294 -12.294 -0.388 -0.388
## 27 games ~~ DB 0.858 5.020 5.020 0.098 0.098
## 28 games ~~ W.pitch 7.524 11.679 11.679 0.279 0.279
## 29 games ~~ balk 0.002 0.067 0.067 0.005 0.005
## 30 games ~~ ER 8.277 -166.009 -166.009 -0.490 -0.490
## 31 DB ~~ W.pitch 0.026 0.109 0.109 0.017 0.017
## 32 DB ~~ balk 0.504 0.167 0.167 0.073 0.073
## 33 DB ~~ ER 4.668 -24.124 -24.124 -0.453 -0.453
## 34 W.pitch ~~ balk 2.126 0.307 0.307 0.164 0.164
## 35 W.pitch ~~ ER 2.641 9.736 9.736 0.223 0.223
## 36 balk ~~ ER 0.373 -1.026 -1.026 -0.068 -0.068これはどこにどのようなパスを仮定すれば，適合度が改善されるかを「機械的に」算出したものです。 機械的にというのを強調したのは，理論的にそこのパスはありえないよ，というものでも出してくるからです。 闇雲に機械のいうことを採用するのではなく，実際的に意味のある提案を採用するようにしましょう。
ところで，この一覧はわかりにくいので，修正指標の大きいじゅんに並べ替えておきましょうか。
mi <- modificationindices(result.fit2)
mi[order(mi$mi,decreasing = T),]## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 23 CG ~~ DB 9.707 1.492 1.492 0.313 0.313
## 30 games ~~ ER 8.277 -166.009 -166.009 -0.490 -0.490
## 28 games ~~ W.pitch 7.524 11.679 11.679 0.279 0.279
## 26 CG ~~ ER 7.091 -12.294 -12.294 -0.388 -0.388
## 17 strength =~ W.pitch 5.436 -3.913 -3.207 -1.381 -1.381
## 33 DB ~~ ER 4.668 -24.124 -24.124 -0.453 -0.453
## 35 W.pitch ~~ ER 2.641 9.736 9.736 0.223 0.223
## 16 strength =~ DB 2.591 3.894 3.191 1.128 1.128
## 34 W.pitch ~~ balk 2.126 0.307 0.307 0.164 0.164
## 25 CG ~~ balk 0.962 -0.148 -0.148 -0.109 -0.109
## 27 games ~~ DB 0.858 5.020 5.020 0.098 0.098
## 32 DB ~~ balk 0.504 0.167 0.167 0.073 0.073
## 19 strength =~ ER 0.501 23.159 18.977 1.009 1.009
## 18 strength =~ balk 0.410 -0.362 -0.297 -0.367 -0.367
## 36 balk ~~ ER 0.373 -1.026 -1.026 -0.068 -0.068
## 24 CG ~~ W.pitch 0.074 0.112 0.112 0.029 0.029
## 31 DB ~~ W.pitch 0.026 0.109 0.109 0.017 0.017
## 29 games ~~ balk 0.002 0.067 0.067 0.005 0.005これを参考に書き直したのが次のモデルです。
model3 <- '
strength =~ CG + games
weakness =~ DB + W.pitch + balk + ER
strength ~ weakness
games ~~ ER
'
result.fit3 <- sem(model3,data=p.data)
summary(result.fit3,fit.measures=T,standardized=T)## lavaan (0.5-20) converged normally after 71 iterations
##
## Number of observations 66
##
## Estimator ML
## Minimum Function Test Statistic 12.786
## Degrees of freedom 7
## P-value (Chi-square) 0.077
##
## Model test baseline model:
##
## Minimum Function Test Statistic 92.123
## Degrees of freedom 15
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.925
## Tucker-Lewis Index (TLI) 0.839
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1051.295
## Loglikelihood unrestricted model (H1) -1044.902
##
## Number of free parameters 14
## Akaike (AIC) 2130.590
## Bayesian (BIC) 2161.245
## Sample-size adjusted Bayesian (BIC) 2117.170
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.112
## 90 Percent Confidence Interval 0.000 0.208
## P-value RMSEA <= 0.05 0.138
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.070
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength =~
## CG 1.000 1.406 0.835
## games -5.307 3.094 -1.715 0.086 -7.461 -0.410
## weakness =~
## DB 1.000 1.942 0.686
## W.pitch 0.591 0.176 3.354 0.001 1.147 0.494
## balk 0.066 0.058 1.142 0.253 0.129 0.160
## ER 7.456 1.834 4.066 0.000 14.482 0.770
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength ~
## weakness 0.476 0.134 3.558 0.000 0.657 0.657
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## games ~~
## ER -119.140 40.429 -2.947 0.003 -119.140 -0.599
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 0.857 1.059 0.809 0.418 0.857 0.302
## games 275.166 56.512 4.869 0.000 275.166 0.832
## DB 4.235 1.073 3.948 0.000 4.235 0.529
## W.pitch 4.072 0.784 5.192 0.000 4.072 0.756
## balk 0.636 0.112 5.703 0.000 0.636 0.974
## ER 143.944 50.504 2.850 0.004 143.944 0.407
## strength 1.122 0.995 1.128 0.259 0.568 0.568
## weakness 3.772 1.417 2.663 0.008 1.000 1.000semPaths(result.fit3,what="stand",style="lisrel")
CFAとEFA
潜在変数を仮定するモデルで因子を想定することができましたが，前回の因子分析とどこが違うのでしょうか？
前回の因子分析では，因子数がいくつなのか，どの項目がどの因子に関係するのかがわからないまま手探りで分析を始めました。そのやり方を探索的因子分析;Exploratory Factor Analysis,EFA と言います。
SEMで因子を想定してモデルを書くときは，どの項目がどの因子に関係するのか，仮説的なモデルを当てはめるのですから，このやり方を検証的因子分析;Conformatory Factor Analysis:CFAと言います。
先行研究などで因子の構造がすでにわかっている場合は，CFAをすると良いでしょう。 先行研究がないとか，先行研究と違う項目が入っている，対象が随分変わっているなどの理由でゼロベースで考え直したい時はEFAをやりましょう。 あるいは，CFAでどうにもモデルの適合度が悪い場合は，EFAで本当に仮説モデルが再現できそうなのか，チェックしてみてもいいでしょう。
同じモデルを違う集団に
SEMを使ってモデルを書くことができる利点の一つは，同じモデルを複数のデータに当てはめることができるという点です。 これはモデルの妥当性を検証する場合に適しています。 同じモデルを当てはめる時に，様々な制約・統制をかけたりすることで，群間の違いがどこにあるかを検証することもできます。
この分析方法を多母集団同時分析と言います。
今回のモデルを，セリーグとパリーグで比較してみましょう。
result.fit4 <- sem(model3,data=p.data,group = "CP")
summary(result.fit4,fit.measures=T,standardized=T)## lavaan (0.5-20) converged normally after 147 iterations
##
## Number of observations per group
## P 29
## C 37
##
## Estimator ML
## Minimum Function Test Statistic 25.352
## Degrees of freedom 14
## P-value (Chi-square) 0.031
##
## Chi-square for each group:
##
## P 17.052
## C 8.299
##
## Model test baseline model:
##
## Minimum Function Test Statistic 119.410
## Degrees of freedom 30
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.873
## Tucker-Lewis Index (TLI) 0.728
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1022.290
## Loglikelihood unrestricted model (H1) -1009.614
##
## Number of free parameters 40
## Akaike (AIC) 2124.581
## Bayesian (BIC) 2212.167
## Sample-size adjusted Bayesian (BIC) 2086.239
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.157
## 90 Percent Confidence Interval 0.047 0.253
## P-value RMSEA <= 0.05 0.053
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.088
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
##
## Group 1 [P]:
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength =~
## CG 1.000 1.072 0.993
## games -9.016 5.659 -1.593 0.111 -9.669 -0.523
## weakness =~
## DB 1.000 2.059 0.682
## W.pitch 0.498 0.224 2.224 0.026 1.025 0.426
## balk 0.074 0.035 2.120 0.034 0.153 0.406
## ER 8.420 2.373 3.548 0.000 17.336 0.913
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength ~
## weakness 0.334 0.108 3.092 0.002 0.641 0.641
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## games ~~
## ER -110.228 53.763 -2.050 0.040 -110.228 -0.901
##
## Intercepts:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 0.724 0.200 3.612 0.000 0.724 0.671
## games 35.621 3.436 10.367 0.000 35.621 1.925
## DB 3.690 0.560 6.583 0.000 3.690 1.222
## W.pitch 2.931 0.447 6.561 0.000 2.931 1.218
## balk 0.172 0.070 2.458 0.014 0.172 0.456
## ER 33.172 3.526 9.407 0.000 33.172 1.747
## strength 0.000 0.000 0.000
## weakness 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 0.015 0.621 0.024 0.981 0.015 0.013
## games 248.893 82.621 3.012 0.003 248.893 0.727
## DB 4.871 1.513 3.220 0.001 4.871 0.535
## W.pitch 4.737 1.266 3.742 0.000 4.737 0.818
## balk 0.119 0.032 3.753 0.000 0.119 0.835
## ER 60.081 62.787 0.957 0.339 60.081 0.167
## strength 0.678 0.605 1.121 0.262 0.589 0.589
## weakness 4.239 2.177 1.947 0.052 1.000 1.000
##
##
## Group 2 [C]:
##
## Latent Variables:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength =~
## CG 1.000 1.720 0.857
## games -3.362 3.468 -0.970 0.332 -5.782 -0.324
## weakness =~
## DB 1.000 1.790 0.672
## W.pitch 0.762 0.271 2.809 0.005 1.363 0.607
## balk 0.070 0.108 0.646 0.518 0.125 0.123
## ER 7.204 2.478 2.907 0.004 12.891 0.706
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## strength ~
## weakness 0.595 0.239 2.492 0.013 0.619 0.619
##
## Covariances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## games ~~
## ER -113.018 51.294 -2.203 0.028 -113.018 -0.517
##
## Intercepts:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 1.243 0.330 3.771 0.000 1.243 0.620
## games 31.486 2.935 10.730 0.000 31.486 1.764
## DB 3.378 0.438 7.712 0.000 3.378 1.268
## W.pitch 2.649 0.369 7.177 0.000 2.649 1.180
## balk 0.351 0.168 2.097 0.036 0.351 0.345
## ER 38.892 3.004 12.948 0.000 38.892 2.129
## strength 0.000 0.000 0.000
## weakness 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
## CG 1.065 2.762 0.385 0.700 1.065 0.265
## games 285.195 73.500 3.880 0.000 285.195 0.895
## DB 3.898 1.291 3.019 0.003 3.898 0.549
## W.pitch 3.180 0.930 3.420 0.001 3.180 0.631
## balk 1.023 0.239 4.279 0.000 1.023 0.985
## ER 167.645 62.506 2.682 0.007 167.645 0.502
## strength 1.824 2.679 0.681 0.496 0.617 0.617
## weakness 3.202 1.658 1.931 0.053 1.000 1.000それぞれのグループでの影響力の違いなどを見ることができます。
本日の課題
- バッターのデータをつかって潜在変数をおいたパスモデルを分析してください。
- できたモデルを，セリーグとパリーグで比較する多母集団同時分析を行ってください。