latent change score model with 2 intervention groups: error with variance and is there change?
307 views
Skip to first unread message
cezre...@gmail.com
Jul 24, 2017, 11:22:47 PM7/24/17
to lavaan
Hello,
I am trying to fit a change score model. I have 2 cognitive indicator variables measured at pre and post for 2 intervention groups. The 2 indicators form a latent cognitive construct (working memory). I have also created a latent variable to model the change (similar to McArdle, 2009). Ultimately, I want to determine if the change from pre to post is significant. I copy my code and output below (and I also attached the path diagram result for both groups, if you can access those). I have
1) When I run the code, I get the following error: Warning message:
In lav_object_post_check(object) : lavaan WARNING: some estimated lv variances are negative. I see that the variance is negative for Group 2 for the variable representing the latent change score. If I split the groups apart and run the model separately, the warning message only emerges with Group 2. Does anyone have a suggestion for addressing this? I tried altering the factor loadings to different values, but did not find a solution there.
2) Assuming the model fits (and I realize at this point I don't have the best model fit yet, but more data is coming), what is the appropriate way to test/conclude whether the change score is significant within a group? And also that it might differ between the 2 groups, say Group 1 changes more (or has a significant change) whereas group 2 changes less (or is not significant).
3) Should I test for measurement invariance for this model, in this case between the 2 groups?
4) Finally, assuming the model fits appropriately again, is there any other output from this model that will help me assess change within and/or across the 2 groups? Should I look at the variances and/or intercepts, for instances, of my change score and, if so, how do I compare the result across the 2 groups.
Thank you.
Chris
mod1<-'
pre.wm=~1*RTotal1 + 1*STotal1 # measurement model at pre
post.wm=~1*RTotal2+equal(pre.wm=~STotal1)*STotal2 #measurement model at post, with the equality constrained factor loadings
post.wm ~ 1*pre.wm # Fixed regression of post on pre
dCOG1 =~ 1*post.wm # Fixed regression of dCOG1 on post
post.wm ~ 0*1 # constrains the intercept of post to 0
post.wm ~~ 0*post.wm # fixes the variance of post to 0
dCOG1 ~ 1 # This estimates the intercept of the change score
pre.wm ~ 1 # This estimates the intercept of pre
dCOG1 ~~ dCOG1 # This estimates the variance of the change scores
pre.wm ~~ pre.wm # This estimates the variance of the pre
dCOG1~pre.wm + Age + Sex + Fitness # This estimates the self-feedback parameter, with covariates
RTotal1~~RTotal2 # This allows residual covariance on indicator X1 across T1 and T2
STotal1~~STotal2 # This allows residual covariance on indicator X2 across T1 and T2
RTotal1~~RTotal1 # This allows residual variance on indicator X1
STotal1~~STotal1 # This allows residual variance on indicator X2
RTotal2~~equal("RTotal1~~RTotal1")*RTotal2 # This allows residual variance on indicator X1 at T2
STotal2~~equal("STotal1~~STotal1")*STotal2 # This allows residual variance on indicator X2 at T2
RTotal1~0*1 # This constrains the intercept of X1 to 0 at T1
STotal1~1 # This estimates the intercept of X2 at T1
RTotal2~0*1 # This constrains the intercept of X1 to 0 at T2
STotal2~equal("STotal1~1")*1 # This estimates the intercept of X2 at T2
'
fitmod1 <- lavaan(mod1, data=all.3, group='Group' , estimator='mlr',fixed.x=T,missing='fiml') #fixed.x=T regards all exogenuous covariates as estimated
summary(fitmod1, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
#And my output:
lavaan (0.5-23.1097) converged normally after 66 iterations
Number of observations per group
1 42
2 45
Number of missing patterns per group
1 2
2 2
Estimator ML Robust
Minimum Function Test Statistic 37.721 39.257
Degrees of freedom 21 21
P-value (Chi-square) 0.014 0.009
Scaling correction factor 0.961
for the Yuan-Bentler correction
Chi-square for each group:
1 21.038 21.894
2 16.684 17.363
Model test baseline model:
Minimum Function Test Statistic 155.747 151.404
Degrees of freedom 36 36
P-value 0.000 0.000
User model versus baseline model:
Comparative Fit Index (CFI) 0.860 0.842
Tucker-Lewis Index (TLI) 0.761 0.729
Robust Comparative Fit Index (CFI) 0.852
Robust Tucker-Lewis Index (TLI) 0.747
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1007.898 -1007.898
Scaling correction factor 0.913
for the MLR correction
Loglikelihood unrestricted model (H1) -989.037 -989.037
Scaling correction factor 0.985
for the MLR correction
Number of free parameters 31 31
Akaike (AIC) 2077.795 2077.795
Bayesian (BIC) 2154.238 2154.238
Sample-size adjusted Bayesian (BIC) 2056.423 2056.423
Root Mean Square Error of Approximation:
RMSEA 0.135 0.141
90 Percent Confidence Interval 0.060 0.204 0.068 0.210
P-value RMSEA <= 0.05 0.035 0.027
Robust RMSEA 0.139
90 Percent Confidence Interval 0.068 0.205
Standardized Root Mean Square Residual:
SRMR 0.123 0.123
Parameter Estimates:
Information Observed
Standard Errors Robust.huber.white
Group 1 [1]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
pre.wm =~
RTotal1 1.000 0.744 0.677
STotal1 1.000 0.744 0.858
post.wm =~
RTotal2 1.000 0.576 0.580
STotal2 (STt1) 1.521 0.340 4.475 0.000 0.876 0.892
dCOG1 =~
post.wm 1.000 0.688 0.688
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
post.wm ~
pre.wm 1.000 1.292 1.292
dCOG1 ~
pre.wm -0.342 0.154 -2.220 0.026 -0.643 -0.643
Age -0.014 0.021 -0.691 0.489 -0.036 -0.187
Sex -0.168 0.177 -0.950 0.342 -0.424 -0.186
Fitness -0.013 0.017 -0.776 0.437 -0.033 -0.199
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.RTotal1 ~~
.RTotal2 0.351 0.138 2.540 0.011 0.351 0.537
.STotal1 ~~
.STotal2 -0.028 0.160 -0.175 0.861 -0.028 -0.142
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.post.wm 0.000 0.000 0.000
.dCOG1 0.884 0.729 1.213 0.225 2.229 2.229
pre.wm 0.000 0.159 0.001 0.999 0.000 0.000
.RTotal1 0.000 0.000 0.000
.STotal1 (.24.) -0.010 0.130 -0.074 0.941 -0.010 -0.011
.RTotal2 0.000 0.000 0.000
.STotal2 (ST1~) -0.010 0.130 -0.074 0.941 -0.010 -0.010
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.post.wm 0.000 0.000 0.000
.dCOG1 0.078 0.046 1.703 0.089 0.496 0.496
pre.wm 0.554 0.166 3.329 0.001 1.000 1.000
.RTotal1 (.19.) 0.654 0.139 4.701 0.000 0.654 0.541
.STotal1 (.20.) 0.198 0.148 1.340 0.180 0.198 0.263
.RTotal2 (RT1~) 0.654 0.139 4.701 0.000 0.654 0.663
.STotal2 (ST1~) 0.198 0.148 1.340 0.180 0.198 0.205
Group 2 [2]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
pre.wm =~
RTotal1 1.000 0.687 0.709
STotal1 1.000 0.687 0.631
post.wm =~
RTotal2 1.000 0.366 0.404
STotal2 1.866 0.896 2.082 0.037 0.684 0.760
dCOG1 =~
post.wm 1.000 0.772 0.772
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
post.wm ~
pre.wm 1.000 1.874 1.874
dCOG1 ~
pre.wm -0.443 0.287 -1.542 0.123 -1.074 -1.074
Age -0.020 0.012 -1.648 0.099 -0.071 -0.363
Sex 0.099 0.120 0.823 0.410 0.349 0.162
Fitness -0.005 0.009 -0.608 0.543 -0.019 -0.096
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.RTotal1 ~~
.RTotal2 -0.032 0.145 -0.218 0.827 -0.032 -0.056
.STotal1 ~~
.STotal2 0.004 0.120 0.037 0.970 0.004 0.009
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.post.wm 0.000 0.000 0.000
.dCOG1 0.322 0.475 0.678 0.498 1.139 1.139
pre.wm 0.007 0.142 0.051 0.959 0.011 0.011
.RTotal1 0.000 0.000 0.000
.STotal1 -0.013 0.162 -0.078 0.938 -0.013 -0.012
.RTotal2 0.000 0.000 0.000
.STotal2 0.151 0.348 0.433 0.665 0.151 0.168
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.post.wm 0.000 0.000 0.000
.dCOG1 -0.027 0.087 -0.306 0.759 -0.334 -0.334
pre.wm 0.471 0.160 2.951 0.003 1.000 1.000
.RTotal1 0.465 0.172 2.710 0.007 0.465 0.497
.STotal1 0.712 0.180 3.962 0.000 0.712 0.602
.RTotal2 0.687 0.153 4.488 0.000 0.687 0.836
.STotal2 0.342 0.255 1.341 0.180 0.342 0.422
Terrence Jorgensen
Jul 25, 2017, 6:23:35 AM7/25/17
to lavaan
1) the variance is negative for Group 2 for the variable representing the latent change score. Does anyone have a suggestion for addressing this?
Because your sample size is so small, an out-of-bounds correlation might occur frequently just due to sampling error. Read Bollen & Kolenikov (2012) about a method to use bootstrap confidence intervals to test whether an out-of-bounds estimate is out-of-bounds in the population:
Given the relative size of your SE for the negative variance, I don't think you can dismiss sampling error as the explanation. But given the poor model fit, you also cannot dismiss model misspecification as the explanation.
2) what is the appropriate way to test/conclude whether the change score is significant within a group?
There is a Wald z test associated with the intercept of the change score
And also that it might differ between the 2 groups, say Group 1 changes more (or has a significant change) whereas group 2 changes less (or is not significant).
You can either constrain that parameter to equality across groups and test the H0 of equality by comparing the nested models with a LRT using anova() or lavTestLRT(), or you can calculate the difference as a user-defined parameter in the model syntax by following the instructions on the mediation tutorial:
The user-defined parameter will also have a Wald z test
3) Should I test for measurement invariance for this model, in this case between the 2 groups?
You don't have enough indicators per factor for that test to be possible. You just have to assume it, which you do by fixing the loadings to 1.
4) Finally, assuming the model fits appropriately again, is there any other output from this model that will help me assess change within and/or across the 2 groups? Should I look at the variances and/or intercepts, for instances, of my change score and, if so, how do I compare the result across the 2 groups.
Depends what hypotheses you want to test.
pre.wm=~1*RTotal1 + 1*STotal1 # measurement model at prepost.wm=~1*RTotal2+equal(pre.wm=~STotal1)*STotal2 #measurement model at post, with the equality constrained factor loadings
Why are you using the equal() operator to set it equal to an already fixed value? Just fix it to 1 like you did the other 3 loadings. The equal() operator is for constraining estimated parameters to equality, not fixed parameters.
RTotal2~~equal("RTotal1~~RTotal1")*RTotal2 # This allows residual variance on indicator X1 at T2STotal2~~equal("STotal1~~STotal1")*STotal2 # This allows residual variance on indicator X2 at T2
Why are you forcing the residual variances to equality across time? That's a testable assumption (if you find it interesting or simply want to make your model more parsimonious given the small sample), but I would not recommend starting out with that assumption.
STotal2~equal("STotal1~1")*1 # This estimates the intercept of X2 at T2
Because you have multiple groups, you should label the parameters in the syntax, so that you can provide a unique label in each group.
That will also make the equal() operator unnecessary because you can simply apply the same label(s) when you want the estimates to be constrained to equality.
fitmod1 <- lavaan(mod1, data=all.3, group='Group' , estimator='mlr',fixed.x=T,missing='fiml') #fixed.x=T regards all exogenuous covariates as estimated
No, that would be fixed.x = FALSE (i.e., do not treat them as fixed, estimate them instead). If you have any missing data in your exogenous covariates, you want to set this to FALSE.
Terrence D. Jorgensen
Postdoctoral Researcher, Methods and Statistics
Research Institute for Child Development and Education, the University of Amsterdam
UvA web page: http://www.uva.nl/profile/t.d.jorgensen
cezre...@gmail.com
Jul 27, 2017, 10:54:39 AM7/27/17
to lavaan
Thanks very much, Terrence, for your helpful responses. A few follow-ups, if I may (which are in red below).
On Tuesday, July 25, 2017 at 5:23:35 AM UTC-5, Terrence Jorgensen wrote:
1) the variance is negative for Group 2 for the variable representing the latent change score. Does anyone have a suggestion for addressing this?Because your sample size is so small, an out-of-bounds correlation might occur frequently just due to sampling error. Read Bollen & Kolenikov (2012) about a method to use bootstrap confidence intervals to test whether an out-of-bounds estimate is out-of-bounds in the population:
I looked at Bollen & Kolenikov (2012). Is there a particular reason you recommend bootstrap CI's over other methods suggested in the paper?
Given the relative size of your SE for the negative variance, I don't think you can dismiss sampling error as the explanation. But given the poor model fit, you also cannot dismiss model misspecification as the explanation.2) what is the appropriate way to test/conclude whether the change score is significant within a group?There is a Wald z test associated with the intercept of the change score
Ah! So I could use lavTestWald() to test the change score intercept. I did not know about this function. I was thinking the z-score of the change score intercept in the summary() output might be sufficient. What is the difference between the z-score from the summary output versus the Wald z-test?
And also that it might differ between the 2 groups, say Group 1 changes more (or has a significant change) whereas group 2 changes less (or is not significant).You can either constrain that parameter to equality across groups and test the H0 of equality by comparing the nested models with a LRT using anova() or lavTestLRT(), or you can calculate the difference as a user-defined parameter in the model syntax by following the instructions on the mediation tutorial:The user-defined parameter will also have a Wald z test
Thanks--makes sense. If I use the anova() approach to compare two models, the result is like an omnibus anova which would tell if there were differences and, for 2 groups, it would be straightforward to infer if there was a difference. But what if I have more than 2 groups? How would the anova() or lavTestLRT() tell me which groups differ significantly?
3) Should I test for measurement invariance for this model, in this case between the 2 groups?You don't have enough indicators per factor for that test to be possible. You just have to assume it, which you do by fixing the loadings to 1.
Okay--thanks. That makes sense.
4) Finally, assuming the model fits appropriately again, is there any other output from this model that will help me assess change within and/or across the 2 groups? Should I look at the variances and/or intercepts, for instances, of my change score and, if so, how do I compare the result across the 2 groups.Depends what hypotheses you want to test.pre.wm=~1*RTotal1 + 1*STotal1 # measurement model at prepost.wm=~1*RTotal2+equal(pre.wm=~STotal1)*STotal2 #measurement model at post, with the equality constrained factor loadingsWhy are you using the equal() operator to set it equal to an already fixed value? Just fix it to 1 like you did the other 3 loadings. The equal() operator is for constraining estimated parameters to equality, not fixed parameters.
Right. The equal() operator is a relic from earlier incarnations of the code when I was not using 1. But yes, that is an overly complicated way to specify loadings when they are all 1.
RTotal2~~equal("RTotal1~~RTotal1")*RTotal2 # This allows residual variance on indicator X1 at T2STotal2~~equal("STotal1~~STotal1")*STotal2 # This allows residual variance on indicator X2 at T2Why are you forcing the residual variances to equality across time? That's a testable assumption (if you find it interesting or simply want to make your model more parsimonious given the small sample), but I would not recommend starting out with that assumption.
Good point.
STotal2~equal("STotal1~1")*1 # This estimates the intercept of X2 at T2Because you have multiple groups, you should label the parameters in the syntax, so that you can provide a unique label in each group.
This is a good idea; but here X1 and X2 denote the indicators (not groups). Right now, I specify group in the lavaan() call, to provide output for each group. Is there a different/better way to set up my code, so as to have group level output specified within the syntax? I am interested in group comparisons sometimes (as discussed above) and, with the model set up this way, I have to resort to a model comparison approach (with anova() for instance) to compare group means.
That will also make the equal() operator unnecessary because you can simply apply the same label(s) when you want the estimates to be constrained to equality.
Nice. Did not know this.
fitmod1 <- lavaan(mod1, data=all.3, group='Group' , estimator='mlr',fixed.x=T,missing='fiml') #fixed.x=T regards all exogenuous covariates as estimatedNo, that would be fixed.x = FALSE (i.e., do not treat them as fixed, estimate them instead). If you have any missing data in your exogenous covariates, you want to set this to FALSE.
Oh--thanks. I misunderstood this call for fixed.x.
Overall, thank you. I am just getting into SEM and I don't yet understand all the nuances yet so this has been quite enlightening.
Terrence Jorgensen
Jul 28, 2017, 6:34:43 AM7/28/17
to lavaan
I looked at Bollen & Kolenikov (2012). Is there a particular reason you recommend bootstrap CI's over other methods suggested in the paper?
I don't, necessarily. Robust CIs work fine too, so you can request se = "robust.sem" or some robust estimator and refer to those. There is a table in the paper that shows what is appropriate under different circumstances, and robust and bootstrap CIs appear in the same box when normality is violated (and I think when the model is not perfectly specified).
Ah! So I could use lavTestWald() to test the change score intercept. I did not know about this function. I was thinking the z-score of the change score intercept in the summary() output might be sufficient. What is the difference between the z-score from the summary output versus the Wald z-test?
I was referring to the Wald z statistic in the summary() output. But if you use lavTestWald() to "constrain" that parameter to zero (without fitting a new model), you will get an identical z statistic and p value.
Thanks--makes sense. If I use the anova() approach to compare two models, the result is like an omnibus anova which would tell if there were differences and, for 2 groups, it would be straightforward to infer if there was a difference. But what if I have more than 2 groups? How would the anova() or lavTestLRT() tell me which groups differ significantly?
No, as you said, it would be an omnibus test. So if there are 3+ groups, you would need to conduct pairwise follow-up tests to see which groups differ, like post-hoc pairwise comparisons of means following a significant F in ANOVA.
Is there a different/better way to set up my code, so as to have group level output specified within the syntax?
As the tutorial above explains, you can add a vector of labels (one per group). So in the simplest 2-group case, you can add different labels to get unique estimates in each group
STotal1 ~ c(int1.g1, int1.g2)*1
STotal2 ~ c(int2.g1, int2.g2)*1To constrain across groups, give each group the same label and do a model comparison with anova()
STotal1 ~ c(int1, int1)*1
STotal2 ~ c(int2, int2)*1Or you can use the same model syntax and specify an equality constraint using the constraints argument when you run the model again.
out1 <- lavaan(syntax, data)
out2 <- lavaan(syntax, data, constraints = c("int1.g1 == int1.g2","int2.g1 == int2.g2"))Or, without running separate models (which is necessary for an omnibus test of all simultaneous constraints), you could test each pair or parameters for equality by defining difference parameters in the unconstrained model's syntax:
int1.diff := int1.g1 - int1.g2
int2.diff := int2.g1 - int2.g2Overall, thank you. I am just getting into SEM and I don't yet understand all the nuances yet so this has been quite enlightening.
You're welcome, good luck.
Reply all
Reply to author
Forward
0 new messages