I am trying to test the measurement invariance in my bi-factor model with CFA tests. I have 23 binary questions coded "correct" and "incorrect" and 3 different groups in my data. The sample sizes are approximately 200 for each group. To test the measurement invariance, I have followed the Van de Schoot et al. (2012) step-by-step process where progressively constrained models are compared. I have estimated the configural, metric and scalar models. The test statistics, fit indices and model results all seem good and reasonable. However, when I estimate the Scalar model, results give smaller degrees of freedom than the metric model. Is there any possible reason for this? I couldn't find anything on the internet although I looked for it carefully. Can someone please assist me in finding the reason or help me to find a reference to have a better understanding?
Here are the codes that I run.
### bi-factor configural model ###
model_mck_bifactor_con_2 <- 'MCK =~ Q1+Q2A+Q3B+Q3C+Q5+Q6+Q8A+Q8B+Q8C+Q10+Q12B+Q13A+Q13B+Q13C+Q17A+Q18B+Q19A+Q19B+Q19C+Q19D+Q20A+Q22+Q23
Q8 =~ Q8A+Q8B+Q8C; Q13 =~ Q13A+Q13B+Q13C; Q19 =~ Q19A+Q19B+Q19C+Q19D;
MCK~~0*Q8; MCK~~0*Q13; MCK~~0*Q19;
Q8~~0*Q13; Q8~~0*Q19; Q13~~0*Q19;
Q3B ~~ Q19A+Q19B+Q19C+Q19D; Q3C ~~ Q19A+Q19B+Q19C+Q19D; Q3B ~~ Q3C'
### Q3 and Q19 are theoretically correlated. ###
cfa.configural.mck.bifactor_2 <- cfa(model_mck_bifactor_con_2, data = data_mck, estimator = "WLSMV", ordered=T, group = "COUNTRY",
missing="pairwise", parameterization="delta")
### metric model ###
cfa.metric.mck.bifactor <- cfa(model_mck_bifactor_con_2, data = data_mck, estimator = "WLSMV", ordered=T, group = "COUNTRY", missing="pairwise",
parameterization="delta", group.equal="loadings")
summary(cfa.metric.mck.bifactor, fit.measures = TRUE, standardized = TRUE)
### scalar model ###
cfa.scalar.mck.bifactor <- cfa(model_mck_bifactor_con_2, data = data_mck, estimator = "WLSMV", ordered=T, group = "COUNTRY", missing="pairwise",
parameterization="delta", group.equal=c("loadings","thresholds"))
Below is the compareFit test results. As you can see, the scalar model has a smaller Df. This causes wrong calculations for chi-square difference tests between models.
Thanks a lot in advance.
Enes.
Dear Dr. Terrence Jorgensen,
Thanks for your response. I really appreciate your time and assistance. I have read both papers although Wu&Estabrook's article is a bit difficult to follow without sophisticated mathematical background as you said. I tried to implement the Sventina et al. (2020) work into my data, but it didn't work. I guess the problem was that my indicators were dichotomous. Sventina's indicators were ordinal with more than 3 categories. The measEq.syntax function didn't yield the required parameter constraints for my data. However, I tried to enter the parameters and constraints manually into CFA and run both the configural and scalar models as you advised. As far as I understand from Wu&Estabrook's article, testing the metric model is not statistically meaningful. In my example, constraining only the item loadings does not change much on the configural model. Additionally, the metric model yields more Degrees of Freedom than the scalar model since the metric model only constrains scale parameters (~*~) in the first group but lets them free for the other groups. This makes the comparison of metric and scalar models impossible. After all, if it's enough and ok to compare only the configural and scalar models in my case, there is no need to worry about this problem. I want to share with you again the codes and constraint parameters to make sure that I am doing everything right. Could you please take a look?
### configural model -------------------
bifactor_configural_model_wout_Q9 <- 'MCK =~ Q1+Q2A+Q3B+Q3C+Q5+Q6+Q8A+Q8B+Q8C+Q10+Q12B+Q13A+Q13B+Q13C+Q17A+Q18B+Q19A+Q19B+Q19C+Q19D+Q20A+Q22+Q23
Q8 =~ Q8B+Q8A+Q8C; Q13 =~ Q13B+Q13A+Q13C; Q19 =~ Q19C+Q19B+Q19A+Q19D;
MCK~~0*Q8; MCK~~0*Q13; MCK~~0*Q19;
Q8~~0*Q13; Q8~~0*Q19; Q13~~0*Q19;
Q3B ~~ Q19A+Q19B+Q19C+Q19D; Q3C ~~ Q19B+Q19C+Q19D; Q3B ~~ Q3C'
model.bifactor.configural.wout.Q9 <- cfa(bifactor_configural_model_wout_Q9, data = data_mck, estimator = "WLSMV", ordered=T, group = "COUNTRY",
missing="pairwise", parameterization="delta")
### scalar model ---------------------
bifactor_scalar_model_allconstrained <- 'MCK =~ c(V1,V1,V1)*Q1+c(V2A,V2A,V2A)*Q2A+c(V3B,V3B,V3B)*Q3B+c(V3C,V3C,V3C)*Q3C+c(V5,V5,V5)*Q5+c(V6,V6,V6)*Q6+c(V8A,V8A,V8A)*Q8A+c(V8B,V8B,V8B)*Q8B+c(V8C,V8C,V8C)*Q8C+c(V10,V10,V10)*Q10+c(V12B,V12B,V12B)*Q12B+c(V13A,V13A,V13A)*Q13A+c(V13B,V13B,V13B)*Q13B+c(V13C,V13C,V13C)*Q13C+c(V17A,V17A,V17A)*Q17A+c(V18B,V18B,V18B)*Q18B+c(V19A,V19A,V19A)*Q19A+c(V19B,V19B,V19B)*Q19B+c(V19C,V19C,V19C)*Q19C+c(V19D,V19D,V19D)*Q19D+c(V20A,V20A,V20A)*Q20A+c(V22,V22,V22)*Q22+c(V23,V23,V23)*Q23;
Q8 =~ c(VV8A,VV8A,VV8A)*Q8A+ c(VV8B,VV8B,VV8B)*Q8B+ c(VV8C,VV8C,VV8C)*Q8C; Q13 =~ c(VV13A,VV13A,VV13A)*Q13A+ c(VV13B,VV13B,VV13B)*Q13B+ c(VV13C,VV13C,VV13C)*Q13C; Q19 =~ c(VV19A,VV19A,VV19A)*Q19A+ c(VV19B,VV19B,VV19B)*Q19B+ c(VV19C,VV19C,VV19C)*Q19C+ c(VV19D,VV19D,VV19D)*Q19D;
### all factor loadings are fixed
MCK~~0*Q8; MCK~~0*Q13; MCK~~0*Q19;
Q8~~0*Q13; Q8~~0*Q19; Q13~~0*Q19;
Q3B ~~ Q19A+Q19B+Q19C+Q19D; Q3C ~~ Q19B+Q19C+Q19D; Q3B ~~ Q3C
### all thresholds are fixed
Q1 | c(tt1,tt1,tt1)*t1; Q2A | c(tt2,tt2,tt2)*t1; Q3B | c(tt3b,tt3b,tt3b)*t1; Q3C | c(tt3c,tt3c,tt3c)*t1; Q5 | c(tt5,tt5,tt5)*t1;
Q6 | c(tt6,tt6,tt6)*t1; Q8A | c(tt8a,tt8a,tt8a)*t1; Q8B | c(tt8b,tt8b,tt8b)*t1; Q8C | c(tt8c,tt8c,tt8c)*t1; Q10 | c(tt10,tt10,tt10)*t1;
Q12B | c(tt12b,tt12b,tt12b)*t1; Q13A | c(tt13a,tt13a,tt13a)*t1; Q13B | c(tt13b,tt13b,tt13b)*t1; Q13C | c(tt13c, tt13c,tt13c)*t1;
Q17A | c(tt17a,tt17a,tt17a)*t1; Q18B | c(tt18b,tt18b,tt18b)*t1; Q19A | c(tt19a,tt19a,tt19a)*t1; Q19B | c(tt19b,tt19b,tt19b)*t1;
Q19C | c(tt19c,tt19c,tt19c)*t1; Q19D | c(tt19d,tt19d,tt19d)*t1; Q20A | c(tt20a,tt20a,tt20a)*t1; Q22 | c(tt22,tt22,tt22)*t1; Q23 | c(tt23,tt23,tt23)*t1
### intercepts are fixed to 0 in the 1. group and freed for the 2. and 3. group
MCK~c(0,NA,NA)*1; Q8~c(0,NA,NA)*1; Q13~c(0,NA,NA)*1; Q19~c(0,NA,NA)*1
### scale parameters are fixed to 1.
Q1 ~*~ c(1,1,1)*Q1; Q2A ~*~ c(1,1,1)*Q2A; Q3B ~*~ c(1,1,1)*Q3B; Q3C ~*~ c(1,1,1)*Q3C; Q5 ~*~ c(1,1,1)*Q5; Q6 ~*~ c(1,1,1)*Q6; Q8A ~*~ c(1,1,1)*Q8A; Q8B ~*~ c(1,1,1)*Q8B; Q8C ~*~ c(1,1,1)*Q8C; Q10 ~*~ c(1,1,1)*Q10; Q12B ~*~ c(1,1,1)*Q12B; Q13A ~*~ c(1,1,1)*Q13A; Q13B ~*~ c(1,1,1)*Q13B; Q13C ~*~ c(1,1,1)*Q13C; Q17A ~*~ c(1,1,1)*Q17A; Q18B ~*~ c(1,1,1)*Q18B; Q19A ~*~ c(1,1,1)*Q19A; Q19B ~*~ c(1,1,1)*Q19B; Q19C ~*~ c(1,1,1)*Q19C; Q19D ~*~ c(1,1,1)*Q19D; Q20A ~*~ c(1,1,1)*Q20A; Q22 ~*~ c(1,1,1)*Q22; Q23 ~*~ c(1,1,1)*Q23'
### partial scalar model ----------------------------------
partial_scalar_model_allconstrained <- 'MCK =~ c(V1,V1,V1)*Q1+c(V2A,V2A,V2A)*Q2A+c(V3B,V3B,V3B)*Q3B+c(V3C,V3C,V3C)*Q3C+c(V5,V5,V5)*Q5+c(V6,V6,V6)*Q6+c(V8A,V8A,V8A)*Q8A+c(V8B,V8B,V8B)*Q8B+c(V8C,V8C,V8C)*Q8C+c(V10,V10,V10)*Q10+c(V12B,V12B,V12B)*Q12B+c(V13A,V13A,V13A)*Q13A+c(V13B,V13B,V13B)*Q13B+c(V13C,V13C,V13C)*Q13C+c(V17A,V17A,V17A)*Q17A+c(V18B,V18B,V18B)*Q18B+c(V19A,V19A,V19A)*Q19A+c(V19B,V19B,V19B)*Q19B+c(V19C,V19C,V19C)*Q19C+c(V19D,V19D,V19D)*Q19D+c(V20A,V20A,V20A)*Q20A+c(V22,V22,V22)*Q22+c(V23,V23,V23)*Q23;
Q8 =~ c(VV8A,VV8A,VV8A)*Q8A+ c(VV8B,VV8B,VV8B)*Q8B+ c(VV8C,VV8C,VV8C)*Q8C; Q13 =~ c(VV13A,VV13A,VV13A)*Q13A+ c(VV13B,VV13B,VV13B)*Q13B+ c(VV13C,VV13C,VV13C)*Q13C; Q19 =~ c(VV19A,VV19A,VV19A)*Q19A+ c(VV19B,VV19B,VV19B)*Q19B+ c(VV19C,VV19C,VV19C)*Q19C+ c(VV19D,VV19D,VV19D)*Q19D;
MCK~~0*Q8; MCK~~0*Q13; MCK~~0*Q19;
Q8~~0*Q13; Q8~~0*Q19; Q13~~0*Q19;
Q3B ~~ Q19A+Q19B+Q19C+Q19D; Q3C ~~ Q19B+Q19C+Q19D; Q3B ~~ Q3C
### thresholds of non-invariant items are freely estimated
Q1 | c(tt1,tt1,tt1)*t1; Q2A | c(NA,NA,NA)*t1; Q3B | c(tt3b,tt3b,tt3b)*t1; Q3C | c(NA,NA,NA)*t1; Q5 | c(NA,NA,NA)*t1;
Q6 | c(tt6,tt6,tt6)*t1; Q8A | c(NA,NA,NA)*t1; Q8B | c(tt8b,tt8b,tt8b)*t1; Q8C | c(tt8c,tt8c,tt8c)*t1; Q10 | c(NA,NA,NA)*t1;
Q12B | c(tt12b,tt12b,tt12b)*t1; Q13A | c(tt13a,tt13a,tt13a)*t1; Q13B | c(NA,NA,NA)*t1; Q13C | c(tt13c, tt13c,tt13c)*t1;
Q17A | c(tt17a,tt17a,tt17a)*t1; Q18B | c(tt18b,tt18b,tt18b)*t1; Q19A | c(tt19a,tt19a,tt19a)*t1; Q19B | c(NA,NA,NA)*t1;
Q19C | c(tt19c,tt19c,tt19c)*t1; Q19D | c(NA,NA,NA)*t1; Q20A | c(NA,NA,NA)*t1; Q22 | c(NA,NA,NA)*t1; Q23 | c(NA,NA,NA)*t1
### intercepts are fixed to 0 in the 1. group and freed for the 2. and 3. group
MCK~c(0,NA,NA)*1; Q8~c(0,NA,NA)*1; Q13~c(0,NA,NA)*1; Q19~c(0,NA,NA)*1
### scale parameters are fixed to 1.
Q1 ~*~ c(1,1,1)*Q1; Q2A ~*~ c(1,1,1)*Q2A; Q3B ~*~ c(1,1,1)*Q3B; Q3C ~*~ c(1,1,1)*Q3C; Q5 ~*~ c(1,1,1)*Q5; Q6 ~*~ c(1,1,1)*Q6; Q8A ~*~ c(1,1,1)*Q8A;
Q8B ~*~ c(1,1,1)*Q8B; Q8C ~*~ c(1,1,1)*Q8C; Q10 ~*~ c(1,1,1)*Q10; Q12B ~*~ c(1,1,1)*Q12B; Q13A ~*~ c(1,1,1)*Q13A; Q13B ~*~ c(1,1,1)*Q13B; Q13C ~*~ c(1,1,1)*Q13C; Q17A ~*~ c(1,1,1)*Q17A; Q18B ~*~ c(1,1,1)*Q18B; Q19A ~*~ c(1,1,1)*Q19A; Q19B ~*~ c(1,1,1)*Q19B; Q19C ~*~ c(1,1,1)*Q19C; Q19D ~*~ c(1,1,1)*Q19D; Q20A ~*~ c(1,1,1)*Q20A; Q22 ~*~ c(1,1,1)*Q22; Q23 ~*~ c(1,1,1)*Q23'
fit.comparison <- compareFit(model.bifactor.configural.wout.Q9, scalar.model.bifactor.allconstrained, partial.scalar.model.bifactor.allconstrained)
summary(fit.comparison)
#### Q2A, Q3C, Q5, Q8A, Q10, Q13B, Q19B, Q19D, Q20A, Q22, Q23 are non-invariant items. So, the thresholds of these items are freely estimated in the partial model.
#### As a result, 11 out of 23 items are non-invariant.
If you don't mind, I have a few questions about my work.
1-) Based on the results above, can I claim that the MCK test (except for the non-invariant 11 questions) is partially invariant? I have also conducted a Generalised Logistic Regression analysis of Magis et al (2011) to detect DIF items among multiple groups. The results of these analyses are quite identical to the CFA results. Non-invariant 11 items have also either moderate or large effect DIF in the GLR analysis results.
2-) Would it statistically make sense to compare CFA and GLR analyses or they are actually the same thing?
3-) Lastly, I want to estimate persons' ability scores and compare group differences in this MCK test. Can I estimate persons' ability scores by for example a 2PL MG-IRT model where invariant items are used as anchor items and non-invariant items are freed? Or do you think there is a better way to estimate persons' ability scores based on the partial invariant model?
I am sorry for the very long answer and for taking too much of your time. I have been struggling a lot with my data for a really long time. I would be very appreciative if could help.
Best Regards, Enes.
