lavaanでMplusのカテゴリカル因子分析を再現する
Mplusで行われたシミュレーションをRで再現しようとした際に調べたことの記録．MplusにおけるESTIMATOR=WLSMVとlavaanでestimator="WLSMV"にした際の挙動について．
データ
Rのpsychのパッケージに入っているbfiの質問紙データの一部（A1-A5の項目まで）を利用した．なぜ一部かというとMplusのデモ版だと使える変数の数に制限があるためである．とりあえずcsvにして書き出す．Mplusだとデータファイルは変数名を含めないようなのでcol_names=FALSEにする．欠損値処理についても考えないために今回はリストワイズ削除．N=2709の5項目（すべて6件法）のデータを使う．
library(readr)
library(lavaan)
bfiA <- psych::bfi[1:5]
bfiA <- na.omit(bfiA)
write_csv(bfiA, "bfiA.csv", col_names = FALSE)
head(bfiA)
# A1 A2 A3 A4 A5
#61617 2 4 3 4 4
#61618 2 4 5 2 5
#1620 5 4 5 4 4
#61621 4 4 6 5 5
#61622 2 3 3 4 5
#61623 6 6 5 6 5MplusでCFA
まずは，Mplusでやってみる．インプットファイルはこんな感じ．
TITLE: test;
DATA:
FILE = "bfiA.csv";
VARIABLE:
NAMES = A1-A5;
USEVARIABLES = A1-A5;
CATEGORICAL = A1-A5;
ANALYSIS:
TYPE = GENERAL;
ESTIMATOR= WLSMV;
MODEL:
F1 BY A1-A5*;
F1@1;
OUTPUT:
SAMPSTAT STDYX MODINDICES(ALL);
結果，モデルフィットの部分は次のような感じ．
MODEL FIT INFORMATION
Number of Free Parameters 30
Chi-Square Test of Model Fit
Value 143.590*
Degrees of Freedom 5
P-Value 0.0000
* The chi-square value for MLM, MLMV, MLR, ULSMV, WLSM and WLSMV cannot be used
for chi-square difference testing in the regular way. MLM, MLR and WLSM
chi-square difference testing is described on the Mplus website. MLMV, WLSMV,
and ULSMV difference testing is done using the DIFFTEST option.
RMSEA (Root Mean Square Error Of Approximation)
Estimate 0.101
90 Percent C.I. 0.087 0.116
Probability RMSEA <= .05 0.000
CFI/TLI
CFI 0.975
TLI 0.950
Chi-Square Test of Model Fit for the Baseline Model
Value 5560.971
Degrees of Freedom 10
P-Value 0.0000
SRMR (Standardized Root Mean Square Residual)
Value 0.022
Optimum Function Value for Weighted Least-Squares Estimator
Value 0.13676150D-01
標準化された推定値は次のような感じ．
STDYX Standardization
Two-Tailed
Estimate S.E. Est./S.E. P-Value
F1 BY
A1 -0.436 0.018 -24.498 0.000
A2 0.718 0.013 56.159 0.000
A3 0.810 0.012 69.528 0.000
A4 0.515 0.017 29.775 0.000
A5 0.668 0.013 50.400 0.000
Thresholds
A1$1 -0.441 0.025 -17.671 0.000
A1$2 0.321 0.025 13.100 0.000
A1$3 0.739 0.027 27.750 0.000
A1$4 1.232 0.032 38.440 0.000
A1$5 1.893 0.049 38.922 0.000
A2$1 -2.112 0.058 -36.114 0.000
A2$2 -1.526 0.038 -40.555 0.000
A2$3 -1.184 0.031 -37.785 0.000
A2$4 -0.475 0.025 -18.922 0.000
A2$5 0.485 0.025 19.301 0.000
A3$1 -1.840 0.047 -39.417 0.000
A3$2 -1.309 0.033 -39.306 0.000
A3$3 -0.952 0.028 -33.407 0.000
A3$4 -0.323 0.025 -13.176 0.000
A3$5 0.610 0.026 23.659 0.000
A4$1 -1.668 0.041 -40.447 0.000
A4$2 -1.143 0.031 -37.143 0.000
A4$3 -0.867 0.028 -31.318 0.000
A4$4 -0.366 0.025 -14.818 0.000
A4$5 0.236 0.024 9.694 0.000
A5$1 -2.018 0.054 -37.454 0.000
A5$2 -1.345 0.034 -39.633 0.000
A5$3 -0.911 0.028 -32.425 0.000
A5$4 -0.245 0.024 -10.077 0.000
A5$5 0.685 0.026 26.104 0.000
Variances
F1 1.000 0.000 999.000 999.000
lavaanでCFA
今度は同じデータをlavaan::cfaしてみる．Mplusの結果を再現するためにはestimator="WLSMV"にすることに加えて，mimic=“Mplus”`にすると良いらしい．
mo1 <- "F1 =~ A1+A2+A3+A4+A5"
res_lav <- cfa(mo1, data=bfiA, ordered=colnames(bfiA), std.lv=TRUE,
estimator="WLSMV", mimic="Mplus")
summary(res_lav, fit.measures=T)まずはモデルフィット．
lavaan 0.6-5 ended normally after 12 iterations
Estimator DWLS
Optimization method NLMINB
Number of free parameters 30
Number of observations 2709
Model Test User Model:
Standard Robust
Test Statistic 74.070 143.589
Degrees of freedom 5 5
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.517
Shift parameter 0.202
for the simple second-order correction (WLSMV)
Model Test Baseline Model:
Test statistic 7447.192 5560.976
Degrees of freedom 10 10
P-value 0.000 0.000
Scaling correction factor 1.340
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.991 0.975
Tucker-Lewis Index (TLI) 0.981 0.950
Robust Comparative Fit Index (CFI) NA
Robust Tucker-Lewis Index (TLI) NA
Root Mean Square Error of Approximation:
RMSEA 0.071 0.101
90 Percent confidence interval - lower 0.058 0.087
90 Percent confidence interval - upper 0.086 0.116
P-value RMSEA <= 0.05 0.006 0.000
Robust RMSEA NA
90 Percent confidence interval - lower NA
90 Percent confidence interval - upper NA
Standardized Root Mean Square Residual:
SRMR 0.031 0.031
Weighted Root Mean Square Residual:
WRMR 1.455 1.455
モデルフィットについてはほぼ一致していて（ずれていても.001程度），実用上問題なさそうである．ただSRMRだけMplusでは0.022でlavaanでは0.031とそこそこずれている．原因はよく分からない．
次にパラメータの推定値．
Parameter Estimates:
Information Expected
Information saturated (h1) model Unstructured
Standard errors Robust.sem
Latent Variables:
Estimate Std.Err z-value P(>|z|)
F1 =~
A1 0.436 0.018 24.488 0.000
A2 -0.718 0.013 -56.137 0.000
A3 -0.810 0.012 -69.503 0.000
A4 -0.515 0.017 -29.764 0.000
A5 -0.668 0.013 -50.381 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.A1 0.000
.A2 0.000
.A3 0.000
.A4 0.000
.A5 0.000
F1 0.000
Thresholds:
Estimate Std.Err z-value P(>|z|)
A1|t1 -0.441 0.025 -17.665 0.000
A1|t2 0.321 0.025 13.095 0.000
A1|t3 0.739 0.027 27.740 0.000
A1|t4 1.232 0.032 38.426 0.000
A1|t5 1.893 0.049 38.908 0.000
A2|t1 -2.112 0.058 -36.101 0.000
A2|t2 -1.526 0.038 -40.540 0.000
A2|t3 -1.184 0.031 -37.771 0.000
A2|t4 -0.475 0.025 -18.915 0.000
A2|t5 0.485 0.025 19.294 0.000
A3|t1 -1.840 0.047 -39.403 0.000
A3|t2 -1.309 0.033 -39.292 0.000
A3|t3 -0.952 0.029 -33.395 0.000
A3|t4 -0.323 0.025 -13.171 0.000
A3|t5 0.610 0.026 23.650 0.000
A4|t1 -1.668 0.041 -40.432 0.000
A4|t2 -1.143 0.031 -37.129 0.000
A4|t3 -0.867 0.028 -31.307 0.000
A4|t4 -0.366 0.025 -14.812 0.000
A4|t5 0.236 0.024 9.690 0.000
A5|t1 -2.018 0.054 -37.440 0.000
A5|t2 -1.345 0.034 -39.619 0.000
A5|t3 -0.911 0.028 -32.413 0.000
A5|t4 -0.245 0.024 -10.073 0.000
A5|t5 0.685 0.026 26.094 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.A1 0.810
.A2 0.485
.A3 0.344
.A4 0.735
.A5 0.554
F1 1.000
こちらもほぼほぼ一致している．
Mplusの昔の挙動
Mplusはversionが5までと6以降ではWLSMVの指す内容が違うらしい．具体的には検定統計量はversion 5までは，mean and variance adjusted test statistic (Satterthwaite type)で算出するのに対して6以降はcaled and shifted test statisticであるようだ．
https://personality-project.org/r/tutorials/summerschool.14/rosseel_sem_cat.pdf
Mplusで古いバージョンの挙動をさせるためにはANALYSISコマンドでSATTERTHWAITE=ONにすると良いようだ．lavaanでこれを再現するにはestimator="WLSMVS"にすると良い．
Mplusでこれを行うと次のようになる．（推定値は変わらないのでモデルフィットの部分のみ）
MODEL FIT INFORMATION
Number of Free Parameters 30
Chi-Square Test of Model Fit
Value 149.421*
Degrees of Freedom 5**
P-Value 0.0000
* The chi-square value for MLM, MLMV, MLR, ULSMV, WLSM and WLSMV cannot be used
for chi-square difference testing in the regular way. MLM, MLR and WLSM
chi-square difference testing is described on the Mplus website. MLMV, WLSMV,
and ULSMV difference testing is done using the DIFFTEST option.
** The degrees of freedom for MLMV, ULSMV and WLSMV are estimated according to
a formula given in the Mplus Technical Appendices at www.statmodel.com.
See degrees of freedom in the index of the Mplus User's Guide.
RMSEA (Root Mean Square Error Of Approximation)
Estimate 0.103
90 Percent C.I. 0.089 0.118
Probability RMSEA <= .05 0.000
CFI/TLI
CFI 0.968
TLI 0.961
Chi-Square Test of Model Fit for the Baseline Model
Value 4469.961
Degrees of Freedom 6
P-Value 0.0000
SRMR (Standardized Root Mean Square Residual)
Value 0.022
Optimum Function Value for Weighted Least-Squares Estimator
Value 0.13676150D-01
一方lavaanでestimator="WLSMVS"でやった結果がこちら．
avaan 0.6-5 ended normally after 12 iterations
Estimator DWLS
Optimization method NLMINB
Number of free parameters 30
Number of observations 2709
Model Test User Model:
Standard Robust
Test Statistic 74.070 149.420
Degrees of freedom 5 5
P-value (Chi-square) 0.000 0.000
Scaling correction factor 0.496
for the mean and variance adjusted correction (WLSMV)
Model Test Baseline Model:
Test statistic 7447.192 4469.965
Degrees of freedom 10 6
P-value 0.000 0.000
Scaling correction factor 1.666
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.991 0.968
Tucker-Lewis Index (TLI) 0.981 0.961
Robust Comparative Fit Index (CFI) NA
Robust Tucker-Lewis Index (TLI) NA
Root Mean Square Error of Approximation:
RMSEA 0.071 0.103
90 Percent confidence interval - lower 0.058 0.084
90 Percent confidence interval - upper 0.086 0.124
P-value RMSEA <= 0.05 0.006 0.000
Robust RMSEA NA
90 Percent confidence interval - lower NA
90 Percent confidence interval - upper NA
Standardized Root Mean Square Residual:
SRMR 0.031 0.031
Weighted Root Mean Square Residual:
WRMR 1.455 1.455
RMSEAの区間に若干のズレはあるものの概ね同等の値が得られている（SRMRは相変わらずずれているが）．