Question regarding analysis proceeding and corresponding lavaan commands

[Isa]

Dear Lavaan Group

 

I would like to test Mediation Effects described in the following MIMIC Model Structure (see attachement). As illustrated by the diagram, the Mediator is a Latent Variable (I omitted its ordinal indicators to keep the diagram simple). Notice that most of the exogenous variables/covariates are binary.

 

As I learned in a SEM introduction course, I first would like to fit the Measurement Part of my Model, which in my case is the single Mediator Latent Variable and its corresponding indicators. I have done this using the cfa() command in lavaan and the WLSMV estimator.

 

Now, I would like to fit the Structural Part of the model, meaning the MIMIC Model itself in order to assess associations between the Latent Variables and covariates (Mediation).  To do so, I used the sem() command in lavaan and again the WLSMV estimator. Is this correct up to this point?

 

In a last step, I would like to explore group differences in my fitted MIMIC Model for the two groups for aetiology: traumatic vs. non-traumatic, dummy coded in variable sci_cause  (I tried to follow this article: https://www.frontiersin.org/articles/10.3389/fpsyg.2018.00740/full). Therefore, I first tested my Measurement Model for Measurement Invariance (MI) in the aetiology variable. I used the commands:

 

config_eti <- cfa(cfa_se,  estimator = "WLSMV", data=data_imputed, group="sci_cause", parameterization = "theta")

 

weak_eti <- cfa(cfa_se,  estimator = "WLSMV", data=data_imputed, group="sci_cause", group.equal=c("loadings"), parameterization = "theta")

lavTestLRT(config_eti, weak_eti, method="satorra.2000")

 

strong_ eti <- cfa(cfa_se,  estimator = "WLSMV", data=data_imputed, group=" sci_cause", group.equal=c("loadings","thresholds"), parameterization = "theta")

lavTestLRT(weak_eti, strong_ eti, method="satorra.2000")

 

In my Model, MI holds on the level of factor loadings. I am now allowed to asses Group Differences in the MIMIC Model by: 0) Fitting again the MIMIC Model without any constrains on regression parameters by the sem() command, and  1) Fitting the MIMIC Model separately for both groups and setting the regression parameters equal across groups by using the command sem(…, group = "sci_cause"). Difference in Model Fit for the Models 0 and 1 can now be assessed by the command:

 

lavTestLRT(sem_fit_mod_se_eti_0,  sem_fit_mod_se_eti_1,  method="satorra.2000")

 

(significant chi-squared difference would indicate group differences in my MIMIC Model on the level of this Fit Index)

 

Are this proceeding and the corresponding lavaan commands correct/reasonable so far? I am unsure about the MI test: should I do it for the Measurement Model only, or also for the MIMIC Model?

 

Thanks for any recommendations and comments!

Best regards, isabel

[Terrence Jorgensen]

I used the sem() command in lavaan and again the WLSMV estimator. Is this correct up to this point?

Yes.  FYI, the cfa() and sem() functions are identical.  They are both wrappers around the general lavaan() function that turn on all the same default options, so they will always give you the same result.

In a last step, I would like to explore group differences in my fitted MIMIC Model for the two groups for aetiology: traumatic vs. non-traumatic, dummy coded in variable sci_cause  

This is confusing, because it makes it sound like you want to use sci_cause as a predictor (cause) in your model.  But you are not.  Instead, you are fitting a multigroup model, which does not explicitly use sci_cause as a predictor, but merely allows (some or all) parameter estimates to differ between sci_cause groups.  This is still a "MIMIC model" in the sense that you have exogenous observed predictors of the latent mediator.  But the MIMIC quality just seems like a distraction from your actual questions, which are about testing invariance of the measurement parameters across sci_cause groups.  (Also, the article you referenced does not mention anything about MIMIC models.)

I am now allowed to asses Group Differences in the MIMIC Model by: 0) Fitting again the MIMIC Model without any constrains on regression parameters by the sem() command, and  1) Fitting the MIMIC Model separately for both groups and setting the regression parameters equal across groups by using the command sem(…, group = "sci_cause"). Difference in Model Fit for the Models 0 and 1 can now be assessed by the command:

 

