SEM - Hamaoka@fbc.keio.ac.jp's home page

構造方程式モデル

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では､一致することを確認したとある｡