Skip to first unread message
Brandon Bretl
Oct 9, 2019, 10:40:43 PM10/9/19
to lavaan
In checking 8 different factors for metric and scalar invariance, one of my factors (with 6 items) came back with a significant difference in the metric invariance test, suggesting it is non-invariant. Is this the end of the game for this factor? Or can I run modification indices, or what can I do to try and correct for metric invariance?
Mauricio Garnier-Villarreal
Oct 10, 2019, 12:58:40 AM10/10/19
to lavaan
The most common next step would be to test for partial invariance. This means to equate most of the items in the factor, but not all of them, allowing the misbehaving items to be different between groups. Modification indices and lavtestscore are good ways to identify which items cannot be equated.
In the Bayesian SEM framework you could use approximate invariance, where the misbehaving items would be allowed to be closer, but not exactly the same, by adjusting the variance in the difference between parameters. You can use blavaan for this that uses similar syntax as lavaan
Terrence Jorgensen
Oct 10, 2019, 6:38:00 AM10/10/19
to lavaan
can I run modification indices
That's only for fixed parameters, not constrained estimates. For the latter, use lavTestScore()
Or do what Mauricio suggested if it seems to be more than just one item with DIF.
Terrence D. Jorgensen
Assistant Professor, Methods and Statistics
Research Institute for Child Development and Education, the University of Amsterdam
Brandon Bretl
Oct 16, 2019, 3:42:50 PM10/16/19
to lavaan
Thanks for all your help. So this is what I did:
I ran a CFA of the factor that showed metric non-invariance:
ce <- '
ce =~ Q1_1 + Q2_1 + Q3_1 + Q4_1 + Q5_1 + Q6_1
'
ce1 <- cfa(ce, data = cfs,std.lv = TRUE, mimic = "mplus", estimator = "MLR")
summary(ce1, fit.measures = TRUE, rsquare=TRUE, standardized = TRUE)
modificationindices(ce1, standardized = TRUE,sort. = TRUE)Then I used your code to fit a metric invariance model and release parameters based my grouping variable "Q19":
## fit a metric-invariance model
fitmg <- update(ce1, group = "Q19", group.equal = "loadings")
## save the full parameter table
PT <- parTable(fitmg)
## print the part with loadings (to see labels) and constraints (to test)
PT[PT$op %in% c("=~","=="), ]
## pass the numbers (in the order they appear in the full PT) of the
## 6 constraints you want to test (i.e., 1:6) one at a time to lavTestScore()
lavTestScore(fitmg, release = 1, epc = TRUE)What am I looking for in the result? How can I tell if that item is causing the metric non-invariance and what will I do with that information?
$test
total score test:
test X2 df p.value
1 score 2.031 1 0.154
$uni
univariate score tests:
lhs op rhs X2 df p.value
1 .p1. == .p21. 2.031 1 0.154
$epc
expected parameter changes (epc) and expected parameter values (epv):
lhs op rhs block group free label plabel est epc epv sepc.lv sepc.all sepc.nox
1 ce =~ Q1_1 1 1 1 .p1. .p1. 0.308 -0.032 0.275 -0.032 -0.046 -0.046
2 ce =~ Q2_1 1 1 2 .p2. .p2. 0.661 0.004 0.665 0.004 0.004 0.004
3 ce =~ Q3_1 1 1 3 .p3. .p3. 0.682 0.005 0.686 0.005 0.005 0.005
4 ce =~ Q4_1 1 1 4 .p4. .p4. 0.677 0.004 0.681 0.004 0.004 0.004
5 ce =~ Q5_1 1 1 5 .p5. .p5. 0.501 0.003 0.504 0.003 0.004 0.004
6 ce =~ Q6_1 1 1 6 .p6. .p6. 0.692 0.005 0.696 0.005 0.005 0.005
7 Q1_1 ~~ Q1_1 1 1 7 .p7. 0.386 0.004 0.390 0.386 0.803 0.803
8 Q2_1 ~~ Q2_1 1 1 8 .p8. 0.450 -0.001 0.449 -0.450 -0.507 -0.507
9 Q3_1 ~~ Q3_1 1 1 9 .p9. 0.524 -0.002 0.522 -0.524 -0.530 -0.530
10 Q4_1 ~~ Q4_1 1 1 10 .p10. 0.769 -0.001 0.768 -0.769 -0.627 -0.627
11 Q5_1 ~~ Q5_1 1 1 11 .p11. 0.392 -0.001 0.392 -0.392 -0.610 -0.610
12 Q6_1 ~~ Q6_1 1 1 12 .p12. 0.432 -0.002 0.430 -0.432 -0.474 -0.474
13 ce ~~ ce 1 1 0 .p13. 1.000 NA NA NA NA NA
14 Q1_1 ~1 1 1 13 .p14. 4.657 0.000 4.657 0.000 0.000 0.000
15 Q2_1 ~1 1 1 14 .p15. 3.617 0.000 3.617 0.000 0.000 0.000
16 Q3_1 ~1 1 1 15 .p16. 3.528 0.000 3.528 0.000 0.000 0.000
17 Q4_1 ~1 1 1 16 .p17. 3.239 0.000 3.239 0.000 0.000 0.000
18 Q5_1 ~1 1 1 17 .p18. 4.470 0.000 4.470 0.000 0.000 0.000
19 Q6_1 ~1 1 1 18 .p19. 3.989 0.000 3.989 0.000 0.000 0.000
20 ce ~1 1 1 0 .p20. 0.000 NA NA NA NA NA
21 ce =~ Q1_1 2 2 19 .p1. .p21. 0.308 0.034 0.342 0.037 0.056 0.056
22 ce =~ Q2_1 2 2 20 .p2. .p22. 0.661 0.004 0.665 0.004 0.004 0.004
23 ce =~ Q3_1 2 2 21 .p3. .p23. 0.682 0.005 0.686 0.005 0.005 0.005
24 ce =~ Q4_1 2 2 22 .p4. .p24. 0.677 0.004 0.681 0.004 0.004 0.004
25 ce =~ Q5_1 2 2 23 .p5. .p25. 0.501 0.003 0.504 0.003 0.003 0.003
26 ce =~ Q6_1 2 2 24 .p6. .p26. 0.692 0.005 0.696 0.005 0.005 0.005
27 Q1_1 ~~ Q1_1 2 2 25 .p27. 0.340 -0.004 0.336 -0.340 -0.752 -0.752
28 Q2_1 ~~ Q2_1 2 2 26 .p28. 0.575 0.002 0.577 0.575 0.526 0.526
29 Q3_1 ~~ Q3_1 2 2 27 .p29. 0.466 0.002 0.468 0.466 0.458 0.458
30 Q4_1 ~~ Q4_1 2 2 28 .p30. 0.936 0.002 0.937 0.936 0.632 0.632
31 Q5_1 ~~ Q5_1 2 2 29 .p31. 0.545 0.001 0.545 0.545 0.646 0.646
32 Q6_1 ~~ Q6_1 2 2 30 .p32. 0.443 0.002 0.446 0.443 0.439 0.439
33 ce ~~ ce 2 2 31 .p33. 1.186 -0.033 1.154 -1.000 -1.000 -1.000
34 Q1_1 ~1 2 2 32 .p34. 4.633 0.000 4.633 0.000 0.000 0.000
35 Q2_1 ~1 2 2 33 .p35. 3.417 0.000 3.417 0.000 0.000 0.000
36 Q3_1 ~1 2 2 34 .p36. 3.506 0.000 3.506 0.000 0.000 0.000
37 Q4_1 ~1 2 2 35 .p37. 3.218 0.000 3.218 0.000 0.000 0.000
38 Q5_1 ~1 2 2 36 .p38. 4.207 0.000 4.207 0.000 0.000 0.000
39 Q6_1 ~1 2 2 37 .p39. 3.671 0.000 3.671 0.000 0.000 0.000
40 ce ~1 2 2 0 .p40. 0.000 NA NA NA NA NAMauricio Garnier-Villarreal
Oct 17, 2019, 12:48:22 AM10/17/19
to lavaan
univariate score tests:
lhs op rhs X2 df p.value
1 .p1. == .p21. 2.031 1 0.154
This shows the hypothesis test if the parameters p1 and p21 are equal, a large X2 and low p-value would indiate that that parameter most likely should not be equated
This case only shows 1 parameter because you set release=1, if you dont it will show you the inivariate test for all the equated parameters lavTestScore(fitmg, epc = TRUE)
This will help you find which indicator should not be held equal between groups
total score test:
test X2 df p.value
1 score 2.031 1 0.154
Shows the overall test, if all the equated parameters (shown in the univariate section) are equal
Brandon Bretl
Oct 17, 2019, 8:49:39 AM10/17/19
to lavaan
Thanks so much for your help with this. So this is what I got:
lhs op rhs X2 df p.value
1 .p1. == .p21. 2.031 1 0.154
2 .p2. == .p22. 7.876 1 0.005
3 .p3. == .p23. 0.030 1 0.863
4 .p4. == .p24. 0.901 1 0.342
5 .p5. == .p25. 3.644 1 0.056
6 .p6. == .p26. 1.425 1 0.233If I am reading this right, .p2. and .p22. should not be equated (and possibly .p5 and .p25.). Do you have any guidance for how I un-equate these in the model then?
I feel like a real bum here because my R and stats knowledge is limited, but this is helping me learn a lot, so I really appreciate it.
Here is the rest of the printout:
lhs op rhs block group free label plabel est epc epv sepc.lv sepc.all sepc.nox
1 ce =~ Q1_1 1 1 1 .p1. .p1. 0.308 -0.032 0.275 -0.032 -0.046 -0.046
2 ce =~ Q2_1 1 1 2 .p2. .p2. 0.661 0.061 0.722 0.061 0.065 0.065
3 ce =~ Q3_1 1 1 3 .p3. .p3. 0.682 0.004 0.686 0.004 0.004 0.004
4 ce =~ Q4_1 1 1 4 .p4. .p4. 0.677 0.027 0.704 0.027 0.025 0.025
5 ce =~ Q5_1 1 1 5 .p5. .p5. 0.501 -0.038 0.463 -0.038 -0.048 -0.048
6 ce =~ Q6_1 1 1 6 .p6. .p6. 0.692 -0.027 0.665 -0.027 -0.028 -0.028
7 Q1_1 ~~ Q1_1 1 1 7 .p7. 0.386 0.004 0.390 0.386 0.803 0.803
8 Q2_1 ~~ Q2_1 1 1 8 .p8. 0.450 -0.023 0.427 -0.450 -0.507 -0.507
9 Q3_1 ~~ Q3_1 1 1 9 .p9. 0.524 -0.001 0.523 -0.524 -0.530 -0.530
10 Q4_1 ~~ Q4_1 1 1 10 .p10. 0.769 -0.008 0.761 -0.769 -0.627 -0.627
11 Q5_1 ~~ Q5_1 1 1 11 .p11. 0.392 0.009 0.402 0.392 0.610 0.610
12 Q6_1 ~~ Q6_1 1 1 12 .p12. 0.432 0.012 0.443 0.432 0.474 0.474
13 ce ~~ ce 1 1 0 .p13. 1.000 NA NA NA NA NA
14 Q1_1 ~1 1 1 13 .p14. 4.657 0.000 4.657 0.000 0.000 0.000
15 Q2_1 ~1 1 1 14 .p15. 3.617 0.000 3.617 0.000 0.000 0.000
16 Q3_1 ~1 1 1 15 .p16. 3.528 0.000 3.528 0.000 0.000 0.000
17 Q4_1 ~1 1 1 16 .p17. 3.239 0.000 3.239 0.000 0.000 0.000
18 Q5_1 ~1 1 1 17 .p18. 4.470 0.000 4.470 0.000 0.000 0.000
19 Q6_1 ~1 1 1 18 .p19. 3.989 0.000 3.989 0.000 0.000 0.000
20 ce ~1 1 1 0 .p20. 0.000 NA NA NA NA NA
21 ce =~ Q1_1 2 2 19 .p1. .p21. 0.308 0.031 0.339 0.034 0.050 0.050
22 ce =~ Q2_1 2 2 20 .p2. .p22. 0.661 -0.078 0.584 -0.085 -0.081 -0.081
23 ce =~ Q3_1 2 2 21 .p3. .p23. 0.682 -0.002 0.679 -0.002 -0.002 -0.002
24 ce =~ Q4_1 2 2 22 .p4. .p24. 0.677 -0.033 0.644 -0.036 -0.029 -0.029
25 ce =~ Q5_1 2 2 23 .p5. .p25. 0.501 0.056 0.557 0.061 0.067 0.067
26 ce =~ Q6_1 2 2 24 .p6. .p26. 0.692 0.031 0.723 0.034 0.034 0.034
27 Q1_1 ~~ Q1_1 2 2 25 .p27. 0.340 -0.004 0.336 -0.340 -0.752 -0.752
28 Q2_1 ~~ Q2_1 2 2 26 .p28. 0.575 0.031 0.606 0.575 0.526 0.526
29 Q3_1 ~~ Q3_1 2 2 27 .p29. 0.466 0.002 0.468 0.466 0.458 0.458
30 Q4_1 ~~ Q4_1 2 2 28 .p30. 0.936 0.012 0.948 0.936 0.632 0.632
31 Q5_1 ~~ Q5_1 2 2 29 .p31. 0.545 -0.014 0.531 -0.545 -0.646 -0.646
32 Q6_1 ~~ Q6_1 2 2 30 .p32. 0.443 -0.015 0.429 -0.443 -0.439 -0.439
33 ce ~~ ce 2 2 31 .p33. 1.186 -0.007 1.180 -1.000 -1.000 -1.000Mauricio Garnier-Villarreal
Oct 17, 2019, 2:24:01 PM10/17/19
to lavaan
ce <- '
ce =~ c(l1,l1)*Q1_1 + c(l21,l22)*Q2_1 + c(l3,l3)*Q3_1 + c(l4,l4)*Q4_1 + c(l5,l5)*Q5_1 + c(l6,l6)*Q6_1
ce ~~ c(1,NA)*cee
'
ce1 <- cfa(ce, data = cfs,std.lv = TRUE, mimic = "mplus", std.lv=T,
estimator = "MLR", group="Q19")
Here I am setting the equating of parameters manually with labels. As you set the same label to a parameter for both groups lavaan will equate this parameter, c(l1,l1)*Q1_1 sets the label l1 for factor loading 1 for both groups. Based on the univariate results, I am un equating the factor loading c(l21,l22)*Q2_1 by switching it different labels for group 1 l21 and group 2 l22Brandon Bretl
Oct 23, 2019, 6:08:18 PM10/23/19
to lavaan
Got it. That worked. Now that I have metric invariance or partial metric invariance for all my factors, I want to compare factor means. Is there any problem with me just entering the code from freeing these parameters into the model syntax I use for comparing factor means? This is the syntax I use:
cfs_model_base <-
'
cpa =~ Q1_2 + Q1_8 + Q2_2 + Q3_2 + Q4_2 + Q5_2
cph =~ Q2_8 + Q3_8 + Q4_8 + Q5_9 + Q6_2 + Q6_9
ce =~ c(l1,l1)*Q1_1 + c(l21,l22)*Q2_1 + c(l3,l3)*Q3_1 + c(l4,l4)*Q4_1 + c(l5,l5)*Q5_1 + c(l6,l6)*Q6_1
ce ~~ c(1,NA)*ce
just =~ Q1_3 + Q2_3 + Q3_3 + Q4_3 + Q5_3 + Q6_3
lib =~ Q1_4 + Q2_4 + Q3_4 + Q4_4 + Q5_4 + Q6_4
auth =~ Q1_5 + Q2_5 + Q3_5 + Q4_5 + Q5_6 + Q6_5
loya =~ Q1_6 + Q2_6 + Q3_6 + Q4_6 + Q5_7 + Q6_7
sanc =~ Q1_7 + Q2_7 + Q3_7 + Q4_7 + Q5_8 + Q6_8
'
factor_compare <- cfa(cfs_model_base, data = cfs,
group.equal = c("loadings", "intercepts"),
group ="Q19", std.lv = TRUE, mimic = "mplus", estimator = "MLR")
summary(factor_compare, fit.measures = TRUE, rsquare=TRUE, standardized = TRUE)And I can compare the difference in intercepts then as the differences in means?:
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
cpa -0.519 0.080 -6.497 0.000 -0.405 -0.405
cph -0.431 0.073 -5.916 0.000 -0.362 -0.362
ce -0.464 0.074 -6.293 0.000 -0.374 -0.374
just -0.085 0.065 -1.322 0.186 -0.082 -0.082
lib -0.401 0.067 -5.967 0.000 -0.396 -0.396
auth -0.231 0.062 -3.718 0.000 -0.224 -0.224
loya -0.225 0.062 -3.598 0.000 -0.219 -0.219
sanc -0.203 0.071 -2.878 0.004 -0.188 -0.188Thus there is a significant difference between these two groups (gender) in all factors except "just" (justice)?
Continued thanks,
Brandon
Mauricio Garnier-Villarreal
Oct 26, 2019, 5:33:41 PM10/26/19
to lavaan
yes, after measurement invariance, you can compare latent parameters, like means, variances, and covariances
Terrence Jorgensen
Oct 28, 2019, 8:29:02 AM10/28/19
to lavaan
Now that I have metric invariance or partial metric invariance for all my factors, I want to compare factor means. Is there any problem with me just entering the code from freeing these parameters into the model syntax I use for comparing factor means?
It is highly dubious to constrain an item's intercepts to equality if its loadings are not constrained. That is analogous to estimating an interaction term (here, ce by Q19) without including a main effect (of Q19), if you had observed ce to include as a predictor in a regression model for item Q2_1.
Your syntax indicates you freed the loadings of Q2_1, so free their intercepts too. Either in the syntax (Q2_1 ~ c(int2.g1, int2.g2)*1 ) or using group.partial = "Q2_1~1"
Brandon Bretl
Oct 28, 2019, 11:13:55 AM10/28/19
to lavaan
Thanks so much. So I freed the intercepts for Q2_1 by adding that code to the syntax. I reran the analysis and this is what I got for intercepts:
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
cpa -0.263 0.089 -2.972 0.003 -0.254 -0.254
cph -0.187 0.083 -2.254 0.024 -0.188 -0.188
ce -0.262 0.089 -2.935 0.003 -0.231 -0.231
just 0.091 0.078 1.169 0.243 0.102 0.102
lib -0.306 0.086 -3.542 0.000 -0.311 -0.311
auth -0.053 0.077 -0.690 0.490 -0.055 -0.055
loya -0.135 0.081 -1.675 0.094 -0.141 -0.141
sanc 0.026 0.093 0.284 0.777 0.027 0.027The way I interpret this is that for factors cpa, cph, ce, and lib, there is a significant difference in factor means between the two groups.
This is essentially saying that boys have less of the factor responsible for ratings of moral intuitions violations related to physical harm to animals (cpa), physical harm to humans (cph), emotional harm to humans (ce), and liberty/autonomy (lib). There is no significant difference in factor means between girls and boys in justice, respect for authority, in-group loyalty, or sanctity/purity.
Brandon
Code I ran:
cfs_metric <-
'
cpa =~ Q1_2 + Q1_8 + Q2_2 + Q3_2 + Q4_2 + Q5_2
cph =~ Q2_8 + Q3_8 + Q4_8 + Q5_9 + Q6_2 + Q6_9
ce =~ c(l1,l1)*Q1_1 + c(l21,l22)*Q2_1 + c(l3,l3)*Q3_1 + c(l4,l4)*Q4_1 + c(l5,l5)*Q5_1 + c(l6,l6)*Q6_1
ce ~~ c(1,NA)*ce
just =~ Q1_3 + Q2_3 + Q3_3 + Q4_3 + Q5_3 + Q6_3
lib =~ Q1_4 + Q2_4 + Q3_4 + Q4_4 + Q5_4 + Q6_4
auth =~ Q1_5 + Q2_5 + Q3_5 + Q4_5 + Q5_6 + Q6_5
loya =~ Q1_6 + Q2_6 + Q3_6 + Q4_6 + Q5_7 + Q6_7
sanc =~ Q1_7 + Q2_7 + Q3_7 + Q4_7 + Q5_8 + Q6_8
Q2_8 ~~ Q3_8
Q1_1 ~~ Q5_1
Q4_3 ~~ Q6_3
Q4_4 ~~ Q6_4
Q1_5 ~~ Q2_5
Q5_8 ~~ Q6_8
auth =~ Q5_7
cph =~ Q6_7
Q4_2 ~~ Q5_9
Q2_1 ~ c(int2.g1, int2.g2)*1
'
factor_compare <- cfa(cfs_metric, data = cfs,
group.equal = c("loadings", "intercepts"),
group ="Q19", std.lv = TRUE, mimic = "mplus", estimator = "MLR")
summary(factor_compare, fit.measures = TRUE, rsquare=TRUE, standardized = TRUE)
Mauricio Garnier-Villarreal
Oct 31, 2019, 12:06:31 AM10/31/19
to lavaan
Assuming you are using an alpha = 0.05, yes that would be the interpretation. For further detail interpretation you can look at std.all, which will be the cohen d effect size , and you can get a more detail interpretation of the magnituded based on this
Brandon Bretl
Jan 22, 2020, 9:52:53 PM1/22/20
to lavaan
Hi Mauricio, I came back to reanalyze some more data, and this time I can't get the chi-square to change in my fit.metric by changing parameters like we did in the other model.
Here I am trying to test metric invariance for 4 groups at once.
Based on results form lavTestScore, I am trying to free the loading of the Q5_2 indicator in group 3. To do that I used this syntax:
cpa <- '
cpa =~ c(l1,l1,l1,l1)*Q1_2 + c(l2,l2,l2,l2)*Q1_8 + c(l31,l31,l31,l31)*Q2_2 +
c(l4,l4,l4,l4)*Q3_2 + c(l51,l51,l51,l51)*Q4_2 + c(l61,l61,l62,l61)*Q5_2
cpa ~~ c(1,1,NA,1)*cpa
'No matter what parameters I try and free, e.g., L61, L62), I keep getting the exact same chi-square change indicating metric non-invariance (25.905, p=.039).
Here is the code I'm using to rerun the analysis after trying to free the parameters as above:
syntax.config <- measEq.syntax(configural.model = mft_fact, data=mft_dat, parameterization="theta",
ID.fac="auto.fix.first", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4))
cat(as.character(syntax.config))
summary(syntax.config)
## threshold invariance
syntax.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
ID.fac = "auto.fix.first", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds")
## notice that constraining 4 thresholds allows intercepts and residual
## variances to be freely estimated in all but the first group & occasion
cat(as.character(syntax.thresh))
## print a summary of model features
summary(syntax.thresh)
## Fit a model to the data either in a subsequent step:
mod.config <- as.character(syntax.config)
fit.config <- cfa(mod.config, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## or in a single step:
fit.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
ID.fac = "auto.fix.first", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds", return.fit = TRUE)
## compare their fit to test threshold invariance
anova(fit.config, fit.thresh)
######
# METRIC INVARIANCE TESTING
######
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
ID.fac = "auto.fix.first", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings"))
summary(syntax.metric) # summarize model features
mod.metric <- as.character(syntax.metric) # save as text
cat(mod.metric) # print/view lavaan syntax
## fit model to data
fit.metric <- cfa(mod.metric, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of loadings, given equivalence of thresholds
anova(fit.thresh, fit.metric)
######
# SCALAR INVARIANCE TESTING
######
## scalar invariance
syntax.scalar <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
ID.fac = "auto.fix.first", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings",
"intercepts"))
summary(syntax.scalar) # summarize model features
mod.scalar <- as.character(syntax.scalar) # save as text
cat(mod.scalar) # print/view lavaan syntax
## fit model to data
fit.scalar <- cfa(mod.scalar, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of intercepts, given equal thresholds & loadings
anova(fit.metric, fit.scalar)
######
# INVARIANCE TESTING TABLE
######
## For a single table with all results, you can pass the models to
## summarize to the compareFit() function
compareFit(fit.config, fit.thresh, fit.metric, fit.scalar)
lavInspect(fit.metric, "cov.lv")Mauricio Garnier-Villarreal
Jan 22, 2020, 10:40:03 PM1/22/20
to lavaan
Brandon
A couple of differences from your previous code. In this case you change the identification method, from fixed variance to marker variable. Would recommend to switch to std.lv intstead of auto.fix.first
syntax.config <- measEq.syntax(configural.model = mft_fact, data=mft_dat, parameterization="theta",
ID.fac="std.lv", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4))
Second, when you add the factor loadings constraints except one (partial invariance), you still need to free the latent variance. This would be your weak invariance model, with equal factor loadings (except one), factor variance fix for group 1, and estimated latent variance for groups 2,3,4
cpa <- '
cpa =~ c(l1,l1,l1,l1)*Q1_2 + c(l2,l2,l2,l2)*Q1_8 + c(l31,l31,l31,l31)*Q2_2 +
c(l4,l4,l4,l4)*Q3_2 + c(l51,l51,l51,l51)*Q4_2 + c(l61,l61,l62,l61)*Q5_2
cpa ~~ c(1,NA,NA,NA)*cpa
'Brandon Bretl
Jan 24, 2020, 12:52:13 PM1/24/20
to lavaan
Thanks Mauricio, but I still get the same results:
My initial model...
cpa <- '
cpa =~ Q1_2 + Q1_8 + Q2_2 + Q3_2 + Q4_2 + Q5_2
'Gives these results...
Chi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit.config 36 15477 15827 46.739
fit.thresh 36 15477 15827 46.739 0.000 0
fit.metric 51 15473 15750 72.645 25.905 15 0.0390291 *
fit.scalar 66 15486 15690 115.471 42.826 15 0.0001675 ***cpa <- '
cpa =~ c(l1,l1,l1,l1)*Q1_2 + c(l2,l2,l2,l2)*Q1_8 + c(l31,l31,l31,l31)*Q2_2 +
c(l4,l4,l4,l4)*Q3_2 + c(l51,l51,l52,l51)*Q4_2 + c(l61,l61,l63,l61)*Q5_2
cpa ~~ c(1,NA,NA,NA)*cpa
'Gives the exact same results...
Chi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit.config 36 15477 15827 46.739
fit.thresh 36 15477 15827 46.739 0.000 0
fit.metric 51 15473 15750 72.645 25.905 15 0.0390291 *
fit.scalar 66 15486 15690 115.471 42.826 15 0.0001675 ***When I run it through the same metric and scalar invariance tests:
# add grouping variable here
g1 <- "dissgroup"
# add factor variable here
mft_fact <- cpa
# add data variable here
mft_dat <- cfs_middle
#Invariance testing
syntax.config <- measEq.syntax(configural.model = mft_fact, data=mft_dat, parameterization="theta",
ID.fac="std.lv", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4))
cat(as.character(syntax.config))
summary(syntax.config)
## threshold invariance
syntax.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds")
## notice that constraining 4 thresholds allows intercepts and residual
## variances to be freely estimated in all but the first group & occasion
cat(as.character(syntax.thresh))
## print a summary of model features
summary(syntax.thresh)
## Fit a model to the data either in a subsequent step:
mod.config <- as.character(syntax.config)
fit.config <- cfa(mod.config, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## or in a single step:
fit.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds", return.fit = TRUE)
## compare their fit to test threshold invariance
anova(fit.config, fit.thresh)
######
# METRIC INVARIANCE TESTING
######
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings"))
summary(syntax.metric) # summarize model features
mod.metric <- as.character(syntax.metric) # save as text
cat(mod.metric) # print/view lavaan syntax
## fit model to data
fit.metric <- cfa(mod.metric, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of loadings, given equivalence of thresholds
anova(fit.thresh, fit.metric)
######
# SCALAR INVARIANCE TESTING
######
## scalar invariance
syntax.scalar <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings",
"intercepts"))
summary(syntax.scalar) # summarize model features
mod.scalar <- as.character(syntax.scalar) # save as text
cat(mod.scalar) # print/view lavaan syntax
## fit model to data
fit.scalar <- cfa(mod.scalar, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of intercepts, given equal thresholds & loadings
anova(fit.metric, fit.scalar)
######
# INVARIANCE TESTING TABLE
######
## For a single table with all results, you can pass the models to
## summarize to the compareFit() function
compareFit(fit.config, fit.thresh, fit.metric, fit.scalar)Terrence Jorgensen
Jan 25, 2020, 4:08:03 AM1/25/20
to lavaan
I am trying to free the loading of the Q5_2 indicator in group 3. To do that I used this syntax:
Nowhere in your syntax do you use the model "cpa <- ..." that you specified. If you are relying on measEq.syntax() to write your syntax for you, it will always assume you are only providing a configural model (which is the name of the argument), so your labels will not be preserved. You can only use the group.partial= argument to free a parameter in all groups/occasions. If you want to free a parameter in only one group (or occasion), you can use the update() method to replace the specification of a particular parameter with your own custom syntax for that parameter (i.e., not the whole model, just 'cpa =~ c(l61,l61,l62,l61)*Q5_2'). See an example of this functionality at the bottom of this thread: https://github.com/simsem/semTools/issues/60
Brandon Bretl
Jan 31, 2020, 12:20:45 PM1/31/20
to lavaan
I tried using the update() function but have the same problem. I assign cpa (now cpa_model) to model_fact at the top (model_fact <- cpa_model).
I ran the following syntax:
g1 <- "dissgroup"
mft_fact <- cpa_model
mft_dat <- cfs_middle
#Invariance testing
syntax.config <- measEq.syntax(configural.model = mft_fact, data=mft_dat, parameterization="theta",
cat(as.character(syntax.config))
summary(syntax.config)
## threshold invariance
syntax.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds")
## notice that constraining 4 thresholds allows intercepts and residual
## variances to be freely estimated in all but the first group & occasion
cat(as.character(syntax.thresh))
## print a summary of model features
summary(syntax.thresh)
## Fit a model to the data either in a subsequent step:
mod.config <- as.character(syntax.config)
fit.config <- cfa(mod.config, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## or in a single step:
fit.thresh <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4), group.equal = "thresholds", return.fit = TRUE)
## compare their fit to test threshold invariance
anova(fit.config, fit.thresh)
######
# METRIC INVARIANCE TESTING
######
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings"))
summary(syntax.metric) # summarize model features
mod.metric <- as.character(syntax.metric) # save as text
cat(mod.metric) # print/view lavaan syntax
switched <- update(syntax.metric, change.syntax = "cpa ~~ c(1,NA,NA,NA)*cpa")
cat(as.character(switched))
cat(as.character(update(switched, change.syntax = "cpa =~ c(l61,l61,l62,l61)*Q5_2")))
## fit model to data
fit.metric <- cfa(mod.metric, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of loadings, given equivalence of thresholds
anova(fit.thresh, fit.metric)And I still get the same Chisq diff and P value:
Chi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit.thresh 36 15477 15827 46.739
fit.metric 51 15473 15750 72.645 25.905 15 0.03903 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Terrence Jorgensen
Feb 4, 2020, 9:08:48 AM2/4/20
to lavaan
cat(as.character(switched))
cat(as.character(update(switched, change.syntax = "cpa =~ c(l61,l61,l62,l61)*Q5_2")))
All you did was print something to the screen. If you want to save the change to the syntax object, you need to assign it to an object. You could also make all changes at once:
switched <- update(syntax.metric, change.syntax = c("cpa ~~ c(1,NA,NA,NA)*cpa",
"cpa =~ c(l61,l61,l62,l61)*Q5_2"))
fit.metric <- cfa(mod.metric, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
Why did you fit your original metric-model syntax to the data, instead of the changes you wanted?
mod.metric2 <- as.character(switched)
fit.metric <- cfa(mod.metric2, ...)Brandon Bretl
Feb 4, 2020, 8:47:07 PM2/4/20
to lavaan
I highly appreciate your patience and persistence in helping me. I am still getting the same output and it seems like my metric-model is not being changed by the update() function:
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mft_fact, data = mft_dat,
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = g1, group.label=c(1,2,3,4),
group.equal = c("thresholds","loadings"))
summary(syntax.metric) # summarize model features
mod.metric <- as.character(syntax.metric) # save as text
cat(mod.metric) # print/view lavaan syntax
switched <- update(syntax.metric, change.syntax =
c("cpa ~~ c(1,NA,NA,NA)*cpa","cpa =~ c(l61,l61,l62,l61)*Q5_2"))
cat(as.character(switched))
# save as text
mod.metric2 <- as.character(switched)
## fit model to data
fit.metric <- cfa(mod.metric2, data = mft_dat, group = g1, group.label=c(1,2,3,4),
parameterization = "theta")
## test equivalence of loadings, given equivalence of thresholds
anova(fit.thresh, fit.metric)
When I print the
cat(as.character(switched))I looks like the same lavaan syntax as before:
> cat(as.character(switched))
## LOADINGS:
cpa =~ c(NA, NA, NA, NA)*Q1_2 + c(lambda.1_1, lambda.1_1, lambda.1_1, lambda.1_1)*Q1_2
cpa =~ c(NA, NA, NA, NA)*Q1_8 + c(lambda.2_1, lambda.2_1, lambda.2_1, lambda.2_1)*Q1_8
cpa =~ c(NA, NA, NA, NA)*Q2_2 + c(lambda.3_1, lambda.3_1, lambda.3_1, lambda.3_1)*Q2_2
cpa =~ c(NA, NA, NA, NA)*Q3_2 + c(lambda.4_1, lambda.4_1, lambda.4_1, lambda.4_1)*Q3_2
cpa =~ c(NA, NA, NA, NA)*Q4_2 + c(lambda.5_1, lambda.5_1, lambda.5_1, lambda.5_1)*Q4_2
cpa =~ c(NA, NA, NA, NA)*Q5_2 + c(lambda.6_1, lambda.6_1, lambda.6_1, lambda.6_1)*Q5_2
## INTERCEPTS:
Q1_2 ~ c(nu.1.g1, nu.1.g2, nu.1.g3, nu.1.g4)*1 + c(NA, NA, NA, NA)*1
Q1_8 ~ c(nu.2.g1, nu.2.g2, nu.2.g3, nu.2.g4)*1 + c(NA, NA, NA, NA)*1
Q2_2 ~ c(nu.3.g1, nu.3.g2, nu.3.g3, nu.3.g4)*1 + c(NA, NA, NA, NA)*1
Q3_2 ~ c(nu.4.g1, nu.4.g2, nu.4.g3, nu.4.g4)*1 + c(NA, NA, NA, NA)*1
Q4_2 ~ c(nu.5.g1, nu.5.g2, nu.5.g3, nu.5.g4)*1 + c(NA, NA, NA, NA)*1
Q5_2 ~ c(nu.6.g1, nu.6.g2, nu.6.g3, nu.6.g4)*1 + c(NA, NA, NA, NA)*1
## UNIQUE-FACTOR VARIANCES:
Q1_2 ~~ c(NA, NA, NA, NA)*Q1_2 + c(theta.1_1.g1, theta.1_1.g2, theta.1_1.g3, theta.1_1.g4)*Q1_2
Q1_8 ~~ c(NA, NA, NA, NA)*Q1_8 + c(theta.2_2.g1, theta.2_2.g2, theta.2_2.g3, theta.2_2.g4)*Q1_8
Q2_2 ~~ c(NA, NA, NA, NA)*Q2_2 + c(theta.3_3.g1, theta.3_3.g2, theta.3_3.g3, theta.3_3.g4)*Q2_2
Q3_2 ~~ c(NA, NA, NA, NA)*Q3_2 + c(theta.4_4.g1, theta.4_4.g2, theta.4_4.g3, theta.4_4.g4)*Q3_2
Q4_2 ~~ c(NA, NA, NA, NA)*Q4_2 + c(theta.5_5.g1, theta.5_5.g2, theta.5_5.g3, theta.5_5.g4)*Q4_2
Q5_2 ~~ c(NA, NA, NA, NA)*Q5_2 + c(theta.6_6.g1, theta.6_6.g2, theta.6_6.g3, theta.6_6.g4)*Q5_2
## LATENT MEANS/INTERCEPTS:
cpa ~ c(alpha.1.g1, alpha.1.g2, alpha.1.g3, alpha.1.g4)*1 + c(0, 0, 0, 0)*1
## COMMON-FACTOR VARIANCES:
cpa ~~ c(1, NA, NA, NA)*cpa + c(psi.1_1.g1, psi.1_1.g2, psi.1_1.g3, psi.1_1.g4)*cpaTerrence Jorgensen
Feb 5, 2020, 3:12:36 PM2/5/20
to lavaan
I am still getting the same output and it seems like my metric-model is not being changed by the update() function:
Well, the factor variance matches, and the fact that loadings are constrained means the cpa~~cpa "change" would not be necessary anyway. But it is odd the factor loading labels aren't changing. Check sessionInfo(), did your version of semTools update to the one that has this functionality (0.5-2.914) or a later version? See my last comment in this thread for syntax to install the update (in a fresh version of R, to be safe):
Brandon Bretl
Feb 13, 2020, 5:31:50 PM2/13/20
to lavaan
Thanks so much, Terrence. That fixed it.
Brandon
Reply all
Reply to author
Forward
0 new messages