Lesson6.CFA,SEM,Multiple Group SEM
Kosugitti
2015年11月27日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
今日のモデルのRQ例
- 風が吹けば桶屋が儲かるらしい。というのもつまり、大風で土ぼこりが立つ→土ぼこりが目に入って、盲人が増える→盲人は三味線を買う（当時の盲人が就ける職に由来）→三味線に使う猫皮が必要になり、ネコが殺される→ネコが減ればネズミが増える→ネズミは桶をかじる→桶の需要が増え桶屋が儲かる，らしい。本当かね？（パス解析）
- 女子力，コミュ力，社会人基礎力・・・いろいろな「力」があるらしいけど，これをどうやって測ればいいの？
- コミュ力があればさ，社会人基礎力だって必然的に身につくし，その結果として年収が高くなるってことだよね。
今日使うパッケージ
- psych
- lavaan
- semPlot
今日の話の準備
お配りしている野球のデータを見てみましょう。
baseball <- read.csv("baseball2016.csv",head=T)
summary(baseball)## name team pay position
## Ｔ−岡田 : 1 ソフトバンク:17 Min. : 4500 外野手:29
## エリアン : 1 日本ハム :15 1st Qu.: 7200 投手 :62
## エルドレッド: 1 巨人 :14 Median :10800 内野手:44
## エルナンデス: 1 西武 :13 Mean :14780 捕手 : 8
## ギャレット : 1 オリックス :12 3rd Qu.:20000
## クルーズ : 1 ヤクルト :12 Max. :60000
## (Other) :137 (Other) :60
## year age height weight
## Min. : 1.00 Min. :20.00 Min. :167 Min. : 66.00
## 1st Qu.: 4.00 1st Qu.:28.00 1st Qu.:177 1st Qu.: 80.00
## Median : 8.00 Median :31.00 Median :180 Median : 87.00
## Mean : 8.28 Mean :30.64 Mean :182 Mean : 87.76
## 3rd Qu.:12.00 3rd Qu.:33.50 3rd Qu.:185 3rd Qu.: 95.00
## Max. :21.00 Max. :41.00 Max. :203 Max. :122.00
##
## games Hit HR
## Min. : 13.00 Min. : 25.0 Min. : 0.00
## 1st Qu.: 27.50 1st Qu.: 60.0 1st Qu.: 3.00
## Median : 70.00 Median :104.0 Median : 8.00
## Mean : 77.93 Mean :102.3 Mean :10.19
## 3rd Qu.:121.00 3rd Qu.:140.0 3rd Qu.:14.00
## Max. :143.00 Max. :195.0 Max. :44.00
## table(baseball$team)##
## DeNA オリックス ソフトバンク ヤクルト ロッテ
## 9 12 17 12 11
## 楽天 巨人 広島 阪神 西武
## 9 14 11 10 13
## 中日 日本ハム
## 10 15table(baseball$position)##
## 外野手 投手 内野手 捕手
## 29 62 44 8ご覧の通り，チームが12球団，ポジションが色々あります。数では投手が一番多いですね。チームは大きく分けるとセントラル・リーグとパシフィック・リーグの二つになります。また，投手にとってのHR,Hitはホームランやヒットを「打たれた」ことですし，それ以外の選手にとっては打ったことになりますので，意味がちょっと違います。なので，これを区別する変数を新たに作ることにします。
# セ・リーグかパ・リーグか
c.league <- c("DeNA","ヤクルト","巨人","広島","阪神","中日")
baseball$CP <- ifelse(is.element(baseball$team,c.league),1,2)
baseball$CP <- factor(baseball$CP,labels=c("C","P"))
# 投手かどうか
baseball$Pos2 <- ifelse(baseball$position=="投手",1,2)
baseball$Pos2 <- factor(baseball$Pos2,labels=c("投手","野手"))
summary(baseball)## name team pay position
## Ｔ−岡田 : 1 ソフトバンク:17 Min. : 4500 外野手:29
## エリアン : 1 日本ハム :15 1st Qu.: 7200 投手 :62
## エルドレッド: 1 巨人 :14 Median :10800 内野手:44
## エルナンデス: 1 西武 :13 Mean :14780 捕手 : 8
## ギャレット : 1 オリックス :12 3rd Qu.:20000
## クルーズ : 1 ヤクルト :12 Max. :60000
## (Other) :137 (Other) :60
## year age height weight
## Min. : 1.00 Min. :20.00 Min. :167 Min. : 66.00
## 1st Qu.: 4.00 1st Qu.:28.00 1st Qu.:177 1st Qu.: 80.00
## Median : 8.00 Median :31.00 Median :180 Median : 87.00
## Mean : 8.28 Mean :30.64 Mean :182 Mean : 87.76
## 3rd Qu.:12.00 3rd Qu.:33.50 3rd Qu.:185 3rd Qu.: 95.00
## Max. :21.00 Max. :41.00 Max. :203 Max. :122.00
##
## games Hit HR CP Pos2
## Min. : 13.00 Min. : 25.0 Min. : 0.00 C:66 投手:62
## 1st Qu.: 27.50 1st Qu.: 60.0 1st Qu.: 3.00 P:77 野手:81
## Median : 70.00 Median :104.0 Median : 8.00
## Mean : 77.93 Mean :102.3 Mean :10.19
## 3rd Qu.:121.00 3rd Qu.:140.0 3rd Qu.:14.00
## Max. :143.00 Max. :195.0 Max. :44.00
## また，この後のエラー回避の為にちょっとだけ細工を。
baseball$pay <- baseball$pay/1000さて，今日は「野手」に限ってモデルを作ります。 勘弁のために，野手用のデータセットを作っておきましょう。
model.data <- subset(baseball,baseball$Pos2=="野手")回帰分析と因子分析
前回「因子分析」の回で，目に見えないものが見えるようになりました。 例えば，年齢と身長と体重を足し合わせることで，「選手力」と呼ぶことができるのです。 そして，「選手力」が年収に影響しているかもしれない，ということを考えるのは自然なことですよね。
この因子分析と回帰分析を足し合わせたような方法として，SEM(構造方程式モデリング，Structural Equation Modeling)を使うことが可能です。
SEMはlavaanパッケージを使います。これを使って，選手力が年俸に影響しているかどうか，調べてみましょう。
library(lavaan)## This is lavaan 0.5-22## lavaan is BETA software! Please report any bugs.model1 <- '
body =~ age + weight + height
pay ~ body
'
result.fit <- sem(model1,data=model.data)
summary(result.fit,fit.measures=T,standardized=T)## lavaan (0.5-22) converged normally after 43 iterations
##
## Number of observations 81
##
## Estimator ML
## Minimum Function Test Statistic 3.893
## Degrees of freedom 2
## P-value (Chi-square) 0.143
##
## Model test baseline model:
##
## Minimum Function Test Statistic 85.274
## Degrees of freedom 6
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.976
## Tucker-Lewis Index (TLI) 0.928
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1039.911
## Loglikelihood unrestricted model (H1) -1037.964
##
## Number of free parameters 8
## Akaike (AIC) 2095.821
## Bayesian (BIC) 2114.977
## Sample-size adjusted Bayesian (BIC) 2089.748
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.108
## 90 Percent Confidence Interval 0.000 0.269
## P-value RMSEA <= 0.05 0.199
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.061
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## age 1.000 0.822 0.200
## weight 12.724 7.485 1.700 0.089 10.461 0.973
## height 4.909 2.779 1.766 0.077 4.036 0.763
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## pay ~
## body 4.452 2.740 1.625 0.104 3.660 0.401
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .age 16.186 2.551 6.344 0.000 16.186 0.960
## .weight 6.154 22.107 0.278 0.781 6.154 0.053
## .height 11.715 3.769 3.108 0.002 11.715 0.418
## .pay 69.926 11.334 6.170 0.000 69.926 0.839
## body 0.676 0.770 0.878 0.380 1.000 1.000library(semPlot)
semPaths(result.fit,what="stand",style="lisrel")
このように，潜在変数を作るときには「=~」で項を結びます。複数の潜在変数を置いてモデルを組んだり，潜在変数同士の関係を回帰的に表現したりすることだって可能です！
model2 <- '
body =~ age + weight + height
performance =~ HR + pay
performance ~ body
'
result.fit2 <- sem(model2,data=model.data)
summary(result.fit2,fit.measures=T,standardized=T)## lavaan (0.5-22) converged normally after 162 iterations
##
## Number of observations 81
##
## Estimator ML
## Minimum Function Test Statistic 25.674
## Degrees of freedom 4
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 146.393
## Degrees of freedom 10
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.841
## Tucker-Lewis Index (TLI) 0.603
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1327.494
## Loglikelihood unrestricted model (H1) -1314.657
##
## Number of free parameters 11
## Akaike (AIC) 2676.988
## Bayesian (BIC) 2703.327
## Sample-size adjusted Bayesian (BIC) 2668.637
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.259
## 90 Percent Confidence Interval 0.169 0.358
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.095
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## age 1.000 0.764 0.186
## weight 13.755 8.397 1.638 0.101 10.505 0.977
## height 5.274 3.203 1.646 0.100 4.027 0.761
## performance =~
## HR 1.000 8.575 0.798
## pay 0.572 0.152 3.757 0.000 4.906 0.537
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## performance ~
## body 8.340 5.164 1.615 0.106 0.743 0.743
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .age 16.279 2.563 6.352 0.000 16.279 0.965
## .weight 5.231 13.545 0.386 0.699 5.231 0.045
## .height 11.783 2.718 4.335 0.000 11.783 0.421
## .HR 41.897 17.823 2.351 0.019 41.897 0.363
## .pay 59.254 10.774 5.500 0.000 59.254 0.711
## body 0.583 0.708 0.824 0.410 1.000 1.000
## .performance 32.958 17.892 1.842 0.065 0.448 0.448library(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 performance =~ age 5.889 -0.277 -2.373 -0.578 -0.578
## 17 performance =~ weight 0.016 -0.138 -1.185 -0.110 -0.110
## 18 performance =~ height 2.352 0.541 4.639 0.877 0.877
## 19 age ~~ weight 2.352 5.576 5.576 0.126 0.126
## 20 age ~~ height 0.016 -0.209 -0.209 -0.010 -0.010
## 21 age ~~ HR 18.158 -16.226 -16.226 -0.368 -0.368
## 22 age ~~ pay 9.042 10.928 10.928 0.292 0.292
## 23 weight ~~ height 5.889 -79.328 -79.328 -1.394 -1.394
## 24 weight ~~ HR 0.012 1.564 1.564 0.014 0.014
## 25 weight ~~ pay 0.029 -1.469 -1.469 -0.015 -0.015
## 26 height ~~ HR 1.093 5.684 5.684 0.100 0.100
## 27 height ~~ pay 0.125 -1.183 -1.183 -0.024 -0.024これはどこにどのようなパスを仮定すれば，適合度が改善されるかを「機械的に」算出したものです。 機械的にというのを強調したのは，理論的にそこのパスはありえないよ，というものでも出してくるからです。 闇雲に機械のいうことを採用するのではなく，実際的に意味のある提案を採用するようにしましょう。
ところで，この一覧はわかりにくいので，修正指標の大きい順に並べ替えておきましょうか。
mi <- modificationindices(result.fit2)
mi[order(mi$mi,decreasing = T),]## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 21 age ~~ HR 18.158 -16.226 -16.226 -0.368 -0.368
## 22 age ~~ pay 9.042 10.928 10.928 0.292 0.292
## 16 performance =~ age 5.889 -0.277 -2.373 -0.578 -0.578
## 23 weight ~~ height 5.889 -79.328 -79.328 -1.394 -1.394
## 19 age ~~ weight 2.352 5.576 5.576 0.126 0.126
## 18 performance =~ height 2.352 0.541 4.639 0.877 0.877
## 26 height ~~ HR 1.093 5.684 5.684 0.100 0.100
## 27 height ~~ pay 0.125 -1.183 -1.183 -0.024 -0.024
## 25 weight ~~ pay 0.029 -1.469 -1.469 -0.015 -0.015
## 20 age ~~ height 0.016 -0.209 -0.209 -0.010 -0.010
## 17 performance =~ weight 0.016 -0.138 -1.185 -0.110 -0.110
## 24 weight ~~ HR 0.012 1.564 1.564 0.014 0.014これを参考に書き直したのが次のモデルです。
model3 <- '
body =~ age + weight + height
performance =~ HR + pay
performance ~ body
age ~~ HR
age ~~ pay
'
result.fit3 <- sem(model3,data=model.data)
summary(result.fit3,fit.measures=T,standardized=T)## lavaan (0.5-22) converged normally after 147 iterations
##
## Number of observations 81
##
## Estimator ML
## Minimum Function Test Statistic 0.224
## Degrees of freedom 2
## P-value (Chi-square) 0.894
##
## Model test baseline model:
##
## Minimum Function Test Statistic 146.393
## Degrees of freedom 10
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.065
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1314.769
## Loglikelihood unrestricted model (H1) -1314.657
##
## Number of free parameters 13
## Akaike (AIC) 2655.538
## Bayesian (BIC) 2686.666
## Sample-size adjusted Bayesian (BIC) 2645.669
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent Confidence Interval 0.000 0.099
## P-value RMSEA <= 0.05 0.912
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.008
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## age 1.000 0.808 0.197
## weight 12.715 7.425 1.713 0.087 10.280 0.956
## height 5.084 2.965 1.714 0.086 4.110 0.777
## performance =~
## HR 1.000 8.639 0.804
## pay 0.564 0.151 3.739 0.000 4.869 0.533
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## performance ~
## body 8.078 5.269 1.533 0.125 0.756 0.756
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .age ~~
## .HR -14.115 4.219 -3.346 0.001 -14.115 -0.549
## .pay 7.214 3.869 1.864 0.062 7.214 0.232
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .age 16.208 2.556 6.342 0.000 16.208 0.961
## .weight 9.902 11.176 0.886 0.376 9.902 0.086
## .height 11.111 2.486 4.470 0.000 11.111 0.397
## .HR 40.791 17.988 2.268 0.023 40.791 0.353
## .pay 59.610 10.782 5.529 0.000 59.610 0.715
## body 0.654 0.761 0.859 0.390 1.000 1.000
## .performance 31.977 17.883 1.788 0.074 0.428 0.428semPaths(result.fit3,what="stand",style="lisrel")
年齢とホームラン，年齢と年収が「相関している」というのはちょっと意味がわからないですね。年齢が年収に別の影響を与えているということでしょうか。
model3 <- '
body =~ weight + height
performance =~ HR + pay
performance ~ body
pay ~ age
'
result.fit3 <- sem(model3,data=model.data)
summary(result.fit3,fit.measures=T,standardized=T)## lavaan (0.5-22) converged normally after 77 iterations
##
## Number of observations 81
##
## Estimator ML
## Minimum Function Test Statistic 16.524
## Degrees of freedom 4
## P-value (Chi-square) 0.002
##
## Model test baseline model:
##
## Minimum Function Test Statistic 146.393
## Degrees of freedom 10
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.908
## Tucker-Lewis Index (TLI) 0.770
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1322.919
## Loglikelihood unrestricted model (H1) -1314.657
##
## Number of free parameters 10
## Akaike (AIC) 2665.838
## Bayesian (BIC) 2689.783
## Sample-size adjusted Bayesian (BIC) 2658.246
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.197
## 90 Percent Confidence Interval 0.105 0.299
## P-value RMSEA <= 0.05 0.007
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.089
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## weight 1.000 10.163 0.945
## height 0.409 0.064 6.428 0.000 4.157 0.786
## performance =~
## HR 1.000 10.300 0.959
## pay 0.463 0.134 3.467 0.001 4.771 0.507
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## performance ~
## body 0.648 0.122 5.305 0.000 0.640 0.640
## pay ~
## age 0.801 0.204 3.927 0.000 0.801 0.349
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .weight 12.293 13.202 0.931 0.352 12.293 0.106
## .height 10.720 2.759 3.885 0.000 10.720 0.383
## .HR 9.328 24.898 0.375 0.708 9.328 0.081
## .pay 55.052 10.162 5.417 0.000 55.052 0.621
## body 103.286 22.286 4.635 0.000 1.000 1.000
## .performance 62.694 27.205 2.304 0.021 0.591 0.591semPaths(result.fit3,what="stand",style="lisrel")
ちょっとはましになったかな？
ちなみに，この最後のモデルは，＜選手力＞が＜パフォーマンス＞に影響し，年齢とともにそれが収入に影響を与える，という経路を説明しています。このようなステップを踏む回帰モデルのことをパス解析と言います。
エラーとモデルの改善
今回はこうして結果が得られましたが，実際のデータを使うと結構エラーがでて困ることがあります。少しエラーの解説をします。
答えが出ないヨゥ
lavaan WARNING: model has NOT converged! という警告は「モデルは収束しませんでした」という意味です。モデルが収束する，というのは「計算結果が出る」ということで，この警告が出たらそもそも答えが出ていない，出ていたとしても間違えた数字になっている，ということです。
理由はいくつか考えられます。
- モデルを識別するのに変数が少なすぎる。例えば一つの項目で一つの潜在変数を作ろうとしても，無理がありますよね。
- モデルを識別するのに条件が少なすぎる，多すぎる。これは識別の問題とよばれ，少し専門的な話になります。とりあえずモデルを書き換える必要があります。
- 計算が思っていた時間内に終わりそうにない。この問題であれば，計算時間の締め切りを伸ばすことで解決できることもあります。そもそもサンプルサイズが小さすぎるという問題があるかもしれません。
理論的におかしな数字になるヨゥ
lavaan WARNING: some estimated ov variances are negative これは「いくつかの分散（の推定値）が負の数になったよ」という意味です。分散は理論的に正の数字にしかならない(二乗和だから)のに，あなたのモデルで計算するとどうしてもいくつかが負の数にならないと理屈が合わない，という警告です。答えは一応示されることがありますが，信用してはいけません。
- モデルを書き換えましょう。
- サンプルサイズが小さすぎてうまくいかないのかもしれません。
単位がおかしくない？
lavaan WARNING: some observed variances are (at least) a factor 1000 times larger than others; use varTable(fit) to investigate これは訳すと「いくつかの観測変数が（少なくとも）他のものより1000倍の単位を持っているよ」という意味です。桁が違いすぎるので大丈夫かおい，ってことです。
今回は年収のデータがそれに当たっていました。ゲームに出る回数はせいぜい3桁なのに，年収は何千万，とかだったりするので，警告が出たのです。心理系の一般的なデータの場合はこういう問題があまり生じないから心配ないんですけど。
- 単位を整えましょう。
この他，SEMをやっていると警告やエラーが出ます。問題の根本的な解決のためには，SEMの数学的側面や理論的側面を学ばねばなりません。授業の範囲を超えますので，それらについてはこの本(http://amzn.to/2g61tlJ)などが参考になろうかと思います。
CFAとEFA
潜在変数を仮定するモデルで因子を想定することができましたが，前回の因子分析とどこが違うのでしょうか？
前回の因子分析では，因子数がいくつなのか，どの項目がどの因子に関係するのかがわからないまま手探りで分析を始めました。そのやり方を探索的因子分析;Exploratory Factor Analysis,EFA と言います。
SEMで因子を想定してモデルを書くときは，どの項目がどの因子に関係するのか，仮説的なモデルを当てはめるのですから，このやり方を検証的因子分析;Conformatory Factor Analysis:CFAと言います。
先行研究などで因子の構造がすでにわかっている場合は，CFAをすると良いでしょう。 先行研究がないとか，先行研究と違う項目が入っている，対象が随分変わっているなどの理由でゼロベースで考え直したい時はEFAをやりましょう。 あるいは，CFAでどうにもモデルの適合度が悪い場合は，EFAで本当に仮説モデルが再現できそうなのか，チェックしてみてもいいでしょう。
同じモデルを違う集団に
SEMを使ってモデルを書くことができる利点の一つは，同じモデルを複数のデータに当てはめることができるという点です。 これはモデルの妥当性を検証する場合に適しています。 同じモデルを当てはめる時に，様々な制約・統制をかけたりすることで，群間の違いがどこにあるかを検証することもできます。
この分析方法を多母集団同時分析と言います。
今回のモデルを，セリーグとパリーグで比較してみましょう。
result.fit4 <- sem(model3,data=baseball,group = "CP")## Warning in lav_object_post_check(lavobject): lavaan WARNING: some estimated
## ov variances are negativesummary(result.fit4,fit.measures=T,standardized=T)## lavaan (0.5-22) converged normally after 376 iterations
##
## Number of observations per group
## P 77
## C 66
##
## Estimator ML
## Minimum Function Test Statistic 16.661
## Degrees of freedom 8
## P-value (Chi-square) 0.034
##
## Chi-square for each group:
##
## P 3.172
## C 13.489
##
## Model test baseline model:
##
## Minimum Function Test Statistic 188.696
## Degrees of freedom 20
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.949
## Tucker-Lewis Index (TLI) 0.872
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2378.357
## Loglikelihood unrestricted model (H1) -2370.026
##
## Number of free parameters 28
## Akaike (AIC) 4812.714
## Bayesian (BIC) 4895.673
## Sample-size adjusted Bayesian (BIC) 4807.077
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.123
## 90 Percent Confidence Interval 0.033 0.207
## P-value RMSEA <= 0.05 0.077
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.062
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
##
## Group 1 [P]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## weight 1.000 9.857 0.966
## height 0.472 0.101 4.656 0.000 4.654 0.714
## performance =~
## HR 1.000 8.027 0.967
## pay 0.555 0.222 2.503 0.012 4.454 0.428
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## performance ~
## body 0.454 0.115 3.954 0.000 0.558 0.558
## pay ~
## age 1.092 0.241 4.530 0.000 1.092 0.418
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .weight 86.143 1.163 74.092 0.000 86.143 8.444
## .height 181.364 0.743 244.062 0.000 181.364 27.813
## .HR 8.974 0.946 9.483 0.000 8.974 1.081
## .pay -18.077 7.270 -2.486 0.013 -18.077 -1.735
## body 0.000 0.000 0.000
## .performance 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .weight 6.924 17.843 0.388 0.698 6.924 0.067
## .height 20.858 5.202 4.009 0.000 20.858 0.491
## .HR 4.531 21.977 0.206 0.837 4.531 0.066
## .pay 69.759 13.121 5.317 0.000 69.759 0.643
## body 97.160 24.439 3.976 0.000 1.000 1.000
## .performance 44.377 23.369 1.899 0.058 0.689 0.689
##
##
## Group 2 [C]:
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## body =~
## weight 1.000 6.127 0.558
## height 1.272 1.048 1.213 0.225 7.792 1.264
## performance =~
## HR 1.000 1.080 0.117
## pay 27.166 319.951 0.085 0.932 29.334 2.714
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## performance ~
## body 0.011 0.137 0.084 0.933 0.065 0.065
## pay ~
## age 0.964 0.253 3.805 0.000 0.964 0.389
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .weight 89.652 1.351 66.380 0.000 89.652 8.171
## .height 182.697 0.759 240.799 0.000 182.697 29.640
## .HR 11.606 1.134 10.235 0.000 11.606 1.260
## .pay -15.345 8.101 -1.894 0.058 -15.345 -1.420
## body 0.000 0.000 0.000
## .performance 0.000 0.000 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .weight 82.842 33.307 2.487 0.013 82.842 0.688
## .height -22.721 48.710 -0.466 0.641 -22.721 -0.598
## .HR 83.702 19.953 4.195 0.000 83.702 0.986
## .pay -761.344 10060.586 -0.076 0.940 -761.344 -6.515
## body 37.545 33.654 1.116 0.265 1.000 1.000
## .performance 1.161 13.730 0.085 0.933 0.996 0.996それぞれのグループでの影響力の違いなどを見ることができます。
本日の課題
「M-plusとRによる構造方程式モデリング入門」のサンプルデータを配布します。これは，友人や周囲の人から情緒的なサポートを受け取っていることと，友人や周囲の人に情緒的なサポートを提供していることが，関係満足感や主観的幸福感に影響を与えるという心理学的な仮説を検証するためのデータです。
データファイルには，ID,relation,V1-V14,statis,swbという変数が入っていますが，ここでIDは個人番号，relationは関係の違い(友人関係か恋愛関係か)，v1-v7はサポートを受け取ったかどうか，v8-v10はサポートを提供しているかどうかについての項目群，statisは関係満足スコア，swbは主観的幸福感スコアです。
これらを使って，次のことを確かめてください。
- v1,v2,v3,v4,v5,v6,v7をつかってサポート受け取り因子を，v8,v9,v10,v11,v12,v13,v14をつかってサポート提供因子を作り，この二つの因子が関係満足スコア，主観的幸福感を回帰しているモデルを作ってください。
- このモデルが関係の種類によって違うかどうかを検証する，多母集団同時分析モデルを作ってください。そのうえで，二群間の違いがどこにどのように見えるか，レポートしてください。
ちょっと高度ですが，頑張ってね。 それともう一つ。
- 何らかのデータを持ち寄って，統計的なモデリングの演習をしてみたいと思います。どういうデータを分析してみたいか，班単位？でディスカッションして，報告してください。
以上です。