構造方程式モデル install.packages(c("lavaan","semTools","semPlot")) #としてインストールしておく 以下､実行例 library(lavaan) library(semTools) library(semPlot) #パス図表 を描く data(HolzingerSwineford1939) #サンプルデータ head(HolzingerSwineford1939) ?HolzingerSwineford1939 plot(HolzingerSwineford1939[,7:15]) cor(HolzingerSwineford1939[,7:15]) factanal(HolzingerSwineford1939[,7:15],factors=3,rotation="promax") #3因子で探索的因子分析 #測定方程式　　因子名は適当につけてok｡それと観測される変数を関連づける HS.model.cfa <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 ' #確認的因子分析 fit <- lavaan(HS.model.cfa, data=HolzingerSwineford1939,　 auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE) summary(fit, fit.measures=TRUE,standardized=T,rsquare=T) #結果出力 modindices(fit,sort=T) #修正指数　大きい順に出力 #参考　　まちがって次のように指定したらどうなるか　x3とx6を入れ替え HS.model.cfa2 <- ' visual =~ x1 + x2 + x6 textual =~ x4 + x5 + x3 speed =~ x7 + x8 + x9 ' #lavaan　関数に　上で定義したモデル名､データ名などを渡して推定 fit2 <- lavaan(HS.model.cfa2, data=HolzingerSwineford1939,auto.var=TRUE,auto.fix.first=TRUE,auto.cov.lv.x=TRUE) summary(fit2, fit.measures=TRUE,standardized=T,rsquare=T) modindices(fit2,sort=T) #修正指数 #上の確認的因子分析についてパスを入れる #測定方程式　　因子名は適当につけてok｡それと観測される変数を関連づける HS.model.sem <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 #構造方程式　　因子間の関係は下記のように指定 speed ~visual + textual ' fit.sem <- lavaan(HS.model.sem, data=HolzingerSwineford1939,　 auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE) summary(fit.sem, fit.measures=TRUE,standardized=T,rsquare=T) #結果出力 modindices(fit.sem,sort=T) #修正指数 #別のモデル　空間認知→言語認知→処理速度 HS.model.sem2 <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 #構造方程式　　因子間の関係は下記のように指定 speed ~ textual textual ~ visual ' fit.sem2 <- lavaan(HS.model.sem2, data=HolzingerSwineford1939,　 auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE) summary(fit.sem2, fit.measures=TRUE,standardized=T,rsquare=T) #結果出力 modindices(fit.sem2,sort=T) #修正指数 #参考}測定方程式のパスをすべて1に固定 HS.model.sem3 <- ' visual =~ 1*x1 + 1*x2 + 1*x3 textual =~ 1*x4 + 1*x5 + 1*x6 speed =~ 1*x7 + 1*x8 + 1*x9 #構造方程式　　因子間の関係は下記のように指定 speed ~ textual textual ~ visual ' fit.sem3 <- lavaan(HS.model.sem3, data=HolzingerSwineford1939,　 auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE) summary(fit.sem3, fit.measures=TRUE,standardized=T,rsquare=T) #結果出力 modindices(fit.sem3,sort=T) #修正指数 !!複数母集団の分析 #このデータは複数のschoolで測定されたもの｡複数母集団で推定してみる　school別に推定 table(HolzingerSwineford1939$school) #探索的因子分析をschool別に factanal(HolzingerSwineford1939[HolzingerSwineford1939$school=="Grant-White",7:15],factors=3,rotation="promax") #探索的因子分析　3因子指定 factanal(HolzingerSwineford1939[HolzingerSwineford1939$school=="Pasteur",7:15],factors=3,rotation="promax") #探索的因子分析　3因子指定 #2母集団でcfa　　指定しなければ同時推定するだけでパラメータは自由 fit.g <- lavaan(HS.model.cfa, data=HolzingerSwineford1939,auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE,group="school") summary(fit.g, fit.measures=TRUE) par(mfrow=c(1,2)) #2グループを並べてプロット ## semPaths(fit.g,whatLabels= "par") #パラメータを ?measurementInvariance #まずは測定の不変性を確認 measurementInvariance(model=HS.model.cfa, data=HolzingerSwineford1939,group="school",strict = TRUE) ***Measurement invariance models ***Model 1 : fit.configural　　パラメーターの制約無し ***Model 2 : fit.loadings　　　測定方程式のパス係数λ等値　(弱い測定の不変性) ***Model 3 : fit.intercepts　　さらに測定方程式(観測変数)の切片(平均値)等値　(強い測定の不変性) ***Model 4 : fit.residuals　　さらに残差を等値 ***Model 5 : fit.means　　　　さらに潜在変数の切片(平均値)等値 　 ***　　　 　Chi Square Difference Test ***　　　 　 Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq) ***fit.configural 48 7484.4 7706.8 115.85 ***fit.loadings 54 7480.6 7680.8 124.04 8.192 6 0.22436 制約してもモデルの適合度は有意には低下しない ***fit.intercepts 60 7508.6 7686.6 164.10 40.059 6 4.435e-07 ***　有意に低下→このように制約しないほうがよい ***fit.residuals 69 7508.1 7652.6 181.51 17.409 9 0.04269 * ***fit.means 72 7541.9 7675.3 221.34 39.824 3 1.161e-08 ***　　すべての等値制約が成立していれば､分けて推定する必要は無い ***--- 　***Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #測定方程式の負荷量λを等値制約 fit.g2 <- lavaan(model=HS.model.cfa, data=HolzingerSwineford1939,auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE,group="school",group.equal=c("loadings")) summary(fit.g2, fit.measures=TRUE,standardized=T,rsquare=T) modindices(fit.g2,sort=T) # par(mfrow=c(1,2)) #2グループを並べてプロット ## semPaths(fit.g2,whatLabels= "par") #パラメータを #λを等値制約したパス図表 で概念間のパスも入れる fit.g.s <- lavaan(model=HS.model.sem, data=HolzingerSwineford1939,auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE,group="school",group.equal=c("loadings")) summary(fit.g.s, fit.measures=TRUE,standardized=T,rsquare=T) #regressionの係数が同じになりそうか確認 modindices(fit.g.s,sort=T) # par(mfrow=c(1,2)) #2グループを並べてプロット ## semPaths(fit.g.s,whatLabels= "par") #エラー､､､､ #λを等値制約したパス図表 で概念間のパスも入れる(等値制約 fit.g.s2 <- lavaan(model=HS.model.sem, data=HolzingerSwineford1939,auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE,group="school",group.equal=c("loadings","regressions")) summary(fit.g.s2, fit.measures=TRUE,standardized=T,rsquare=T) modindices(fit.g.s2,sort=T) # par(mfrow=c(1,2)) #2グループを並べてプロット ## semPaths(fit.g.s2,whatLabels= "par") #エラー､､､､ #二つのグループで一部の係数が異なるとしたいとき HS.model.sem2 <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 #構造方程式　　因子間の関係は下記のように指定 #group.equal=c("loadings","regressions"))としたので回帰部分の係数はグループ間で等値とされる｡ #ただし､下のように指定すると､visualへの係数はグループで異なると推定される｡ speed ~c(b11,b12)*visual + textual ' fit.g.2 <- lavaan(model=HS.model.sem2, data=HolzingerSwineford1939,auto.var=TRUE, auto.fix.first=TRUE,auto.cov.lv.x=TRUE,group="school",group.equal=c("loadings","regressions")) summary(fit.g.2, fit.measures=TRUE,standardized=T,rsquare=T) #regressionの出力部分を比較して､そのようになっていることを確認｡AIC,BICなどで､適合度の高いモデルを選択 modindices(fit.g.2,sort=T) # par(mfrow=c(1,2)) #2グループを並べてプロット ## semPaths(fit.g.2,whatLabels= "par") #エラー､､､､関数廃止 **推定結果がEQS　Mplusなどと一致するかは確認していないが､同プロジェクトの[HP|http://lavaan.ugent.be/]では､一致することを確認したとある｡