井出草平の研究ノート

トップ > R > lavaanを用いてWLSMVによる測定の不変性を計算する[R]

lavaanを用いてWLSMVによる測定の不変性を計算する[R]

R lavaan 因子分析 構造方程式モデリング 計量

WLSMV(adjusted diagonally weighted least squares)での測定の不変性の方法が確立しているらしい。

https://www.tandfonline.com/doi/abs/10.1080/10705511.2019.1602776

最尤法での測定の不変性についてはこちらを参照のこと。 ides.hatenablog.com

サンプルのデータセットやコマンドはこちらに掲載されている。 https://figshare.com/s/3f1d195da6c78195dd70

データは2011年TIMSS 4thのbullyingスケールから4項目を利用(Mullis et al., 2009)。恣意的に選ばれた3つの国（31＝アゼルバイジャン、40＝オーストリア、246＝フィンランド）のデータを使用。すべての項目は0（まったくない）から3（少なくとも週に一度はある）までの4段階のリッカート尺度で測定されている。この尺度の項目では1年の間に学校で以下のようなことがどのくらいの頻度で起こったかを尋ねている。サンプルサイズはアゼルバイジャン、オーストリア、フィンランドでそれぞれ3808、4457、4520だった。

library("lavaan")
library(semTools)
dat<-read.table("BULLY.dat", header=FALSE)
names(dat) <- c("IDCNTRY", "R09A", "R09B", "R09C", "R09D")
head(dat)

データの頭部分。

  IDCNTRY R09A R09B R09C R09D
1      31    3    3    0    0
2      31    0    0    0    0
3      31    3    2    1    3
4      31    0    0    3    0
5      31    0    0    0    0
6      31    0    0    0    0

グループ分けは"IDCNTRY"を用いる。

結果を格納するマトリックスを作成

all.results<-matrix(NA, ncol = 6, nrow=4)

ベースラインのモデルを作成

mod.cat <- 'F1 =~ R09A + R09B + R09C + R09D'
baseline <- measEq.syntax(configural.model = mod.cat,
                 data = dat, #結果に違いが出るので必ず記載
                 ordered = c("R09A", "R09B", "R09C", "R09D"),
                 parameterization = "delta",
                 ID.fac = "std.lv",
                 ID.cat = "Wu.Estabrook.2016",
                 group = "IDCNTRY",
                 group.equal = "configural")
summary(baseline)

cat(as.character(baseline)) # モデル内のすべての制約条件を見るために
model.baseline <- as.character(baseline) #lavaanで走らせるためにはas.characterを指定する必要がある。

fit.baseline <- cfa(model.baseline, data = dat, group = "IDCNTRY", ordered = c("R09A", "R09B", "R09C", "R09D")) 
summary(fit.baseline)

結果の格納

all.results[1,]<-round(
  data.matrix(
    fitmeasures(
      fit.baseline,fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled", "rmsea.scaled", "cfi.scaled", "tli.scaled")
               )
             ), digits=3
)

結果の格納マトリックスはこのような感じ。

       [,1] [,2] [,3]  [,4]  [,5]  [,6]
[1,] 50.944    6    0 0.042 0.997 0.991
[2,]     NA   NA   NA    NA    NA    NA
[3,]     NA   NA   NA    NA    NA    NA

閾値普遍性(Threshold invariance)

group.equal = c("thresholds")を設定する。

prop4 <- measEq.syntax(configural.model = mod.cat,
                        data = dat,
                        ordered = c("R09A", "R09B", "R09C", "R09D"),
                        parameterization = "delta",
                        ID.fac = "std.lv",
                        ID.cat = "Wu.Estabrook.2016",
                        group = "IDCNTRY",
                        group.equal = c("thresholds"))
model.prop4 <- as.character(prop4)
fit.prop4 <- cfa(model.prop4, data = dat, group = "IDCNTRY", ordered = c("R09A", "R09B", "R09C", "R09D")) 
summary(fit.prop4)

結果の格納

all.results[2,]<-round(
  data.matrix(
    fitmeasures(
      fit.prop4,fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled", "rmsea.scaled", "cfi.scaled", "tli.scaled")
               )
             ), digits=3
)

尤度比検定　ベースラインvs.閾値不変性モデル

lavTestLRT(fit.baseline,fit.prop4)

結果。

             Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
fit.baseline  6         26.942                                  
fit.prop4    14         42.170     46.132       8  2.243e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

閾値と負荷の不変性

group.equalに"loadings"を追加する。

prop7 <- measEq.syntax(configural.model = mod.cat,
                       data = dat,
                       ordered = c("R09A", "R09B", "R09C", "R09D"),
                       parameterization = "delta",
                       ID.fac = "std.lv",
                       ID.cat = "Wu.Estabrook.2016",
                       group = "IDCNTRY",
                       group.equal = c("thresholds", "loadings"))