lavTestLRT(sem_fit_mod_se_eti_0,  sem_fit_mod_se_eti_1,  method="satorra.2000")

 

(significant chi-squared difference would indicate group differences in my MIMIC Model on the level of this Fit Index)

 

Are this proceeding and the corresponding lavaan commands correct/reasonable so far?

Yes.

I am unsure about the MI test: should I do it for the Measurement Model only, or also for the MIMIC Model?

The whole model is a mimic model, so please refer to specific parameters you are curious about.  Testing measurement invariance only involves the equality of the mediator's factor loadings across sci_cause groups.  (You do not need to establish strong invariance because your questions are all about covariance structure, not means/intercepts.)  Given that you have weak invariance, you can next constrain slopes pointing to and from the mediator to equality across groups, and compare that model to the weak-invariance model to see if those effects are moderated by sci_cause.  In my mind, those are the "MIMIC model" parameters, distinct from the measurement-model parameters, but I'm not sure if that is what you meant.

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Isa]

Dear Terrence

Thank you for your critical reading and for the clarification of some points!  

Please find in the attachment another picture in which I wanted to show what I meant by Measurement Model and MIMIC Model and which analyses I planned to do with which part of the model . (Should be consistent with the description in my first post)

In a last step, I would like to explore group differences in my fitted MIMIC Model for the two groups for aetiology: traumatic vs. non-traumatic, dummy coded in variable sci_cause  

This is confusing, because it makes it sound like you want to use sci_cause as a predictor (cause) in your model.  But you are not.  Instead, you are fitting a multigroup model, which does not explicitly use sci_cause as a predictor, but merely allows (some or all) parameter estimates to differ between sci_cause groups.  This is still a "MIMIC model" in the sense that you have exogenous observed predictors of the latent mediator.  But the MIMIC quality just seems like a distraction from your actual questions, which are about testing invariance of the measurement parameters across sci_cause groups.  (Also, the article you referenced does not mention anything about MIMIC models.)

Yes, you are right, the MIMIC Quality is distracting from the question of Measurement Invariance Testing.

 

My idea was the following:

I would like to assess Measurement Invariance for sci_cause in my Measurement Model (the Mediator Construct only) using Multiple Group CFA:

 

cfa_se <- "SE =~ item1 + item2 + item3 + item4 + item5 + item6"

config_eti <- cfa(cfa_se,  estimator = "WLSMV", data=data_imputed, group="sci_cause", parameterization = "theta")

weak_eti <- cfa(cfa_se,  estimator = "WLSMV", data=data_imputed, group="sci_cause", group.equal=c("loadings"), parameterization = "theta")

lavTestLRT(config_eti, weak_eti, method="satorra.2000")

And then I would be able to test Group Differences for sci_cause in the bigger Structural Model (the MIMIC Model) by doing:

 

sem_mod_se_eti_0 <- "

# Latent constructs

SE =~ item1 + item2 + item3 + item4 + item5 + item6

# Regression

# Direct effect

 scim_location ~ decubitus + urin_infection  + pulmonary + cardiac_function + bowel_care + pain + sex + age_sci_group + partet + compl

# Mediation

SE ~ decubitus + urin_infection  + pulmonary + cardiac_function + bowel_care + pain + partet + compl

scim_location ~ SE

"

sem_fit_mod_se_eti_0 <- sem(sem_mod_se_eti_0,  estimator = "WLSMV", data = data_imputed, group = "sci_cause")

 

sem_mod_se_eti_1 <- "

# Latent constructs

SE =~ item1 + item2 + item3 + item4 + item5 + item6

# Regression

# Direct effect

scim_location ~ c(a,a)*decubitus + c(b,b)*urin_infection  + c(c,c)*pulmonary + c(d,d)*cardiac_function + c(e,e)*bowel_care + c(f,f)*pain + c(g,g)*sex + c(h,h)*age_sci_group + c(i,i)*partet + c(j,j)*compl

# Mediation

SE ~ c(z,z)*decubitus + c(y,y)*urin_infection  + c(x,x)*pulmonary + c(w,w)*cardiac_function + c(v,v)*bowel_care + c(u,u)*pain + c(t,t)*partet + c(s,s)*compl