model.prop7 <- as.character(prop7)
fit.prop7 <- cfa(model.prop7, data = dat, group = "IDCNTRY", ordered = c("R09A", "R09B", "R09C", "R09D")) 
summary(fit.prop7)

結果の格納

all.results[3,]<-round(
  data.matrix(
    fitmeasures(
      fit.prop7,fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled", "rmsea.scaled", "cfi.scaled", "tli.scaled")
               )
             ), digits=3
)

尤度比検定　閾値不変性 vs. 閾値と負荷の不変性

lavTestLRT(fit.prop4, fit.prop7)
lavTestLRT(fit.prop7, fit.baseline)

結果。

Scaled Chi-Squared Difference Test (method = “satorra.2000”)

lavaan NOTE:
    The “Chisq” column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
          Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
fit.prop4 14         42.170                                  
fit.prop7 20         93.115      54.67       6  5.404e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Scaled Chi-Squared Difference Test (method = “satorra.2000”)

fit.prop4とfit.prop7が有意であるため、制約を解放する。

修正指標を表示する。

mi <- modindices(fit.prop7, free.remove = FALSE)
mi[mi$op == "=~",]

結果。

   lhs op  rhs block group level    mi    epc sepc.lv sepc.all sepc.nox
1   F1 =~ R09A     1     1     1 0.493  0.012   0.012    0.012    0.012
2   F1 =~ R09B     1     1     1 0.080  0.005   0.005    0.005    0.005
3   F1 =~ R09C     1     1     1 2.698 -0.027  -0.027   -0.027   -0.027
4   F1 =~ R09D     1     1     1 0.662  0.014   0.014    0.014    0.014
31  F1 =~ R09A     2     2     1 0.162 -0.004  -0.002   -0.004   -0.004
32  F1 =~ R09B     2     2     1 0.366  0.007   0.004    0.007    0.007
33  F1 =~ R09C     2     2     1 0.008  0.001   0.001    0.001    0.001
34  F1 =~ R09D     2     2     1 0.177 -0.008  -0.004   -0.006   -0.006
61  F1 =~ R09A     3     3     1 0.001  0.000   0.000    0.000    0.000
62  F1 =~ R09B     3     3     1 0.695 -0.010  -0.004   -0.010   -0.010
63  F1 =~ R09C     3     3     1 1.164  0.013   0.005    0.013    0.013
64  F1 =~ R09D     3     3     1 0.219 -0.009  -0.004   -0.007   -0.007

miの値のもっとも大きいパスの制約を解放する。

prop7.part <- measEq.syntax(configural.model = mod.cat,
                       data = dat,
                       ordered = c("R09A", "R09B", "R09C", "R09D"),
                       parameterization = "delta",
                       ID.fac = "std.lv",
                       ID.cat = "Wu.Estabrook.2016",
                       group = "IDCNTRY",
                       group.equal = c("thresholds", "loadings"),
                       group.partial = "F1  =~ R09C")

model.prop7.part <- as.character(prop7.part)
fit.prop7.part <- cfa(model.prop7.part, data = dat, group = "IDCNTRY", ordered = c("R09A", "R09B", "R09C", "R09D")) 
summary(fit.prop7.part)

結果の格納

all.results[4,]<-round(
  data.matrix(
    fitmeasures(
      fit.prop7.part,fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled", "rmsea.scaled", "cfi.scaled", "tli.scaled")
               )
             ), digits=3
)

尤度比検定　閾値不変性 vs. 閾値と負荷の不変性パーシャルモデル

lavTestLRT(fit.prop7.part, fit.prop4)

結果。

Scaled Chi-Squared Difference Test (method = “satorra.2000”)

lavaan NOTE:
    The “Chisq” column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
               Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit.prop4      14         42.17                              
fit.prop7.part 18         47.73     6.8127       4     0.1461

結果の整理

column.names<-c("chisq.scaled","df.scaled","pvalue.scaled", "rmsea.scaled", "cfi.scaled", "tli.scaled" )
row.names<-c("baseline","prop4","prop7","prop7.partial" )

colnames(all.results)<-column.names
rownames(all.results)<-row.names

all.results
## 結果をcsvで出力する。
write.csv(all.results, file = "results.csv", row.names=TRUE, fileEncoding ="UTF-8")

結果。

              chisq.scaled df.scaled pvalue.scaled rmsea.scaled cfi.scaled tli.scaled
baseline            50.944         6             0        0.042      0.997      0.991
prop4              111.985        14             0        0.041      0.993      0.992
prop7              210.644        20             0        0.047      0.987      0.989
prop7.partial      110.487        18             0        0.035      0.994      0.994