scim_location ~ c(q,q)*SE

"

sem_fit_mod_se_eti_1 <- sem(sem_mod_se_eti_1,  estimator = "WLSMV", data = data_imputed, group = "sci_cause")

 

lavTestLRT(sem_fit_mod_se_eti_0, sem_fit_mod_se_eti_1, method="satorra.2000")

 

I hope this is now a better description of what I would like to do. Do you think this proceeding is correct?

 

I am unsure about the MI test: should I do it for the Measurement Model only, or also for the MIMIC Model?

The whole model is a mimic model, so please refer to specific parameters you are curious about.  Testing measurement invariance only involves the equality of the mediator's factor loadings across sci_cause groups.  (You do not need to establish strong invariance because your questions are all about covariance structure, not means/intercepts.)  Given that you have weak invariance, you can next constrain slopes pointing to and from the mediator to equality across groups, and compare that model to the weak-invariance model to see if those effects are moderated by sci_cause. 

I am not sure if I understood this part correctly. Two questions:

1. What do you mean by “weak-invariance model”: In the code above, would this correspond to the model sem_mod_se_eti_0?

2. As far as I planned, I test the whole Structural Model as such for Group Differences only, meaning in a first step I don’t impose any constrains regarding the regression parameters of my Model and in the second step I constrain all corresponding regression parameters to be the same across the two Models for the two groups of sci_cause (see code above). Regarding your command, it seems that I also could test each Mediation Path separately for Group Differences, is this correct?

 

Thanks again and best wishes,

Isabel

[Terrence Jorgensen]

Do you think this proceeding is correct?

Yes.

1. What do you mean by “weak-invariance model”:

The model with equality constraints only on loadings.

In the code above, would this correspond to the model sem_mod_se_eti_0?

Yes.

2. As far as I planned, I test the whole Structural Model as such for Group Differences only, meaning in a first step I don’t impose any constrains regarding the regression parameters of my Model and in the second step I constrain all corresponding regression parameters to be the same across the two Models for the two groups of sci_cause (see code above). Regarding your command, it seems that I also could test each Mediation Path separately for Group Differences, is this correct?

Yes, if you reject the null-hypothesized sem_mod_se_eti_1, then groups have different regression slopes in sem_mod_se_eti_1 .  So you would provide labels for the parameters, like in sem_mod_se_eti_1, except that the labels differ across groups, e.g.:  SE ~ c(z1, z2)*decubitus

Then you define indirect and total paths the same way, but separate parameters per group.  

via_decubitus_1 := z1*q1
via_decubitus_2 := z2*q2

You could also define group differences between the same parameters, to test whether each individual (direct or indirect) effect 

via_decubitus_diff :=  via_decubitus_1  - via_decubitus_2

[Isa]

Dear Terrence,

thank you very much for your support and this nice idea to test same parameters for the different groups!

May I ask you another questions that appeared? I would like to show in a summary of my analysis also the correlation matrix of my data frame (all the variables used in the model).

I used the following command:

cor_full <- lavCor(data_imputed[,c("t4_scim_location",
  "t4_decubitus","t4_urin_infection","t4_pulmonary","t4_cardiac_function","t4_bowel_care","t4_pain",
  "age_sci_group","sex","t4_partet","t4_compl",
  "t4_belief_qualities_trait","t4_belief_useless_trait", "t4_belief_useful_trait", "t4_belief_positive_trait")])

And I get the following warning message:

Error in lav_data_full(data = data, group = group, cluster = cluster,  : 
  lavaan ERROR: unordered factor(s) detected; make them numeric or ordered: t4_decubitus t4_urin_infection t4_pulmonary t4_cardiac_function t4_bowel_care t4_pain age_sci_group sex t4_partet t4_compl

As you can see in the picture of my model which I posted some weeks ago, I am using also binary varibles for the analysis. I guess the command lavCor() is not able to calculate tetrachoric correlations?

Is there another option to get the correlation matrix?

Best wishes and many thanks again!

Isabel

[Terrence Jorgensen]

lavaan ERROR: unordered factor(s) detected; make them numeric or ordered: t4_decubitus t4_urin_infection t4_pulmonary t4_cardiac_function t4_bowel_care t4_pain age_sci_group sex t4_partet t4_compl

As you can see in the picture of my model which I posted some weeks ago, I am using also binary varibles for the analysis. I guess the command lavCor() is not able to calculate tetrachoric correlations?

It can, just follow the instructions in the message.  Binary variables should also be declared as ordered, not nominal factors.  (For binary data, the distinction between nominal and ordinal is arbitrary.)

Is there another option to get the correlation matrix?

Yes, it is already available from a fitted model.

lavInspect(fit, "cor.ov")

[Terrence Jorgensen]

Is there another option to get the correlation matrix?

Yes, it is already available from a fitted model.

lavInspect(fit, "cor.ov")

 

Sorry, those are model-implied correlations, so they will differ across models.  For the unrestricted polychoric correlations, use

lavInspect(fit, "sampstat")

[Isa]

ahhh, now it works for the correlation matrix. Thank you!

But your answer that the distinction for nominal and ordered is arbitrary for binary data raises another question:

In my model, I declared my binary variables (t4_decubitus,  t4_urin_infection, t4_pulmonary, t4_cardiac_function, t4_bowel_care, t4_pain, sex, t4_partet, t4_compl) as factors and the model worked without any warning.

When I declare them now as ordered and I fit the model, then I get the following warning:

sem_med_se <- "
# Latent construct
SE =~ t4_belief_qualities_trait + t4_belief_useless_trait + t4_belief_useful_trait + t4_belief_positive_trait

# Regression
# Direct effect

t4_scim_location ~ c1*t4_decubitus + c2*t4_urin_infection  + c3*t4_pulmonary + c4*t4_cardiac_function + c5*t4_bowel_care + c6*t4_pain
+ age_sci_group + sex + c7*t4_partet + c8*t4_compl

# Mediation
SE ~ a1*t4_decubitus + a2*t4_urin_infection  + a3*t4_pulmonary + a4*t4_cardiac_function + a5*t4_bowel_care + a6*t4_pain
+ a7*t4_partet + a8*t4_compl

t4_scim_location ~ b*SE

# Indirect effect (a*b)
a1b := a1*b
a2b := a2*b
a3b := a3*b
a4b := a4*b
a5b := a5*b
a6b := a6*b
a7b := a7*b
a8b := a8*b

# total effect
total1 := c1 + (a1*b)
total2 := c2 + (a2*b)
total3 := c3 + (a3*b)
total4 := c4 + (a4*b)
total5 := c5 + (a5*b)
total6 := c6 + (a6*b)
total7 := c7 + (a7*b)
total8 := c8 + (a8*b)
"
sem_fit_med_se <- sem(sem_med_se,  estimator = "WLSMV", data = data_imputed, mimic="Mplus")

Warning message:
In lav_data_full(data = data, group = group, cluster = cluster,  :
  lavaan WARNING: exogenous variable(s) declared as ordered in data: t4_decubitus t4_urin_infection t4_pulmonary t4_cardiac_function t4_bowel_care t4_pain age_sci_group sex t4_partet t4_compl

Do you know what is the issue here? I thought the model is working fine...

Best wishes, Isabel

[Terrence Jorgensen]

  lavaan WARNING: exogenous variable(s) declared as ordered in data: t4_decubitus t4_urin_infection t4_pulmonary t4_cardiac_function t4_bowel_care t4_pain age_sci_group sex t4_partet t4_compl

Do you know what is the issue here?

Exogenous categorical variables should only be represented with numeric codes.

http://lavaan.ugent.be/tutorial/cat.html

[Isa]

Thanks again a lot for your answer!

I now classified all my binary exogenous variables as "numeric" and made them dummy variables coded with 0 and 1.

The warning message disappeared. Just out of curiosity and to understand what is going on: When I check the fit and Regression estimates now, they are the same like when I classified the binary exogenous variables wrongly as "factors". Did lavaan automatically treat them as numerics?

Best regards, isabel

[Terrence Jorgensen]

Did lavaan automatically treat them as numerics?

Sounds like it does.  But the warning is probably to make sure users understand that, and to use fixed.x=TRUE so that nothing is assumed about their distribution (which would invalidate your test statistics).