Re: Latent Growth Curve Model

[Terrence Jorgensen]

1. Am I doing this right?

Looks right. 

1. How do I interpret the output?

https://lavaan.ugent.be/tutorial/growth.html

Try Googling for more tutorials about growth curves in lavaan, like this one: https://quantdev.ssri.psu.edu/tutorials/growth-modeling-chapter-10-growth-models-nonlinearity-time

1. Why are all intercepts 0?

They aren't.  Only the indicators have intercepts == 0 because their means are reproduces as functions of the growth parameters.  The growth parameter means are merely very close to zero (try adding the argument nd=7 for the summary() to print up to 7 digits of precision).  Slopes nearly zero seems very consistent with the dark blue line in your plot being almost flat.  The plot indicates your (unconditional) average intercept should be somewhere between 2 and 3; however,  you have several predictors of the growth factors.  Thus, the intercepts are the mean level and change when all those predictors == 0.  Does that make sense for your data?

1. What does it mean that intercept, slope, and quadratic slope are not significant?

It implies those parameters are consistent with a population parameter == 0.

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

 

[Manuel]

Dear Terrence,

Thank you very much for your answer. It helped a lot.

Is there a way to run a double latent growth model in Lavaan?

I would like to use boredom factor scores instead of manifest boredom means. In my mind it should look like this:

latent_boredom_growth <- '

boredom_0 =~ bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0; #latent boredom

boredom_1 =~ bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1;

boredom_2 =~ bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2;

    

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2; # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2; # slope

'

boredom_growth_model = growth(latent_boredom_growth, data=lavan)

But this results in an error:

In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    Could not compute standard errors! The information matrix could

    not be inverted. This may be a symptom that the model is not

    identified.

Thanks again and kind regards,

Manuel

[Mauricio Garnier-Villarreal]

Manuel

If you want to do a Growth curve with factors, you first have to establsih longitudinal strong measurement invariance for the Boredom_0, boredom_1, and boredom_2 factors.

Steps:

- Configural model: The factor model with no equality constrint over time. Using fixed variance identification method

- Weak invariance: Equate factor loadings across time, freely estimate the variance at time 2 and 3, keep the variance of time 1 fixed to 1

- Strong invariance: equate item intercepts over time, freely estimate the mean of time 2 and 3, keep the mean at time 1 fixed at 0

Between each step you compare the model to the previous model, to see if the constraints decrease fit substantially

After that you can use the factors from the Strong invariance model for the growth curve

[Patrick (Malone Quantitative)]

Manuel,

I suspect you're running into issues with the factor intercepts.

I'm not sure whether the growth() function can do this, but the sem()
function could. Then if you specify scalar invariance on the boredom
factors as Mauricio mentioned, then

boredom_0 ~ 0*1
boredom_1 ~ 0*1
boredom_2 ~ 0*1
i ~ 1
s ~ 1

That should get you on the right track.

> --
> You received this message because you are subscribed to the Google Groups "lavaan" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/2ac441be-e959-4a3b-b72e-f72997d94bebn%40googlegroups.com.



--
Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

[Manuel]

Dear Magua, dear Patrick,

Thank you for your help!

I have managed to establish scalar longitudinal measurement invariance for boredom and added the scalar boredom factor scores to my data frame using this code:

idx <- lavInspect(fit_strong_boredom, "case.idx")

fscores <- lavPredict(fit_strong_boredom, assemble = TRUE)

## loop over factors

for (fs in colnames(fscores)) {

  df1[idx, fs] <- fscores[, fs]

}

Now I am trying to extract the slope parameters using the boredom factor scores:

boredom_growth <- '

    i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

    s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope'

growthCurveModel = growth(boredom_growth, data=df1)

But this again results in an error:

Error in lav_data_full(data = data, group = group, cluster = cluster,  : 

  lavaan ERROR: unordered factor(s) detected; make them numeric or ordered: boredom_0 boredom_1 boredom_2

The factor scores range from -0,000154379130794116 to 2,78171374926366. Using as.numeric(df1$boredom_0) strangely results in values from 1 to 1186.

Do you have any idea how I can solve this?

Thanks again,

Manuel

[Patrick (Malone Quantitative)]

Manuel,

It sounds like a different usage of "factor."

What do you get from class(df1$boredom_0) ? It should be either numeric (double) or ordered, not factor. Factor in that sense would mean a nominal variable, which can't be used as an endogenous variable in lavaan. Using as.numeric will give a value for each unique value of a nominal variable.

lavPredict() should be returning class numeric. I'm not sure where they're getting converted to class factor.

Pat

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/5f9bb835-84c4-4390-b659-33e6fef4d105n%40googlegroups.com.

[Mauricio Garnier-Villarreal]

Manuel

It is not recommended for you to add the factor scores as measured variables, and then run the growth model with them. This goes back into the porblem of factor score indeteminancy and other issues of factor score variability and method.

The recommend steps would be t estimate the growth curve on top of the longitudinal scalar invariance model. This would required you to add 2 pieces to the measurement inavriance model

- The intercept and slope factors define by the first order factors, as you have it

- fix the first order factors means to 0. This way the mean change over time would be only define by the growth curve

Hope this helps

[Manuel]

Dear Mauricio,

I see. Then I misunderstood you saying I should "use the factors from the Strong invariance model for the growth curve". It also felt unusual.

This is my model for longitudinal scalar invariance and I added the growth curve at the end:

strong_boredom_growth <- '

 # measurement model

 boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

 boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

 boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

 # intercepts

 bos1_0 ~ i1*1

 bos2_0 ~ i2*1

 bos3_0 ~ i3*1

 bos4_0 ~ i4*1

 bos5_0 ~ i5*1

 bos6_0 ~ i6*1

 bos1_1 ~ i1*1

 bos2_1 ~ i2*1

 bos3_1 ~ i3*1

 bos4_1 ~ i4*1

 bos5_1 ~ i5*1

 bos6_1 ~ i6*1

 bos1_2 ~ i1*1

 bos2_2 ~ i2*1

 bos3_2 ~ i3*1

 bos4_2 ~ i4*1

 bos5_2 ~ i5*1

 bos6_2 ~ i6*1

 # unique variances and covariances

 bos1_0 ~~ bos1_0

 bos2_0 ~~ bos2_0

 bos3_0 ~~ bos3_0

 bos4_0 ~~ bos4_0

 bos5_0 ~~ bos5_0

 bos6_0 ~~ bos6_0

 bos1_1 ~~ bos1_1

 bos2_1 ~~ bos2_1

 bos3_1 ~~ bos3_1

 bos4_1 ~~ bos4_1

 bos5_1 ~~ bos5_1

 bos6_1 ~~ bos6_1

 bos1_2 ~~ bos1_2

 bos2_2 ~~ bos2_2

 bos3_2 ~~ bos3_2

 bos4_2 ~~ bos4_2

 bos5_2 ~~ bos5_2

 bos6_2 ~~ bos6_2

 bos1_0 ~~ bos1_1 + bos1_2

 bos2_0 ~~ bos2_1 + bos2_2

 bos3_0 ~~ bos3_1 + bos3_2

 bos4_0 ~~ bos4_1 + bos4_2

 bos5_0 ~~ bos5_1 + bos5_2

 bos6_0 ~~ bos6_1 + bos6_2 

 # latent variables means

 boredom_0 ~ 0*1

 boredom_1 ~ 1

 boredom_2 ~ 1

 # latent variables variances and covariances

 boredom_0 ~~ 1*boredom_0

 boredom_1 ~~ boredom_1

 boredom_2 ~~ boredom_2

 boredom_0 ~~ boredom_1

 boredom_0 ~~ boredom_2

 boredom_1 ~~ boredom_2

 # adding the growth model

 i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

 s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

'

# model estimation

boredom_growth_model = growth(strong_boredom_growth, data=df1)

summary(boredom_growth_model)

Unfortunately, this again results in an error:

1: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    The variance-covariance matrix of the estimated parameters (vcov)

    does not appear to be positive definite! The smallest eigenvalue

    (= 1.853961e-18) is close to zero. This may be a symptom that the

    model is not identified.

2: In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

I looked at the summary anyways and attached it as pdf, if this is of any help.

Kind regards,

Manuel

[Mauricio Garnier-Villarreal]

Manuel

Change these ines of code

You want the first order factor means fix to 0, so the mean change is represented by the intercepte and slope

Want to have the facto relation represented by the growth curve

 # latent variables means

 boredom_0 ~ 0*1

 boredom_1 ~ 0*1

 boredom_2 ~ 0*1

 # latent variables variances and covariances

 boredom_0 ~~ 1*boredom_0

 boredom_1 ~~ boredom_1

 boredom_2 ~~ boredom_2

 boredom_0 ~~ 0*boredom_1

 boredom_0 ~~ 0*boredom_2

 boredom_1 ~~ 0*boredom_2

[Manuel]

Unfortunately changing means to 0 still leads resulted in the error. Adding missing = "FIML" solved it though.

Just to double check I am doing it right, here is my code:

# 2. Define Models --------------------------------------------------------

# exploring measurement invariance over time

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# unique variances

bos1_0 ~~ bos1_0

bos2_0 ~~ bos2_0

bos3_0 ~~ bos3_0

bos4_0 ~~ bos4_0

bos5_0 ~~ bos5_0

bos6_0 ~~ bos6_0

bos1_1 ~~ bos1_1

bos2_1 ~~ bos2_1

bos3_1 ~~ bos3_1

bos4_1 ~~ bos4_1

bos5_1 ~~ bos5_1

bos6_1 ~~ bos6_1

bos1_2 ~~ bos1_2

bos2_2 ~~ bos2_2

bos3_2 ~~ bos3_2

bos4_2 ~~ bos4_2

bos5_2 ~~ bos5_2

bos6_2 ~~ bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML")

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

weak_boredom <- ' 

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# model estimation

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML")

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# model estimation

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML")

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

# 3. Model Comparison ----------------------------------------------------

# * 3.1 Steinmetz, 2014 ---------------------------------------------------

# corrected chi-square differences test

CorrDiff <- function(fit.baseline, fit.restricted) # basline = model with less degrees of freedom

{

  Corrc2Base <- fitMeasures(fit.baseline)[[6]] # corrected chi-square value

  MLc2Base <- fitMeasures(fit.baseline)[[3]] # uncorrected chi-square value

  dfBase <- fitMeasures(fit.baseline)[[4]] # degrees of freedom

  

  Corrc2Rest <- fitMeasures(fit.restricted)[[6]]

  MLc2Rest <- fitMeasures(fit.restricted)[[3]]

  dfRest <- fitMeasures(fit.restricted)[[4]]

  

  c0 <- MLc2Rest / Corrc2Rest

  c1 <- MLc2Base / Corrc2Base

  cd <- (dfRest*c0-dfBase*c1)/(dfRest-dfBase)

  CorrDiff <- (MLc2Rest-MLc2Base)/cd # corrected difference test statistic T_d^*

  dfDiff <- dfRest - dfBase

  pDiff <- 1-pchisq(CorrDiff, dfDiff)

  result <- list(Corr.chisq.difference = CorrDiff,

                 df.difference=dfDiff, p.value=pDiff)

  return(result)

}

# model comparison

CorrDiff(fit_configural_boredom, fit_weak_boredom)

CorrDiff(fit_weak_boredom, fit_strong_boredom)

CorrDiff(fit_strong_boredom, fit_strict_boredom)

# 3.2 PennState Alternative -----------------------------------------------

# according to: https://quantdev.ssri.psu.edu/tutorials/intro-basics-longitudinal-measurement-invariance

# compare model fit statistics

round(cbind(configural=inspect(fit_configural_boredom, 'fit.measures'), 

            weak=inspect(fit_weak_boredom, 'fit.measures'),

            strong=inspect(fit_strong_boredom, 'fit.measures'),

            strict=inspect(fit_strict_boredom, 'fit.measures')),3)

# chi-square difference test for nested models 

anova(fit_configural_boredom, fit_weak_boredom)

anova(fit_weak_boredom, fit_strong_boredom)

anova(fit_strong_boredom, fit_strict_boredom)

# 4. Fitting Boredom Growth in Lavaan using Strong Invariance -------------

# model specification (first order factor means fixed to 0)

boredom_1 ~ 0*1

boredom_2 ~ 0*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ 0*boredom_1

boredom_0 ~~ 0*boredom_2

boredom_1 ~~ 0*boredom_2

# adding the growth model

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

'

# model estimation

fit_strong_boredom_growth <- growth(strong_boredom_growth, data=df1, missing = "FIML")

summary(fit_strong_boredom_growth, fit.measures = TRUE, standardized = TRUE)

# * 4.1 Adding Boredom Growth Parameters to Dataframe ---------------------

idx <- lavInspect(fit_strong_boredom_growth, "case.idx")

fscores <- lavPredict(fit_strong_boredom_growth, assemble = TRUE)

## loop over factors

for (fs in colnames(fscores)) {

  df1[idx, fs] <- fscores[, fs]

}

[Patrick (Malone Quantitative)]

I haven't been following this thread closely, but what's the benefit of the first right-hand-side term in

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 +

(and similar for _1 and _2). I can't figure out what it would do, versus

boredom_0 =~ lambda1*bos1_0 +

Also, I'd recommend using the sem() wrapper instead of the growth() wrapper, since growth() is based on, as far as I know, growth of observed variables, not latents. You specify so much else in your model, I think sem() would work fine and be less likely to run into unintended defaults. Though you may need to add:

i ~ 1

s ~ 1

Pat

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/8c2bc5a4-e34e-4526-a29b-0b662e94b644n%40googlegroups.com.

[Manuel]

Dear Pat,

As commented in the code (in the note directly above that line of code) this is done to force this factor loading to be free. I therefore pre-multiply it with NA, as a hint to Lavaan that the value of this parameter is still unknown.

I tried the growth and sem function and attached the two outputs as pdf. Using the sem function without i ~ 1 and s ~ 1 intercept and slope would be 0.

I find it a bit spooky that the results differ when using different functions. Are you certain that I can trust the results of the sem function?

Thanks,

Manuel

[Patrick (Malone Quantitative)]

Manuel,

Interleaved:

On Mon, Dec 14, 2020 at 5:44 AM Manuel <mailsc...@gmail.com> wrote:

Dear Pat,

As commented in the code (in the note directly above that line of code) this is done to force this factor loading to be free. I therefore pre-multiply it with NA, as a hint to Lavaan that the value of this parameter is still unknown.

I understand that, but you shouldn't need to do *both* that *and* lambda1*bos1_0 in the same measurement equation. The latter makes the former redundant.

I tried the growth and sem function and attached the two outputs as pdf. Using the sem function without i ~ 1 and s ~ 1 intercept and slope would be 0.

I find it a bit spooky that the results differ when using different functions. Are you certain that I can trust the results of the sem function?

That's why I said you'd need to specify those for the sem() wrapper. As far as which to trust, it looks like there's a slight shift in the observed-variable handling--the intercepts all differ by .018, which is reflected in a difference of .018 in the other direction in the intercept of i. The model appears to be otherwise identical, and the fit statistics are identical. The difference appears to be an arbitrary artifact that won't affect your interpretation.

Pat

 

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/8668c7e6-8035-4d1d-8d8e-13034d1d61cbn%40googlegroups.com.

[Manuel]

Dear Pat,

You are absolutely right! I am getting there…

I changed my code accordingly to establish longitudinal strong measurement invariance as follows:

# exploring measurement invariance over time

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ lambda0*bos1_1 + bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ lambda0*bos1_2 + bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

# same factor loadings for each item over time

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

Unfortunately, if I am comparing the models now using the code from Steinmetz (2014) the first p-value is NA:

# corrected chi-square differences test

  CorrDiff <- function(fit.baseline, fit.restricted) # basline = model with less degrees of freedom

  {

  Corrc2Base <- fitMeasures(fit.baseline)[[6]] # corrected chi-square value

  MLc2Base <- fitMeasures(fit.baseline)[[3]] # uncorrected chi-square value

  dfBase <- fitMeasures(fit.baseline)[[4]] # degrees of freedom

  

  Corrc2Rest <- fitMeasures(fit.restricted)[[6]]

  MLc2Rest <- fitMeasures(fit.restricted)[[3]]

  dfRest <- fitMeasures(fit.restricted)[[4]]

  

  c0 <- MLc2Rest / Corrc2Rest

  c1 <- MLc2Base / Corrc2Base

  cd <- (dfRest*c0-dfBase*c1)/(dfRest-dfBase)

  CorrDiff <- (MLc2Rest-MLc2Base)/cd # corrected difference test statistic T_d^*

  dfDiff <- dfRest - dfBase

  pDiff <- 1-pchisq(CorrDiff, dfDiff)

  result <- list(Corr.chisq.difference = CorrDiff,

                 df.difference=dfDiff, p.value=pDiff)

  return(result)

  }

# model comparison

CorrDiff(fit_configural_boredom, fit_weak_boredom)

CorrDiff(fit_weak_boredom, fit_strong_boredom)

Warning:

In pchisq(CorrDiff, dfDiff) : NaNs produced

Furthermore, if I am then using the sem or growth function to calculate my growth model using the longitudinal scalar invariance model with the addition by Mauricio I get negative lv variances:

strong_boredom_growth <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

i ~ 1

s ~ 1

'

# model estimation

fit_strong_boredom_growth <- sem(strong_boredom_growth, data=df1, missing = "FIML", estimator = "MLR")

Warning:

In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

Thanks,

Manuel

[Mauricio Garnier-Villarreal]

Manuel

I think you are mixing contrainst based on identification methods. You are setting the marker variable identification , and using constraints like fixed variance.

Would make these changes, changes are highlighted

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ bos1_0 + bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ bos1_1 + bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ bos1_2 + bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

# same factor loadings for each item over time

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

'

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2 

boredom_0 ~0*1

boredom_1 ~NA*1

boredom_2 ~NA*1

'

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

############

Why are you building that function for the LRT between models? Becuase you are using the MLR estimator?

lavaan methods of nested model comparison do the correct adjust for MLR, MLM, DWLS, etc, estimator. WIth the lavTestLRT() for nested model comparison. If you go into details, you will some options for adjustsments based on the estimator

#########

Last, when you have a bagative LV variance, for which latent variable is it? If it for the Intecept or slope might be an indicationg that this is not a random parameter but should be treated as fix, which can be odd on the SEM, but its standard in growth curves in MLM.

[Manuel]

Dear Mauricio,

Thanks for your help! I did not know that and your solution looks much more elegant. I was building my models based on different tutorials as for example this one from Pennsylvania State University.
I have a few questions to your approach.

configural_boredom <- '

# measurement model

boredom_0 =~ bos1_0 + bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ bos1_1 + bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ bos1_2 + bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

Do we have to use the first item bos1 from all three timepoints twice? If not it would for example be:

boredom_0 =~bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

 

I am wondering the same about the weak and strong invariance models.

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

'

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

 

strong_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2 

boredom_0 ~0*1

boredom_1 ~NA*1

boredom_2 ~NA*1

'

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

 

Could we strike lambda0*bos1_0 for all time points and just use as for example:

boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

Moreover, the configural and weak models are working fine now. but the strong model still results in a warning:

In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    The variance-covariance matrix of the estimated parameters (vcov)

    does not appear to be positive definite! The smallest eigenvalue

    (= 2.594362e-17) is close to zero. This may be a symptom that the

    model is not identified.

[Mauricio Garnier-Villarreal]

Manuel

I didnt notice you had item1 repated. No need to do that, with this method you include each item only once

When you said that the strong invariance model is giving that error:

- You still have all the intercept equality constraints, right?

- When you do these model comparisons, what would you conclude about the model constraints?

- You can use the function lavTestScore() to look at univariate LRT tests for each equality constraint. This can help identify if a specific item should not be constraint equal

- Can you share theoutput from the summary(strong_model, fit.measures=T, standardized=T)?

[Manuel]

Dear Mauricio,

I have indeed deleted the intercept constraints. I have just put them back into the model (I hope the right way now) and it looks as follows:

# * 2.1 Configural Invariance ---------------------------------------------

configural_boredom <- '

# measurement model

boredom_0 =~ bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

fit_configural_boredom <- cfa(configural_boredom, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_configural_boredom, fit.measures=T, standardized=T)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

# same factor loadings for each item over time

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

'

# model estimation

fit_weak_boredom <- cfa(weak_boredom, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_weak_boredom, fit.measures=T, standardized=T)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_1 ~ NA*1

boredom_2 ~ NA*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

'

# model estimation

fit_strong_boredom <- cfa(strong_boredom, data=df1, missing="FIML", estimator="MLR", std.lv=T)

lavTestScore(fit_strong_boredom)

summary(fit_strong_boredom, fit.measures=T, standardized=T, modindices=T)

Now all models run without error. But afterwards I am using this code to test measurement invariance:

# * 3.1 Steinmetz, 2014 ---------------------------------------------------

# corrected chi-square differences test

  CorrDiff <- function(fit.baseline, fit.restricted) # basline=model with less degrees of freedom

  {

  Corrc2Base <- fitMeasures(fit.baseline)[[6]] # corrected chi-square value

  MLc2Base <- fitMeasures(fit.baseline)[[3]] # uncorrected chi-square value

  dfBase <- fitMeasures(fit.baseline)[[4]] # degrees of freedom

  

  Corrc2Rest <- fitMeasures(fit.restricted)[[6]]

  MLc2Rest <- fitMeasures(fit.restricted)[[3]]

  dfRest <- fitMeasures(fit.restricted)[[4]]

  

  c0 <- MLc2Rest / Corrc2Rest

  c1 <- MLc2Base / Corrc2Base

  cd <- (dfRest*c0-dfBase*c1)/(dfRest-dfBase)

  CorrDiff <- (MLc2Rest-MLc2Base)/cd # corrected difference test statistic T_d^*

  dfDiff <- dfRest - dfBase

  pDiff <- 1-pchisq(CorrDiff, dfDiff)

  result <- list(Corr.chisq.difference=CorrDiff,

                 df.difference=dfDiff, p.value=pDiff)

  return(result)

  }

# model comparison

CorrDiff(fit_configural_boredom, fit_weak_boredom)

CorrDiff(fit_weak_boredom, fit_strong_boredom)

This results in the error:

In pchisq(CorrDiff, dfDiff) : NaNs produced

And I am not getting a p-value for the chi square difference test of fit_configural_boredom and fit_weak_boredom.

How are you usually testing measurement invariance?

[Patrick (Malone Quantitative)]

Manuel,

Mauricio already asked you why you're writing your own model comparison function instead of using the established and well-debugged lavTestLRT() ?

Pat

--

You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/192ebed6-1314-4448-ad81-7ed220246ea3n%40googlegroups.com.

[Mauricio Garnier-Villarreal]

Manuel

Here is the code with some modifications. And as Pat mentioned, I use the lavaan function for model comparison, it will use the appropriate adjustement for using MLR. You can also use the compareFit() function from semTools that present multiple indices comparison

# * 2.1 Configural Invariance ---------------------------------------------

configural_boredom <- '
# measurement model
boredom_0 =~ bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0
boredom_1 =~ bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1
boredom_2 =~ bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

# intercepts

bos1_0 ~ 1

bos2_0 ~ 1
bos3_0 ~ 1
bos4_0 ~ 1
bos5_0 ~ 1
bos6_0 ~ 1

bos1_1 ~ 1

bos2_1 ~ 1
bos3_1 ~ 1
bos4_1 ~ 1
bos5_1 ~ 1
bos6_1 ~ 1

bos1_2 ~ 1

bos1_1 ~~ bos1_2
bos2_0 ~~ bos2_1 + bos2_2
bos2_1 ~~ bos2_2
bos3_0 ~~ bos3_1 + bos3_2
bos3_1 ~~ bos3_2
bos4_0 ~~ bos4_1 + bos4_2
bos4_1 ~~ bos4_2
bos5_0 ~~ bos5_1 + bos5_2
bos5_1 ~~ bos5_2

bos6_0 ~~ bos6_1 + bos6_2

bos6_1 ~~ bos6_2

# latent variables means
boredom_0 ~ 0*1

boredom_1 ~ 0*1
boredom_2 ~ 0*1

# latent variables variances and covariances
boredom_0 ~~ 1*boredom_0

boredom_1 ~~ 1*boredom_1
boredom_2 ~~ 1*boredom_2

'

# model estimation
fit_configural_boredom <- cfa(configural_boredom, data=df1,
                              missing="FIML", estimator="MLR", std.lv=T)
summary(fit_configural_boredom, fit.measures=T, standardized=T)


# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)
# same factor loadings for each item over time
weak_boredom <- '
# measurement model
boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0
boredom_1 =~ lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1
boredom_2 =~ lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

# intercepts

bos1_0 ~ 1

bos2_0 ~ 1
bos3_0 ~ 1
bos4_0 ~ 1
bos5_0 ~ 1
bos6_0 ~ 1

bos1_1 ~ 1

bos2_1 ~ 1
bos3_1 ~ 1
bos4_1 ~ 1
bos5_1 ~ 1
bos6_1 ~ 1

bos1_2 ~ 1

bos1_1 ~~ bos1_2
bos2_0 ~~ bos2_1 + bos2_2
bos2_1 ~~ bos2_2
bos3_0 ~~ bos3_1 + bos3_2
bos3_1 ~~ bos3_2
bos4_0 ~~ bos4_1 + bos4_2
bos4_1 ~~ bos4_2
bos5_0 ~~ bos5_1 + bos5_2
bos5_1 ~~ bos5_2

bos6_0 ~~ bos6_1 + bos6_2

bos6_1 ~~ bos6_2

# latent variables means
boredom_0 ~ 0*1

boredom_1 ~ 0*1
boredom_2 ~ 0*1

bos1_1 ~~ bos1_2
bos2_0 ~~ bos2_1 + bos2_2
bos2_1 ~~ bos2_2
bos3_0 ~~ bos3_1 + bos3_2
bos3_1 ~~ bos3_2
bos4_0 ~~ bos4_1 + bos4_2
bos4_1 ~~ bos4_2
bos5_0 ~~ bos5_1 + bos5_2
bos5_1 ~~ bos5_2

bos6_0 ~~ bos6_1 + bos6_2

bos6_1 ~~ bos6_2

# latent variables means
boredom_0 ~ 0*1

boredom_1 ~ NA*1
boredom_2 ~ NA*1

# latent variables variances and covariances
boredom_0 ~~ 1*boredom_0
boredom_1 ~~ NA*boredom_1
boredom_2 ~~ NA*boredom_2
'

# model estimation
fit_strong_boredom <- cfa(strong_boredom, data=df1,
                          missing="FIML", estimator="MLR", std.lv=T)
lavTestScore(fit_strong_boredom)
summary(fit_strong_boredom, fit.measures=T, standardized=T, modindices=T)

lavTestLRT(fit_configural_boredom,
           fit_weak_boredom)

lavTestLRT(fit_weak_boredom,
           fit_strong_boredom)

semTools::compareFit(fit_configural_boredom,
                     fit_weak_boredom,
                     fit_strong_boredom)

[Manuel]

Dear Pat and Mauricio,

I must have overseen this. lavTestLRT() and semTools::compareFit work great! Thank you.

I have talked to my professor about this today and we came up with a code that is similar to yours Mauricio. In both codes the models on their own seem to work fine but there is measurement non-invariance. I highlighted the differences to your code below and attached the model comparison outputs from your and my code.

A reason for non-invariance could be that at time point 0 the questionnaire was filled out at home and not at school like at time point 1 and 2.

Can I still use the intercept and slope of the strong_boredom model to test if the trajectories of boredom linked with health?

bos1_0 ~~ bos1_1 + bos1_2

bos1_1 ~~ bos1_2

bos2_0 ~~ bos2_1 + bos2_2

bos2_1 ~~ bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos3_1 ~~ bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos4_1 ~~ bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos5_1 ~~ bos5_2

bos6_0 ~~ bos6_1 + bos6_2

bos6_1 ~~ bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_strong_boredom <- cfa(strong_boredom, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom, fit.measures=T, standardized=T)

[Mauricio Garnier-Villarreal]

Manuel

I took these constriant out, because you dont need to constraint interceots until you are in the string invariance model. I dont see why these constraints are there for the configural and weak models?

bos1_0 ~ i1*1

For the same these lines. In the condigural and weak models, these means are basically the item means for which you have the intercepts constrained. These means are not taking into account all the items until the strong invariance model. Only when you have multiple intercepts equated you can say that the factor mean is different than a single indicator mean

boredom_1 ~ 1

boredom_2 ~ 1

On the same topic these lines. Until you have equated factor loadings across time, these factor variances are not meaningful. This because I used factor variance method of identification (std.lv=T). So, these are not properly specified for the configural model

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

#####

About the model comaprison (I ony looked at the output from my code). You could argue that the differences have a small magnitude, based on the CFI and SRMR model differences. There is a disagreement in the literature if these comaprison shdould be done with LRT or CFI/RMSEA...

Cheung, G. W., & Rensvold, R. B. (2002). Evaluating goodness-of-fit indexes for testing measurement invariance. Structural Equation Modeling, 9(2), 233–255.

Jorgensen, T. D., Kite, B. A., Chen, P.-Y., & Short, S. D. (2018). Permutation randomization methods for testing measurement equivalence and detecting differential item functioning in multiple-group confirmatory factor analysis. Psychological Methods, 23(4), 708–728. https://doi.org/10.1037/met0000152

Liang, X., & Luo, Y. (2020). A Comprehensive Comparison of Model Selection Methods for Testing Factorial Invariance. Structural Equation Modeling: A Multidisciplinary Journal, 27(3), 380–395. https://doi.org/10.1080/10705511.2019.1649983

Meade, A. W., Johnson, E. C., & Braddy, P. W. (2008). Power and sensitivity of alternative fit indices in tests of measurement invariance. Journal of Applied Psychology, 93(3), 568–592. https://doi.org/10.1037/0021-9010.93.3.568

You should use the function lavTestScore(weakmodel), these will do a LRT for each specific equality constraint in the model. So, you can identify if a single item may be problematic, if these is the case you can relase the equality constraint for only that 1 item. And if you are satisfied with that model, you can continue, this is called partial-measurement invariance. Its fine to continue if you have partial-invariance. The recommensation I make is to seek to have at least half the items being invariant to be able to make reliable inferences for latent parameters (like your growth curve). This comes from personla experience, I havent seen research suggesting a minimum threshold of items to be invariant

[Manuel]

Dear Mauricio,

Thank you for the explanation. It makes perfect sense to me now.

Since I have a large sample and the differences in CFI (-0.002) and RMSEA (-0.001) are really small I would like to refrain from releasing constraints for single items.

In the next step, I added the growth model to the strict_boredom as follows:

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

i ~ 1

s ~ 1

'

fit_strong_boredom_growth <- cfa(strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth, fit.measures=T, standardized=T)

Unfortunately, this results in the error:

1: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    The variance-covariance matrix of the estimated parameters (vcov)

    does not appear to be positive definite! The smallest eigenvalue

    (= 4.907922e-18) is close to zero. This may be a symptom that the

    model is not identified.

2: In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

After solving this, I would add the boredom growth parameters to my data frame:

idx <- lavInspect(fit_strong_boredom_growth, "case.idx")

fscores <- lavPredict(fit_strong_boredom_growth, assemble=T)

## loop over factors

for (fs in colnames(fscores)) {

  df1[idx, fs] <- fscores[, fs]

}

And then investigate the growth parameters of boredom are linked with health:

boredom_health <- '

scale1 + scale2 + scale3 + scale4 ~ i

scale1 + scale2 + scale3 + scale4 ~ s

'

fit_boredom_health <- sem(model=boredom_health, data=df1, missing ="FIML", estimator="MLR")

summary(fit_boredom_health, fit.measures=T, standardized=T)

[Mauricio Garnier-Villarreal]

Manuel

I think that is a defensible argument, would still recomend you to look at lavtestscore() to make sure there is no item that is misbehaving.

When you add the latent growth curve, you also need to fix all the first order factor means to 0, and estimate the variances and means of the growth curves. This because you want to estimate mean bhavior on the growth curve.

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

i ~NA*1

s ~ NA*1

i ~~ NA*i

  s~~ NA*s

boredom_0 ~ 0*1

boredom_1 ~0*1

boredom_2 ~0*1

###########

About combining the factor scores for the growth curve factors into your data set and use that in posterior analysis. This is NOT recommended. Regressions with factor scores will be bias, unless you do some method to account for correction, like Croon's bias correction, or plausible values. Croons correction is still an experimental feature on lavaan, and you could use plausible values with the semTools functions. But any of this options require extra steps. Just adding your estimated factor scores into your data set will give you bias regressions.

Or do the latent regression, you could just add this regressions into the larger model, and keep adding to the model syntax. Not sure why you decide to go with factor scores automatically. Just add those regressions to same growth curve model syntax

Devlieger, I., Mayer, A., & Rosseel, Y. (2016). Hypothesis Testing Using Factor Score Regression: A Comparison of Four Methods. Educational and Psychological Measurement, 76(5), 741–770. https://doi.org/10.1177/0013164415607618

Devlieger, I., Talloen, W., & Rosseel, Y. (2019). New Developments in Factor Score Regression: Fit Indices and a Model Comparison Test. Educational and Psychological Measurement, 79(6), 1017–1037. https://doi.org/10.1177/0013164419844552

Asparouhov, T. & Muthen, B. O. (2010). Plausible values for latent variables using Mplus. Technical Report. Retrieved from www.statmodel.com/download/Plausible.pdf

[Manuel]

Dear Mauricio,

I looked at lavTestScore for fit_weak_boredom and fit_strong_boredom and no epc was even close to being greater than 0.1

Is it recommended to report this in a paper as well?

I added your suggestions to the code that now looks as follows and works flawlessly:

strong_boredom_growth <- '

# measurement model

boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ i2*1

bos3_0 ~ i3*1

bos4_0 ~ i4*1

bos5_0 ~ i5*1

bos6_0 ~ i6*1

bos1_1 ~ i1*1

bos2_1 ~ i2*1

bos3_1 ~ i3*1

bos4_1 ~ i4*1

bos5_1 ~ i5*1

bos6_1 ~ i6*1

bos1_2 ~ i1*1

bos2_2 ~ i2*1

bos3_2 ~ i3*1

bos4_2 ~ i4*1

bos5_2 ~ i5*1

bos6_2 ~ i6*1

# unique variances and covariances

# latent variables means (first order factor means fixed to 0)

boredom_0 ~ 0*1

boredom_1 ~ 0*1

boredom_2 ~ 0*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

# adding the growth model

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

# estimate variances and means of the growth curve

i ~ NA*1

s ~ NA*1 

i ~~ NA*i

s ~~ NA*s

'

# model estimation

fit_strong_boredom_growth <- cfa(strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth, fit.measures=T, standardized=T)

I was not aware that the factor scores would be biased. Adding the regressions to same growth curve model resulted in an error:

strong_boredom_growth_health <- '

# measurement model

boredom_0 =~ lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ i2*1

bos3_0 ~ i3*1

bos4_0 ~ i4*1

bos5_0 ~ i5*1

bos6_0 ~ i6*1

bos1_1 ~ i1*1

bos2_1 ~ i2*1

bos3_1 ~ i3*1

bos4_1 ~ i4*1

bos5_1 ~ i5*1

bos6_1 ~ i6*1

bos1_2 ~ i1*1

bos2_2 ~ i2*1

bos3_2 ~ i3*1

bos4_2 ~ i4*1

bos5_2 ~ i5*1

bos6_2 ~ i6*1

# unique variances and covariances

# latent variables means (first order factor means fixed to 0)

boredom_0 ~ 0*1

boredom_1 ~ 0*1

boredom_2 ~ 0*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

# adding the growth model

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

# estimate variances and means of the growth curve

i ~ NA*1

s ~ NA*1 

i ~~ NA*i

s ~~ NA*s

# multiple regression with health dimensions

Dimension1 + Dimension2 + Dimension3 + Dimension4 + Dimension5 ~ i

Dimension1 + Dimension2 + Dimension3 + Dimension4 + Dimension5 ~ s

'

fit_strong_boredom_growth_health <- sem(strong_boredom_growth_health, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth_health, fit.measures=T, standardized=T)

In lav_model_estimate(lavmodel = lavmodel, lavpartable = lavpartable,  :

  lavaan WARNING: the optimizer warns that a solution has NOT been found!

[Mauricio Garnier-Villarreal]

Manuel

That depends on how you want to present it, but its talking about the univariate tests can help the argument that weak and strong invariance hold

You are etimaing a large model now. The estimator might need longer to find a solution. You can add the argument verbose=T to see how far is the the estimation method. Or later you can look into th control options of the optimizier to have it run for longer to reach convergence. I dont remeber the arguments for controling the optimizer

[Manuel]

Hi Mauricio,

based on Yves slides I added bounds=T and the estimation ends normally after 400 iterations. Unfortunately, this way non of the regressions get significant as I would expect from previous correlation analysis or when using the growth curve factors in posterior analysis (I know its bad practice).

Allowing 20000 iterations to reach a solution as described here results in the error:

lavdata@nlevels > 1L is not TRUE

[Mauricio Garnier-Villarreal]

Manuel

The second error happens because those are arguments for the Multilevel Modeling. You its looking for multilevel data

When you have a convergent result with bounds=T. Is there any reason you think that solution is worng or have some estimation problem? Do you have a large sample? As this is now a farily large model.

It just might be that those regression are "not signifficant" (I hate this term). And that is the best generalization you can do. You can still describe the effect sizes based on this.

If you post the model output i can see if there is any notorious issue

[Manuel]

Dear Mauricio,

Here is the model output including a correlation matrix as well as the model output of an inverted regression.

Thank you so much for your help!

[Mauricio Garnier-Villarreal]

Manuel

I dont see any particular problem in the model parameters. Both match the "observed" correlation matrix, as the correlations are small and negative. The main reason you get "signifficant" correlations is your large sample.

Regressions in either direction match this, by showing small effect sizes in general. Now, the regressions will also be different by estimating the conditional relation (controlling for the other predcitors effect), while the correlations are looking at the bivariate relation

Now, the direction of the regressions are fundamentally different research questions as you are stating the directionality. I obviously dont know the theoreticall background here, but be very concious about having a strong theoretical defense for which ever directio you use.

If you are not interested in the direction, and only looking for the "relation" between them you can estimate the correlations between them in lavaan.

But this depends on the question and theory you are looking to test

[Manuel]

Dear Mauricio,

Thanks again for your help and a happy new year 🙂

After going through literature in the past weeks I decided to look only at the relation between the variables. How can this be done in Lavaan?

I could not find any information so I tried adding all variables using i ~~ KC27ph_T etc. but this results in the error:

lav_model(lavpartable = lavpartable, lavoptions = lavoptions,  : 

  lavaan ERROR: parameter is not defined: i ~~ KC27ph_T

I am also interested in partial correlations.

Thanks,

Manuel

[Mauricio Garnier-Villarreal]

Manuel

Can you share the code that creates this error. and the output of str() for your data?

For partial correlations you cn estimate them by adding the covariate of interest predicting both variables of the correlation

[Manuel]

Hi Mauricio,

I added the correlations to my strong boredom model as follows:

# model specification (first order factor means fixed to 0)

health_strong_boredom_growth <- '

# correlation with health dimensions

i ~~ KC27ph_T

i ~~ KC27pw_T

i ~~ KC27pa_T

i ~~ KC27pe_T

i ~~ KC27sc_T

s ~~ KC27ph_T

s ~~ KC27pw_T

s ~~ KC27pa_T

s ~~ KC27pe_T

s ~~ KC27sc_T

'

# model estimation

fit_health_strong_boredom_growth <- lavaan::cfa(health_strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_health_strong_boredom_growth, fit.measures=T, standardized=T)

Here is the str of my data frame:

> str(df1)

'data.frame': 1229 obs. of  69 variables:

 $ codfbs1_2: Factor w/ 1229 levels "104-2868","104-2892",..: 1 2 3 4 5 6 8 9 10 11 ...

 $ bos1_0   : int  2 2 1 2 1 3 2 2 2 1 ...

 $ bos2_0   : int  3 1 1 1 3 2 2 1 1 1 ...

 $ bos3_0   : int  1 1 2 1 3 2 1 2 1 1 ...

 $ bos4_0   : int  1 2 4 2 4 5 4 2 4 1 ...

 $ bos5_0   : int  3 2 3 2 1 4 4 2 3 1 ...

 $ bos6_0   : int  4 2 2 1 1 4 4 2 2 1 ...

 $ bos1_1   : int  2 2 2 3 2 2 4 4 2 1 ...

 $ bos2_1   : int  2 2 1 2 4 1 2 1 1 1 ...

 $ bos3_1   : int  2 1 1 1 3 1 2 2 1 1 ...

 $ bos4_1   : int  3 1 2 2 4 5 5 4 3 1 ...

 $ bos5_1   : int  2 2 2 2 2 3 5 3 4 1 ...

 $ bos6_1   : int  2 2 1 2 1 3 1 2 3 1 ...

 $ bos1_2   : int  2 3 2 3 2 1 3 3 2 1 ...

 $ bos2_2   : int  2 1 3 1 4 3 2 3 1 1 ...

 $ bos3_2   : int  2 2 2 1 4 3 2 3 1 1 ...

 $ bos4_2   : int  2 3 5 2 5 5 5 4 3 1 ...

 $ bos5_2   : int  2 2 3 2 4 4 4 3 2 2 ...

 $ bos6_2   : int  3 2 1 1 2 5 2 3 1 1 ...

 $ KC27ph_T : num  47.1 44.7 40.5 44.7 59.4 ...

 $ KC27pw_T : num  36.7 48.5 46.5 48.5 53.1 ...

 $ KC27pa_T : num  53.3 74.4 46.5 55.8 59.1 ...

 $ KC27pe_T : num  66.3 66.3 66.3 49.8 66.3 ...

 $ KC27sc_T : num  48.1 48.1 51.1 48.1 51.1 ...

Could you provide an example of how I would estimate partial correlations by "adding the covariate of interest predicting both variables of the correlation"?

Thanks

Manuel

[Mauricio Garnier-Villarreal]

Manuel

Actually I dont know why you are getting that error. I dont see a clear mistake on the syntax. Was thinking might have been a misspelled of a variable name. Are you sure you have the latest lavaan version installed?

Could try changing the name of the health variables, or standardizing them with the function scale(). But these are just guesses

About a partial correlation. For example I have variables A, B, and C. And I want to see the partial correlation between A and B, partialing out C. In the case C would be the covariate. Would look like this in lavaan syntax

'

A ~~ B

A ~ C

 B ~ C

'

So the artial correlation would be between the residuals of A and B after regressing C on them

[Manuel]

Hi,

as always, you are right and helped a lot 🙂

Updating Lavaan solved the issue.

In a first step I added the correlations to my strong boredom model and compared the results to a manifest correlation matrix.

health_boredom_correlation <- '

# correlation with health dimensions (Rasch t-values)

boredom_0 ~~ KC27ph_T

boredom_0 ~~ KC27pw_T

boredom_0 ~~ KC27pa_T

boredom_0 ~~ KC27pe_T

boredom_0 ~~ KC27sc_T

boredom_1 ~~ KC27ph_T

boredom_1 ~~ KC27pw_T

boredom_1 ~~ KC27pa_T

boredom_1 ~~ KC27pe_T

boredom_1 ~~ KC27sc_T

boredom_2 ~~ KC27ph_T

boredom_2 ~~ KC27pw_T

boredom_2 ~~ KC27pa_T

boredom_2 ~~ KC27pe_T

boredom_2 ~~ KC27sc_T

'

# model estimation

fit_health_boredom_correlation <- lavaan::cfa(health_boredom_correlation, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_health_boredom_correlation, fit.measures=T, standardized=T)

When looking at the estimates they sometimes have different operators (+/-) than in the manifest correlation matrix (see attached pdf). Also different constructs get significant. Is this because of the latent construct modeling or am I doing something wrong / misinterpreting the estimates?

In a second step I added the correlations to the boredom growth model:

Looking at the output (see attached pdf) the intercept now only correlates significantly with KC27sc_T.

Thanks

Manuel

[Mauricio Garnier-Villarreal]

Manuel

##

When looking at the estimates they sometimes have different operators (+/-) than in the manifest correlation matrix (see attached pdf). Also different constructs get significant. Is this because of the latent construct modeling or am I doing something wrong / misinterpreting the estimates?

##

I dont see any syntax issues here. So, I would tend to think that the differences are due to looking at relations between observed and latent variables. Also, after you have strong invariance

##

In a second step I added the correlations to the boredom growth model:

Looking at the output (see attached pdf) the intercept now only correlates significantly with KC27sc_T.

##

Now, these relations have no similar with manifest variables, as you dont have a measure intercept/slope. But relating them to characteristics of change. So, wouldnt be sorprise if the relations are "different"

Glad this has help

[Manuel]

Dear Mauricio,

I don't know how to explain this properly, but to me it seems like adding all the variables using ~~ at once is not the same as a "simple" correlation matrix. Maybe it is testing if not one but all variables correlate and therefore the results differ that much from the results that I get using manifest variables?

If I only add one of the variables to the strong boredom model as follows, the results match the ones I get using manifest variables:

strong_boredom_pw <- '

boredom_0 ~ 0*1

boredom_1 ~ NA*1

boredom_2 ~ NA*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

# correlation with KC27pw_T

boredom_0 ~~ KC27pw_T

boredom_1 ~~ KC27pw_T

boredom_2 ~~ KC27pw_T

'

# model estimation

fit_strong_boredom_pw <- lavaan::cfa(strong_boredom_pw, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_ph, fit.measures=T, standardized=T)

After running this model I would copied the code and only change the # correlation with KC27pw_T part to the next variable. This way my latent results match the manifest ones.

Strangely the covariances of the boredom time points are varying aroind +/- 0.001 between when I am running the models with only one variable one after another.

Thanks,

Manuel

[Mauricio Garnier-Villarreal]

Manuel

If you are just adding them by ~~ the correlations (std.all) should not change when adding more variables.

The only thing that comes to mind is that it is handling the covariates as fix

So, try with all the variables, and add the argument fixed.x=F to the cfa() function. And in the syntax add all the covariances between the covariates, like

KC27pw_T ~~ KC27pw_F +  KC27pw_l 

KC27pw_F ~~  KC27pw_l

This way you will be sure it is estimating the full covariance matrix between all the variables

[Manuel]

Great. You mean like this right?

strong_boredom_correlation <- '

# correlation with health dimensions

boredom_0 ~~ KC27ph_T

boredom_0 ~~ KC27pw_T

boredom_0 ~~ KC27pa_T

boredom_0 ~~ KC27pe_T

boredom_0 ~~ KC27sc_T

boredom_0 ~~ KC10in_T

boredom_1 ~~ KC27ph_T

boredom_1 ~~ KC27pw_T

boredom_1 ~~ KC27pa_T

boredom_1 ~~ KC27pe_T

boredom_1 ~~ KC27sc_T

boredom_1 ~~ KC10in_T

boredom_2 ~~ KC27ph_T

boredom_2 ~~ KC27pw_T

boredom_2 ~~ KC27pa_T

boredom_2 ~~ KC27pe_T

boredom_2 ~~ KC27sc_T

boredom_2 ~~ KC10in_T

# covariances between the covariates

KC27ph_T ~~ KC27pw_T + KC27pa_T + KC27pe_T + KC27sc_T + KC10in_T

KC27pw_T ~~ KC27pa_T + KC27pe_T + KC27sc_T + KC10in_T

KC27pa_T ~~ KC27pe_T + KC27sc_T + KC10in_T

KC27pe_T ~~ KC27sc_T + KC10in_T

KC27sc_T ~~ KC10in_T

'

# model estimation

fit_strong_boredom_correlation <- lavaan::cfa(strong_boredom_correlation, data=df1, missing="FIML", estimator="MLR", std.lv=T, fixed.x=F)

summary(fit_strong_boredom_correlation, fit.measures=T, standardized=T)

[Mauricio Garnier-Villarreal]

Yes, that is equivalent to estimating all the covariances without conditional relations

[Manuel]

Perfect. Thank you.

Is there a way to add Bonferroni correction to the p-values? Since all of them are 0.000 and options(digits=10) does not work I don't know how to account for multiple testing.

[Terrence Jorgensen]

Is there a way to add Bonferroni correction to the p-values?

You could just use an adjusted alpha level.  Or to adjust p values, you can use the p.adjust() function.  To access the vector of p values, use parameterEstimates()

PE <- parameterEstimates(fit_strong_boredom_correlation)

PE$p.bonf <- p.adjust(PE$pvalue, method = "bonferroni")

Since all of them are 0.000 and options(digits=10) does not work 

Use the nd= argument to lavaan's summary() method

summary(fit_strong_boredom_correlation, fit.measures=T, standardized=T, nd = 6)

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

[Manuel]

Dear all,

I hope you are doing well. I have another question :)

To recapitalize:

1. I established longitudinal strong measurement invariance for three time points of boredom
   
2. I added growth parameters (intercept and slope) to the strong boredom model to investigate the trajectories of boredom
3. I added health variables to this model that were raised at time point 3 to explore if changes in boredom across the three time points were correlatively linked with those variables

The results of 3. show that stronger increases of boredom (slope) are linked with two health variables, which I want to report.

But I am wondering – since at 2. the variance of this slope is not significantly different from zero (tau2 = .864, p = .055) – does it even make sense to explain the variance of the slope in 3.? On the other hand, the constraints of the strong model are limiting the variance.

I attached the results and will add my R script in the next post. I am looking forward to hearing your thoughts.

Thanks,

Manuel

[Manuel]

# model estimation

fit_strong_boredom <- lavaan::cfa(strong_boredom, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom, fit.measures=T, standardized=T)

# further insights

lavaan::modindices(fit_strong_boredom, sort.=T)

lavaan::lavTestScore(fit_strong_boredom, epc=T)

# * 2.4 model comparison --------------------------------------------------

library(semTools)

lavTestLRT(fit_configural_boredom, fit_weak_boredom)

lavTestLRT(fit_weak_boredom, fit_strong_boredom)

modelcomparison <- semTools::compareFit(fit_configural_boredom, fit_weak_boredom, fit_strong_boredom)

modelcomparison

# 3. Boredom ~~ (is correlated with) Health -------------------------------

# model specification (additional constraints on manifest intercepts)

summary(fit_strong_boredom_correlation, fit.measures=T, standardized=T, nd = 3)

# holm correction

PE <- parameterEstimates(fit_strong_boredom_correlation) # extract parameter estimates

PE$p.holm <- p.adjust(PE$pvalue, method = "holm") # adjust p values

options(scipen=999) # prevent scientific notation of p-values

PE

# 4. Fitting Boredom Growth in Lavaan using Strong Invariance -------------

# model specification (first order factor means fixed to 0)

fit_strong_boredom_growth <- lavaan::cfa(strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth, fit.measures=T, standardized=T, ci=T)

# 5. Boredom Growth ~~ (is correlated with) Health ------------------------

i ~~ KC10in_T

s ~~ KC27ph_T

s ~~ KC27pw_T

s ~~ KC27pa_T

s ~~ KC27pe_T

s ~~ KC27sc_T

s ~~ KC10in_T

# covariances between the covariates

KC27ph_T ~~ KC27pw_T + KC27pa_T + KC27pe_T + KC27sc_T + KC10in_T

KC27pw_T ~~ KC27pa_T + KC27pe_T + KC27sc_T + KC10in_T

KC27pa_T ~~ KC27pe_T + KC27sc_T + KC10in_T

KC27pe_T ~~ KC27sc_T + KC10in_T

KC27sc_T ~~ KC10in_T

'

# model estimation

fit_health_strong_boredom_growth <- lavaan::cfa(health_strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_health_strong_boredom_growth, fit.measures=T, standardized=T)

# holm correction

pe_growth <- parameterEstimates(fit_health_strong_boredom_growth) # extract parameter estimates

pe_growth$p.holm <- p.adjust(pe_growth$pvalue, method = "holm") # adjust p values

options(scipen=999) # prevent scientific notation of p-values

pe_growth

[Patrick (Malone Quantitative)]

Manuel,

The variance being not significant does not mean exactly zero -- especially with p = .055. Some users (or reviewers) think you shouldn't predict such a slope, but it's not uncommon in my experience to find what you found.

Pat

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/da70b0b3-991a-44ab-93f1-4da1c29f1e26n%40googlegroups.com.

--

Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

[Mauricio Garnier-Villarreal]

manuel

I agree with Patrick. Also, would add that that p-value assumes that the parameter is normally distributed, which can be hard to defend for variances (specially small ones)

[Manuel]

Thank you.

Do you have any arguments against users (or reviewers) why I should predict this slope anyways? And is it common to report the variance including its p-value if I want to report the results of 3.?

@Patrick

could you elaborate what you mean by "it's not uncommon in my experience to find what you found"?

[Terrence Jorgensen]

Do you have any arguments against users (or reviewers) why I should predict this slope anyways?

The fixed effect of an exogenous predictor on the latent slope might be (significantly) nonzero, even if the variance of random effects around that latent slope is not (significantly) different from zero. 

In a different but related situation:

I have fitted a model whose chi-squared is not significant, yet when I freed a parameter and conducted the chi-squared different test, the less restricted model fit significantly better.

 

And is it common to report the variance including its p-value if I want to report the results of 3.?

Not sure about what is common, but you can do so.  THE ML estimator in lavaan does not (by default) restrict variances to be non-negative, so their sampling distributions are in fact asymptotically normal and can include negative values (particularly in smaller samples, especially when the true parameter is not substantially larger than zero relative to its sampling error).  ML estimates of variances will have non-normal sampling distributions only when using constrained estimation to guarantee non-negative estimates (e.g., the default in multilevel modeling software).

[Michal Levinsky]

Dear Manuel,

Can you please share the code of the graph, or specify how you created it? I will be grateful.

[Manuel]

Dear Michael,

Sure. For this plot I just used manifest means and ggplot:

library(dplyr)

by_time <- group_by(df3, TimePoint)

s <- summarise(by_time, Boredom = mean(Boredom, na.rm=T))

myplot1 <- ggplot(df3, aes(x= TimePoint, y = Boredom)) +

  geom_point(size = 1, alpha=0.1, position=position_jitter(width=.3,height=.03)) +

  geom_line(alpha = .01, aes(group=codfbs1)) +

  stat_summary(alpha = .5, aes(group = class), fun="mean", geom="line", color = "#63BEEE", size = 1) + # add boredom class means to plot

  geom_line(data=s, aes(x = TimePoint, y = Boredom , group = 1), alpha = .9, color="#0F3FA0", size = 5) + # add boredom mean to plot

  theme_classic()

Does anybody have a code to plot latent growth curves including confidence intervals?

Thanks

[Michal Levinsky]

Thank you so much!!! 

בתאריך יום ג׳, 2 במרץ 2021, 10:35, מאת Manuel ‏<mailsc...@gmail.com>:

--
You received this message because you are subscribed to a topic in the Google Groups "lavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lavaan/dNLXNaRFhzA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lavaan+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/34e0f0bf-8c86-4349-9c65-63f39a2df260n%40googlegroups.com.

[Manuel]

Dear All,

I read a paper doing similar analyses and they added classes as clusters. Since my data is also from different classrooms I wanted to add cluster="class" which results in a warning. Would you recommend using clusters on all models?

fit_strong_boredom_correlation <- lavaan::cfa(strong_boredom_correlation, data=df1, missing="FIML", estimator="MLR", std.lv=T, cluster="class")

Warning:

In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    The variance-covariance matrix of the estimated parameters (vcov)

    does not appear to be positive definite! The smallest eigenvalue

    (= 2.464946e-14) is close to zero. This may be a symptom that the

    model is not identified.

In a next step I wanted to added gender as a grouping variable using group = "gender", group.equal = c("loadings", "intercepts"). Adding this to the model above solved the warning. But if I add this to my growth correlation model it resulted in a warning:

fit_health_strong_boredom_growth_gender <- lavaan::cfa(health_strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T, group = "s01_y01", group.equal = c("loadings", "intercepts"), cluster="class")

Warning:

In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

If I decide to add gender as a grouping variable, would I add this to the measurement invariance analysis as well?

Also, do you have a code to plot latent growth curves including confidence intervals?

Kind regards,

Manuel

[Manuel]

Dear all,

I am slowly but surely getting forward. Can you please help me with two more things?

First, am I right to believe that the means and intercepts of latent variables are arbitrary (in my case with the variance of the intercept and slope fixed to 1) and not directly interpretable (e.g. the intercept is not the mean of time point one)?

Second, I have been asked to compare my linear growth model to intercept-only and curvilinear models to provide a rationale for this approach.

I started by running the linear growth model, added a quadratic model and an intercept-only model, and compared the fits using lavTestLRT.

Is this the "correct" procedure?

Unfortunately, the quadratic model (2.) results in a warning:

In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

And this results in a warning when comparing the fits:

In lavTestLRT(fit_strong_boredom_growth, fit_strong_boredom_growth_q,  :

  lavaan WARNING: some models have the same degrees of freedom

Is there an error in my quadratic model?

Thank you very much and kind regards,

Manuel 

Here is my code:

# 1. Fitting Boredom Growth in Lavaan using Strong Invariance ------------------

# model specification (first order factor means fixed to 0)

# model estimation

fit_strong_boredom_growth <- lavaan::cfa(strong_boredom_growth, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth, fit.measures=T, standardized=T, ci=T)

# 2. Fitting Boredom Quadratic Growth in Lavaan using Strong Invariance --------

# model specification (first order factor means fixed to 0)

strong_boredom_growth_q <- '

q =~ 0*boredom_0 + 1*boredom_1 + 4*boredom_2 # q slope

# estimate variances and means of the growth curve

i ~ NA*1

q ~ NA*1 

i ~~ NA*i

q ~~ NA*q

'

# model estimation

fit_strong_boredom_growth_q <- lavaan::cfa(strong_boredom_growth_q, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_growth_q, fit.measures=T, standardized=T, ci=T)

# 3. Fitting Boredom Intercept in Lavaan using Strong Invariance ---------------

# model specification (first order factor means fixed to 0)

strong_boredom_intercept <- '

# estimate variances and means of the growth curve

i ~ NA*1

i ~~ NA*i

'

# model estimation

fit_strong_boredom_intercept <- lavaan::cfa(strong_boredom_intercept, data=df1, missing="FIML", estimator="MLR", std.lv=T)

summary(fit_strong_boredom_intercept, fit.measures=T, standardized=T, ci=T)

# * 4. Model Comparison --------------------------------------------------

lavTestLRT(fit_strong_boredom_growth, fit_strong_boredom_growth_q, fit_strong_boredom_intercept)

[Pat Malone]

Manuel,

The first warning is likely not an issue, if your results look generally plausible (specifically, none of the latent variable variances/disturbances is negative and large relative to its standard error).

The second suggests something wrong. The most likely issue is that, broadly speaking, a quadratic LGM needs at least four occasions.

Pat

You received this message because you are subscribed to the Google Groups "lavaan" group.

To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/eeff7dbf-1647-41e1-951b-16a10be002a0n%40googlegroups.com.

--

Patrick S. Malone, PhD

Sr Research Statistician, FAR HARBΦR

+1 803.553.4181  |  pat@  |  farharbor.com

This message may contain confidential information; if you are not the intended recipient please notify the sender and delete the message.

[Pat Malone]

Ah, I looked more closely.

The quadratic model should also include the linear slope. If you do that, your model will be under-identified because of too few occasions.

[Mauricio Garnier-Villarreal]

###

First, am I right to believe that the means and intercepts of latent variables are arbitrary (in my case with the variance of the intercept and slope fixed to 1) and not directly interpretable (e.g. the intercept is not the mean of time point one)?

####

the factor means are arbitrary based on the identification method you used. With fixed variance/mean method the means witll be 0 at time 1 and (after strong invariance) for the next times points will be in reference to the mean change compared to time 1. The intercept of growth curve is the average score at which time point you set the 0 in the time scale. This is also arbitrary as you can choose where to set your time 0. But this growth curve is interpretable, most commonly as the mean at the first time point, or start point of the longitudinal change

###

Second, I have been asked to compare my linear growth model to intercept-only and curvilinear models to provide a rationale for this approach.

###

The intercept only model will be done by excluding any other factor slopes, just one growth factor.

The linear model, would have the intercept factor, and the linear slope factors

The quadratic model would have, the intercept, linear and quadratic factors

Would recommend this book for these other extensions of models

Grimm, K. J., Ram, N., & Estabrook, R. (2016). Growth Modeling: Structural Equation and Multilevel Modeling Approaches.

A quadratic factor is consider a higher order term. This means that for estimating any higher order terms, you need to include all previous lower order terms, in this case, linear and intercept factors

[Manuel]

Dear Mauricio,

This helped a lot. Thank you also for the literature recommendation.

Kind regards,

Manuel 

[Philip]

Dear Manuel,

thank you for making the post so detailed and descriptive to discuss.

Do you already have an answer to this?

> If I decide to add gender as a grouping variable, would I add this to the measurement invariance analysis as well?

Thanks and regards,

Philip

[Terrence Jorgensen]

Do you already have an answer to this?

> If I decide to add gender as a grouping variable, would I add this to the measurement invariance analysis as well?

(Changes in) latent means are only comparable across groups when measurement-invariance constraints are in place.  And it is best to test the validity of those constraints before drawing conclusions.

[sura liu]

Hi Mauricio,

Thank you very much for your posts. As a psych student, I learned a lot from this thread when conducting my own analysis. Though, I still have some questions, could you possibly help?

I wondered: 

1. Why should I include std.lv = T? I thought the right way is to fix the first loading to 1.0 and freely estimate others?
2. I excluded estimator = “MLR” because it gave me a warning, but why? When should I use it? (I read the function description and still didn’t understand)

For a brief background, I'm estimating the level and slope of well-being across 3 waves of data, where well-being was measured by a group of same 5 manifest variables at each wave. I encountered the same lv variances are negative error as Manuel, but adding std.lv = T solved it.

I'd greatly appreciate it if you could help clarify!

Tingshu

On Wednesday, December 16, 2020 at 3:22:23 PM UTC-6 mauga...@gmail.com wrote:

Manuel

I think you are mixing contrainst based on identification methods. You are setting the marker variable identification , and using constraints like fixed variance.

Would make these changes, changes are highlighted

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ bos1_0 + bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ bos1_1 + bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ bos1_2 + bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

# same factor loadings for each item over time

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2

'

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

boredom_0 ~~ 1* boredom_0 

boredom_1 ~~ NA*boredom_1

boredom_2 ~~ NA*boredom_2 

boredom_0 ~0*1

boredom_1 ~NA*1

boredom_2 ~NA*1

'

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML", estimator = "MLR", std.lv=T)

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

############

Why are you building that function for the LRT between models? Becuase you are using the MLR estimator?

lavaan methods of nested model comparison do the correct adjust for MLR, MLM, DWLS, etc, estimator. WIth the lavTestLRT() for nested model comparison. If you go into details, you will some options for adjustsments based on the estimator

#########

Last, when you have a bagative LV variance, for which latent variable is it? If it for the Intecept or slope might be an indicationg that this is not a random parameter but should be treated as fix, which can be odd on the SEM, but its standard in growth curves in MLM.

On Wednesday, December 16, 2020 at 6:00:26 PM UTC+1 mailsc...@gmail.com wrote:

Dear Pat,

You are absolutely right! I am getting there…

I changed my code accordingly to establish longitudinal strong measurement invariance as follows:

# exploring measurement invariance over time

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ lambda0*bos1_1 + bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ lambda0*bos1_2 + bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# unique variances

bos1_0 ~~ bos1_0

bos2_0 ~~ bos2_0

bos3_0 ~~ bos3_0

bos4_0 ~~ bos4_0

bos5_0 ~~ bos5_0

bos6_0 ~~ bos6_0

bos1_1 ~~ bos1_1

bos2_1 ~~ bos2_1

bos3_1 ~~ bos3_1

bos4_1 ~~ bos4_1

bos5_1 ~~ bos5_1

bos6_1 ~~ bos6_1

bos1_2 ~~ bos1_2

bos2_2 ~~ bos2_2

bos3_2 ~~ bos3_2

bos4_2 ~~ bos4_2

bos5_2 ~~ bos5_2

bos6_2 ~~ bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

# same factor loadings for each item over time

weak_boredom <- ' 

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# unique variances and covariances

bos1_0 ~~ bos1_0

bos2_0 ~~ bos2_0

bos3_0 ~~ bos3_0

bos4_0 ~~ bos4_0

bos5_0 ~~ bos5_0

bos6_0 ~~ bos6_0

bos1_1 ~~ bos1_1

bos2_1 ~~ bos2_1

bos3_1 ~~ bos3_1

bos4_1 ~~ bos4_1

bos5_1 ~~ bos5_1

bos6_1 ~~ bos6_1

bos1_2 ~~ bos1_2

bos2_2 ~~ bos2_2

bos3_2 ~~ bos3_2

bos4_2 ~~ bos4_2

bos5_2 ~~ bos5_2

bos6_2 ~~ bos6_2

bos1_0 ~~ bos1_1 + bos1_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML", estimator = "MLR")

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

Unfortunately, if I am comparing the models now using the code from Steinmetz (2014) the first p-value is NA:

# corrected chi-square differences test

  CorrDiff <- function(fit.baseline, fit.restricted) # basline = model with less degrees of freedom

  {

  Corrc2Base <- fitMeasures(fit.baseline)[[6]] # corrected chi-square value

  MLc2Base <- fitMeasures(fit.baseline)[[3]] # uncorrected chi-square value

  dfBase <- fitMeasures(fit.baseline)[[4]] # degrees of freedom

  

  Corrc2Rest <- fitMeasures(fit.restricted)[[6]]

  MLc2Rest <- fitMeasures(fit.restricted)[[3]]

  dfRest <- fitMeasures(fit.restricted)[[4]]

  

  c0 <- MLc2Rest / Corrc2Rest

  c1 <- MLc2Base / Corrc2Base

  cd <- (dfRest*c0-dfBase*c1)/(dfRest-dfBase)

  CorrDiff <- (MLc2Rest-MLc2Base)/cd # corrected difference test statistic T_d^*

  dfDiff <- dfRest - dfBase

  pDiff <- 1-pchisq(CorrDiff, dfDiff)

  result <- list(Corr.chisq.difference = CorrDiff,

                 df.difference=dfDiff, p.value=pDiff)

  return(result)

  }

# model comparison

CorrDiff(fit_configural_boredom, fit_weak_boredom)

CorrDiff(fit_weak_boredom, fit_strong_boredom)

Warning:

In pchisq(CorrDiff, dfDiff) : NaNs produced

Furthermore, if I am then using the sem or growth function to calculate my growth model using the longitudinal scalar invariance model with the addition by Mauricio I get negative lv variances:

strong_boredom_growth <- '

# measurement model

boredom_0 =~ lambda0*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0 + lambda6*bos6_0

boredom_1 =~ lambda0*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1 + lambda6*bos6_1

boredom_2 =~ lambda0*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2 + lambda6*bos6_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 0*1

boredom_2 ~ 0*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ 0*boredom_1

boredom_0 ~~ 0*boredom_2

boredom_1 ~~ 0*boredom_2

# adding the growth model

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

i ~ 1

s ~ 1

'

# model estimation

fit_strong_boredom_growth <- sem(strong_boredom_growth, data=df1, missing = "FIML", estimator = "MLR")

Warning:

In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

Thanks,

Manuel

mal...@malonequantitative.com schrieb am Montag, 14. Dezember 2020 um 16:20:03 UTC+1:

Manuel,

Interleaved:

On Mon, Dec 14, 2020 at 5:44 AM Manuel <mailsc...@gmail.com> wrote:

Dear Pat,

As commented in the code (in the note directly above that line of code) this is done to force this factor loading to be free. I therefore pre-multiply it with NA, as a hint to Lavaan that the value of this parameter is still unknown.

I understand that, but you shouldn't need to do *both* that *and* lambda1*bos1_0 in the same measurement equation. The latter makes the former redundant.

I tried the growth and sem function and attached the two outputs as pdf. Using the sem function without i ~ 1 and s ~ 1 intercept and slope would be 0.

I find it a bit spooky that the results differ when using different functions. Are you certain that I can trust the results of the sem function?

That's why I said you'd need to specify those for the sem() wrapper. As far as which to trust, it looks like there's a slight shift in the observed-variable handling--the intercepts all differ by .018, which is reflected in a difference of .018 in the other direction in the intercept of i. The model appears to be otherwise identical, and the fit statistics are identical. The difference appears to be an arbitrary artifact that won't affect your interpretation.

Pat

 

Thanks,

Manuel

mal...@malonequantitative.com schrieb am Freitag, 11. Dezember 2020 um 17:33:05 UTC+1:

I haven't been following this thread closely, but what's the benefit of the first right-hand-side term in

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 +

(and similar for _1 and _2). I can't figure out what it would do, versus

boredom_0 =~ lambda1*bos1_0 +

Also, I'd recommend using the sem() wrapper instead of the growth() wrapper, since growth() is based on, as far as I know, growth of observed variables, not latents. You specify so much else in your model, I think sem() would work fine and be less likely to run into unintended defaults. Though you may need to add:

i ~ 1

s ~ 1

Pat

On Fri, Dec 11, 2020 at 11:17 AM Manuel <mailsc...@gmail.com> wrote:

Unfortunately changing means to 0 still leads resulted in the error. Adding missing = "FIML" solved it though.

Just to double check I am doing it right, here is my code:

# 2. Define Models --------------------------------------------------------

# exploring measurement invariance over time

library(lavaan)

# * 2.1 Configural Invariance ---------------------------------------------

# note: to force this factor loading to be free, we pre-multiply it with NA, as a hint to lavaan that the value of this parameter is still unknown.

configural_boredom <- '

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# unique variances

bos1_0 ~~ bos1_0

bos2_0 ~~ bos2_0

bos3_0 ~~ bos3_0

bos4_0 ~~ bos4_0

bos5_0 ~~ bos5_0

bos6_0 ~~ bos6_0

bos1_1 ~~ bos1_1

bos2_1 ~~ bos2_1

bos3_1 ~~ bos3_1

bos4_1 ~~ bos4_1

bos5_1 ~~ bos5_1

bos6_1 ~~ bos6_1

bos1_2 ~~ bos1_2

bos2_2 ~~ bos2_2

bos3_2 ~~ bos3_2

bos4_2 ~~ bos4_2

bos5_2 ~~ bos5_2

bos6_2 ~~ bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_configural_boredom <- cfa(configural_boredom, data = df1, missing = "FIML")

summary(fit_configural_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.2 Weak Invariance ---------------------------------------------------

# model specification (additional constraints on latent factors)

weak_boredom <- ' 

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

# intercepts

bos1_0 ~ i1*1

bos2_0 ~ 1

bos3_0 ~ 1

bos4_0 ~ 1

bos5_0 ~ 1

bos6_0 ~ 1

bos1_1 ~ i1*1

bos2_1 ~ 1

bos3_1 ~ 1

bos4_1 ~ 1

bos5_1 ~ 1

bos6_1 ~ 1

bos1_2 ~ i1*1

bos2_2 ~ 1

bos3_2 ~ 1

bos4_2 ~ 1

bos5_2 ~ 1

bos6_2 ~ 1

# unique variances and covariances

bos1_0 ~~ bos1_0

bos2_0 ~~ bos2_0

bos3_0 ~~ bos3_0

bos4_0 ~~ bos4_0

bos5_0 ~~ bos5_0

bos6_0 ~~ bos6_0

bos1_1 ~~ bos1_1

bos2_1 ~~ bos2_1

bos3_1 ~~ bos3_1

bos4_1 ~~ bos4_1

bos5_1 ~~ bos5_1

bos6_1 ~~ bos6_1

bos1_2 ~~ bos1_2

bos2_2 ~~ bos2_2

bos3_2 ~~ bos3_2

bos4_2 ~~ bos4_2

bos5_2 ~~ bos5_2

bos6_2 ~~ bos6_2

bos1_0 ~~ bos1_1 + bos1_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_weak_boredom <- cfa(weak_boredom, data = df1, missing = "FIML")

summary(fit_weak_boredom, fit.measures = TRUE, standardized = TRUE)

# * 2.3 Strong Invariance -------------------------------------------------

# model specification (additional constraints on manifest intercepts)

strong_boredom <- '

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 1

boredom_2 ~ 1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ boredom_1

boredom_0 ~~ boredom_2

boredom_1 ~~ boredom_2

'

# model estimation

fit_strong_boredom <- cfa(strong_boredom, data = df1, missing = "FIML")

summary(fit_strong_boredom, fit.measures = TRUE, standardized = TRUE)

# 3. Model Comparison ----------------------------------------------------

# * 3.1 Steinmetz, 2014 ---------------------------------------------------

# corrected chi-square differences test

CorrDiff <- function(fit.baseline, fit.restricted) # basline = model with less degrees of freedom

{

  Corrc2Base <- fitMeasures(fit.baseline)[[6]] # corrected chi-square value

  MLc2Base <- fitMeasures(fit.baseline)[[3]] # uncorrected chi-square value

  dfBase <- fitMeasures(fit.baseline)[[4]] # degrees of freedom

  

  Corrc2Rest <- fitMeasures(fit.restricted)[[6]]

  MLc2Rest <- fitMeasures(fit.restricted)[[3]]

  dfRest <- fitMeasures(fit.restricted)[[4]]

  

  c0 <- MLc2Rest / Corrc2Rest

  c1 <- MLc2Base / Corrc2Base

  cd <- (dfRest*c0-dfBase*c1)/(dfRest-dfBase)

  CorrDiff <- (MLc2Rest-MLc2Base)/cd # corrected difference test statistic T_d^*

  dfDiff <- dfRest - dfBase

  pDiff <- 1-pchisq(CorrDiff, dfDiff)

  result <- list(Corr.chisq.difference = CorrDiff,

                 df.difference=dfDiff, p.value=pDiff)

  return(result)

}

# model comparison

CorrDiff(fit_configural_boredom, fit_weak_boredom)

CorrDiff(fit_weak_boredom, fit_strong_boredom)

CorrDiff(fit_strong_boredom, fit_strict_boredom)

# 3.2 PennState Alternative -----------------------------------------------

# according to: https://quantdev.ssri.psu.edu/tutorials/intro-basics-longitudinal-measurement-invariance

# compare model fit statistics

round(cbind(configural=inspect(fit_configural_boredom, 'fit.measures'), 

            weak=inspect(fit_weak_boredom, 'fit.measures'),

            strong=inspect(fit_strong_boredom, 'fit.measures'),

            strict=inspect(fit_strict_boredom, 'fit.measures')),3)

# chi-square difference test for nested models 

anova(fit_configural_boredom, fit_weak_boredom)

anova(fit_weak_boredom, fit_strong_boredom)

anova(fit_strong_boredom, fit_strict_boredom)

# 4. Fitting Boredom Growth in Lavaan using Strong Invariance -------------

# model specification (first order factor means fixed to 0)

strong_boredom_growth <- '

# measurement model

boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

bos2_0 ~~ bos2_1 + bos2_2

bos3_0 ~~ bos3_1 + bos3_2

bos4_0 ~~ bos4_1 + bos4_2

bos5_0 ~~ bos5_1 + bos5_2

bos6_0 ~~ bos6_1 + bos6_2

# latent variables means

boredom_0 ~ 0*1

boredom_1 ~ 0*1

boredom_2 ~ 0*1

# latent variables variances and covariances

boredom_0 ~~ 1*boredom_0

boredom_1 ~~ boredom_1

boredom_2 ~~ boredom_2

boredom_0 ~~ 0*boredom_1

boredom_0 ~~ 0*boredom_2

boredom_1 ~~ 0*boredom_2

# adding the growth model

i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

'

# model estimation

fit_strong_boredom_growth <- growth(strong_boredom_growth, data=df1, missing = "FIML")

summary(fit_strong_boredom_growth, fit.measures = TRUE, standardized = TRUE)

# * 4.1 Adding Boredom Growth Parameters to Dataframe ---------------------

idx <- lavInspect(fit_strong_boredom_growth, "case.idx")

fscores <- lavPredict(fit_strong_boredom_growth, assemble = TRUE)

## loop over factors

for (fs in colnames(fscores)) {

  df1[idx, fs] <- fscores[, fs]

}

mauga...@gmail.com schrieb am Donnerstag, 10. Dezember 2020 um 11:01:00 UTC+1:

Manuel

Change these ines of code

You want the first order factor means fix to 0, so the mean change is represented by the intercepte and slope

Want to have the facto relation represented by the growth curve

 # latent variables means

 boredom_0 ~ 0*1

 boredom_1 ~ 0*1

 boredom_2 ~ 0*1

 # latent variables variances and covariances

 boredom_0 ~~ 1*boredom_0

 boredom_1 ~~ boredom_1

 boredom_2 ~~ boredom_2

 boredom_0 ~~ 0*boredom_1

 boredom_0 ~~ 0*boredom_2

 boredom_1 ~~ 0*boredom_2

On Thursday, December 10, 2020 at 8:58:10 AM UTC+1 mailsc...@gmail.com wrote:

Dear Mauricio,

I see. Then I misunderstood you saying I should "use the factors from the Strong invariance model for the growth curve". It also felt unusual.

This is my model for longitudinal scalar invariance and I added the growth curve at the end:

strong_boredom_growth <- '

 # measurement model

 boredom_0 =~ NA*bos1_0 + lambda1*bos1_0 + lambda2*bos2_0 + lambda3*bos3_0 + lambda4*bos4_0 + lambda5*bos5_0+ lambda6*bos6_0

 boredom_1 =~ NA*bos1_1 + lambda1*bos1_1 + lambda2*bos2_1 + lambda3*bos3_1 + lambda4*bos4_1 + lambda5*bos5_1+ lambda6*bos6_1

 boredom_2 =~ NA*bos1_2 + lambda1*bos1_2 + lambda2*bos2_2 + lambda3*bos3_2 + lambda4*bos4_2 + lambda5*bos5_2+ lambda6*bos6_2

 bos2_0 ~~ bos2_1 + bos2_2

 bos3_0 ~~ bos3_1 + bos3_2

 bos4_0 ~~ bos4_1 + bos4_2

 bos5_0 ~~ bos5_1 + bos5_2

 bos6_0 ~~ bos6_1 + bos6_2 

 # latent variables means

 boredom_0 ~ 0*1

 boredom_1 ~ 1

 boredom_2 ~ 1

 # latent variables variances and covariances

 boredom_0 ~~ 1*boredom_0

 boredom_1 ~~ boredom_1

 boredom_2 ~~ boredom_2

 boredom_0 ~~ boredom_1

 boredom_0 ~~ boredom_2

 boredom_1 ~~ boredom_2

 # adding the growth model

 i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

 s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope

'

# model estimation

boredom_growth_model = growth(strong_boredom_growth, data=df1)

summary(boredom_growth_model)

Unfortunately, this again results in an error:

1: In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats,  :

  lavaan WARNING:

    The variance-covariance matrix of the estimated parameters (vcov)

    does not appear to be positive definite! The smallest eigenvalue

    (= 1.853961e-18) is close to zero. This may be a symptom that the

    model is not identified.

2: In lav_object_post_check(object) :

  lavaan WARNING: some estimated lv variances are negative

I looked at the summary anyways and attached it as pdf, if this is of any help.

Kind regards,

Manuel

mauga...@gmail.com schrieb am Donnerstag, 10. Dezember 2020 um 00:49:39 UTC+1:

Manuel

It is not recommended for you to add the factor scores as measured variables, and then run the growth model with them. This goes back into the porblem of factor score indeteminancy and other issues of factor score variability and method.

The recommend steps would be t estimate the growth curve on top of the longitudinal scalar invariance model. This would required you to add 2 pieces to the measurement inavriance model

- The intercept and slope factors define by the first order factors, as you have it

- fix the first order factors means to 0. This way the mean change over time would be only define by the growth curve

Hope this helps

On Wednesday, December 9, 2020 at 4:54:55 PM UTC+1 mal...@malonequantitative.com wrote:

Manuel,

It sounds like a different usage of "factor."

What do you get from class(df1$boredom_0) ? It should be either numeric (double) or ordered, not factor. Factor in that sense would mean a nominal variable, which can't be used as an endogenous variable in lavaan. Using as.numeric will give a value for each unique value of a nominal variable.

lavPredict() should be returning class numeric. I'm not sure where they're getting converted to class factor.

Pat

On Wed, Dec 9, 2020 at 10:43 AM Manuel <mailsc...@gmail.com> wrote:

Dear Magua, dear Patrick,

Thank you for your help!

I have managed to establish scalar longitudinal measurement invariance for boredom and added the scalar boredom factor scores to my data frame using this code:

idx <- lavInspect(fit_strong_boredom, "case.idx")

fscores <- lavPredict(fit_strong_boredom, assemble = TRUE)

## loop over factors

for (fs in colnames(fscores)) {

  df1[idx, fs] <- fscores[, fs]

}

Now I am trying to extract the slope parameters using the boredom factor scores:

boredom_growth <- '

    i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2 # intercept

    s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2 # slope'

growthCurveModel = growth(boredom_growth, data=df1)

But this again results in an error:

Error in lav_data_full(data = data, group = group, cluster = cluster,  : 

  lavaan ERROR: unordered factor(s) detected; make them numeric or ordered: boredom_0 boredom_1 boredom_2

The factor scores range from -0,000154379130794116 to 2,78171374926366. Using as.numeric(df1$boredom_0) strangely results in values from 1 to 1186.

Do you have any idea how I can solve this?

Thanks again,

Manuel

mal...@malonequantitative.com schrieb am Mittwoch, 2. Dezember 2020 um 16:40:33 UTC+1:

Manuel,

I suspect you're running into issues with the factor intercepts.

I'm not sure whether the growth() function can do this, but the sem()
function could. Then if you specify scalar invariance on the boredom
factors as Mauricio mentioned, then

boredom_0 ~ 0*1
boredom_1 ~ 0*1
boredom_2 ~ 0*1

i ~ 1
s ~ 1

That should get you on the right track.

On Wed, Dec 2, 2020 at 5:07 AM Manuel <mailsc...@gmail.com> wrote:
>
> Dear Terrence,
>
> Thank you very much for your answer. It helped a lot.
>
> Is there a way to run a double latent growth model in Lavaan?
> I would like to use boredom factor scores instead of manifest boredom means. In my mind it should look like this:
>
> latent_boredom_growth <- '
> boredom_0 =~ bos1_0 + bos2_0 + bos3_0 + bos4_0 + bos5_0 + bos6_0; #latent boredom

> boredom_1 =~ bos1_1 + bos2_1 + bos3_1 + bos4_1 + bos5_1 + bos6_1;

> boredom_2 =~ bos1_2 + bos2_2 + bos3_2 + bos4_2 + bos5_2 + bos6_2;

>
> i =~ 1*boredom_0 + 1*boredom_1 + 1*boredom_2; # intercept
> s =~ 0*boredom_0 + 1*boredom_1 + 2*boredom_2; # slope

> '
> boredom_growth_model = growth(latent_boredom_growth, data=lavan)
>
> But this results in an error:

>
> In lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, :
> lavaan WARNING:

> Could not compute standard errors! The information matrix could
> not be inverted. This may be a symptom that the model is not
> identified.
>
> Thanks again and kind regards,
> Manuel
> Terrence Jorgensen schrieb am Freitag, 13. November 2020 um 11:59:12 UTC+1:
>>>
>>> Am I doing this right?
>>
>> Looks right.
>>>
>>> How do I interpret the output?
>>
>> https://lavaan.ugent.be/tutorial/growth.html
>> Try Googling for more tutorials about growth curves in lavaan, like this one: https://quantdev.ssri.psu.edu/tutorials/growth-modeling-chapter-10-growth-models-nonlinearity-time
>>
>>> Why are all intercepts 0?
>>
>> They aren't. Only the indicators have intercepts == 0 because their means are reproduces as functions of the growth parameters. The growth parameter means are merely very close to zero (try adding the argument nd=7 for the summary() to print up to 7 digits of precision). Slopes nearly zero seems very consistent with the dark blue line in your plot being almost flat. The plot indicates your (unconditional) average intercept should be somewhere between 2 and 3; however, you have several predictors of the growth factors. Thus, the intercepts are the mean level and change when all those predictors == 0. Does that make sense for your data?
>>>
>>> What does it mean that intercept, slope, and quadratic slope are not significant?
>>
>> It implies those parameters are consistent with a population parameter == 0.

>>
>> 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
>>
>

> --
> You received this message because you are subscribed to the Google Groups "lavaan" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/2ac441be-e959-4a3b-b72e-f72997d94bebn%40googlegroups.com.

--
Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

--

You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/5f9bb835-84c4-4390-b659-33e6fef4d105n%40googlegroups.com.

--

Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

--

You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/8c2bc5a4-e34e-4526-a29b-0b662e94b644n%40googlegroups.com.

--

Patrick S. Malone, Ph.D., Malone Quantitative
NEW Service Models: http://malonequantitative.com

He/Him/His

--

You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/8668c7e6-8035-4d1d-8d8e-13034d1d61cbn%40googlegroups.com.

[Mauricio Garnier-Villarreal]

Tingshu

Here are my comments

####

Why should I include std.lv = T? I thought the right way is to fix the first loading to 1.0 and freely estimate others?

####

There are a couple different ways to identify the factors, as you mention the most common (default) method is marker variable (fixing a factor loading to 1), the second most common is fixed factor variance (std.lv=T). The identification method does not affect the model fit, or standardized parameters. But it affects the metric of interpretation for the unstandardized paramaters. With marker variable, the factor mean and variances are in function of the items that has the factor loading to 1. Whith fixed factor variance, the metric and interpretation of the factors is as an standardized metric (mean = 0, variance = 1). I strongly prefer fix factor variance, as the metric does not depend on a single variable, but still some people prefer marker variable.

####

I excluded estimator = “MLR” because it gave me a warning, but why? When should I use it? (I read the function description and still didn’t understand)

####

MLR adjusts the chi-square and standard errors for non-normal data distribituions (it doesnt change the parameter estimates). So, it is not a required to use it, in case you dont need to adjust for non-normality, using the default ML is fine. Now depending on what the warning is, might represent more serious problems. But is it works well with ML I dont think there is mayor concerns

####

For a brief background, I'm estimating the level and slope of well-being across 3 waves of data, where well-being was measured by a group of same 5 manifest variables at each wave. I encountered the same lv variances are negative error as Manuel, but adding std.lv = T solved it.

####

With 3 time points you are at the minimum number of time points, so be aware of possible limitations becasue of this. The latent variances were negative because they were defined in function of the choosen indicator, when you changed to std.lv=T the factor variance was in function of the standardized metric. Another possible solution for that, with marker variable is to use another variane as the one with the fixed factor loading. But again, I prefer and recommend std.lv=T

Hope this helps

[Zachary Salander]

Hi Mauricio,

This thread has been incredibly helpful with my own analyses, so thank you!

If I have found the best fitting LGM for multiple constructs and I want to explore parallel processes, would it be as simple as just combining each construct's code (see below) or are there further adjustments/edits that are recommended? Any/All input is greatly appreciated!

Code for Parallel Processes (cognitive competence=CC, value=VAL, anxiety=ANX):

predict <- '
##### MEASUREMENT MODEL: #####
CC_1 =~ lambda1*CC1 + lambda2*CC2 + lambda3*CC4 + lambda4*CC5
CC_2 =~ lambda1*CC1_2 + lambda2*CC2_2 + lambda3*CC4_2 + lambda4*CC5_2
CC_3 =~ lambda1*CC1_3 + lambda2*CC2_3 + lambda3*CC4_3 + lambda4*CC5_3

VAL_1 =~ lambda1.2*VAL1 + lambda2.2*VAL2 + lambda3.2*VAL3 + lambda4.2*VAL4 + lambda5.2*VAL5
VAL_2 =~ lambda1.2*VAL1_2 + lambda2.2*VAL2_2 + lambda3.2*VAL3_2 + lambda4.2*VAL4_2 + lambda5.2*VAL5_2
VAL_3 =~ lambda1.2*VAL1_3 + lambda2.2*VAL2_3 + lambda3.2*VAL3_3 + lambda4.2*VAL4_3 + lambda5.2*VAL5_3

ANX_1 =~ lambda1.3*ANX1 + lambda2.3*ANX3 + lambda3.3*ANX4 + lambda4.3*ANX6 + lambda5.3*ANX7 + lambda6.3*ANX8
ANX_2 =~ lambda1.3*ANX1_2 + lambda2.3*ANX3_2 + lambda3.3*ANX4_2 + lambda4.3*ANX6_2 + lambda5.3*ANX7_2 + lambda6.3*ANX8_2
ANX_3 =~ lambda1.3*ANX1_3 + lambda2.3*ANX3_3 + lambda3.3*ANX4_3 + lambda4.3*ANX6_3 + lambda5.3*ANX7_3 + lambda6.3*ANX8_3


##### INTERCEPTS: #####
CC1 ~ i1*1
CC2 ~ i2*1
CC4 ~ i3*1
CC5 ~ i4*1
CC1_2 ~ i5*1
CC2_2 ~ i2*1
CC4_2 ~ i7*1
CC5_2 ~ i4*1
CC1_3 ~ i9*1
CC2_3 ~ i2*1
CC4_3 ~ i11*1
CC5_3 ~ i4*1

VAL1 ~ i1.2*1
VAL2 ~ i2.2*1
VAL3 ~ i3.2*1
VAL4 ~ i4.2*1
VAL5 ~ i5.2*1
VAL1_2 ~ i1.2*1
VAL2_2 ~ i2.2*1
VAL3_2 ~ i3.2*1
VAL4_2 ~ i4.2*1
VAL5_2 ~ i10.2*1
VAL1_3 ~ i1.2*1
VAL2_3 ~ i2.2*1
VAL3_3 ~ i3.2*1
VAL4_3 ~ i4.2*1
VAL5_3 ~ i15.2*1

ANX1 ~ i1.3*1
ANX3 ~ i2.3*1
ANX4 ~ i3.3*1
ANX6 ~ i4.3*1
ANX7 ~ i5.3*1
ANX8 ~ i6.3*1
ANX1_2 ~ i1.3*1
ANX3_2 ~ i8.3*1
ANX4_2 ~ i9.3*1
ANX6_2 ~ i10.3*1
ANX7_2 ~ i11.3*1
ANX8_2 ~ i6.3*1
ANX1_3 ~ i1.3*1
ANX3_3 ~ i14.3*1
ANX4_3 ~ i15.3*1
ANX6_3 ~ i16.3*1
ANX7_3 ~ i17.3*1
ANX8_3 ~ i6.3*1


##### UNIQUE VARIANCES & COVARIANCES: #####
CC1 ~~ CC1 + CC1_2 + CC1_3
CC2 ~~ CC2 + CC2_2 + CC2_3
CC4 ~~ CC4 + CC4_2 + CC4_3
CC5 ~~ CC5 + CC5_2 + CC5_3
CC1_2 ~~ CC1_2 + CC1_3
CC2_2 ~~ CC2_2 + CC2_3
CC4_2 ~~ CC4_2 + CC4_3
CC5_2 ~~ CC5_2 + CC5_3
CC1_3 ~~ CC1_3
CC2_3 ~~ CC2_3
CC4_3 ~~ CC4_3
CC5_3 ~~ CC5_3

VAL1 ~~ VAL1 + VAL1_2 + VAL1_3
VAL2 ~~ VAL2 + VAL2_2 + VAL2_3
VAL3 ~~ VAL3 + VAL3_2 + VAL3_3
VAL4 ~~ VAL4 + VAL4_2 + VAL4_3
VAL5 ~~ VAL5 + VAL5_2 + VAL5_3
VAL1_2 ~~ VAL1_2 + VAL1_3
VAL2_2 ~~ VAL2_2 + VAL2_3
VAL3_2 ~~ VAL3_2 + VAL3_3
VAL4_2 ~~ VAL4_2 + VAL4_3
VAL5_2 ~~ VAL5_2 + VAL5_3
VAL1_3 ~~ VAL1_3
VAL2_3 ~~ VAL2_3
VAL3_3 ~~ VAL3_3
VAL4_3 ~~ VAL4_3
VAL5_3 ~~ VAL5_3

ANX1 ~~ ANX1 + ANX1_2 + ANX1_3
ANX3 ~~ ANX3 + ANX3_2 + ANX3_3
ANX4 ~~ ANX4 + ANX4_2 + ANX4_3
ANX6 ~~ ANX6 + ANX6_2 + ANX6_3
ANX7 ~~ ANX7 + ANX7_2 + ANX7_3
ANX8 ~~ ANX8 + ANX8_2 + ANX8_3
ANX1_2 ~~ ANX1_2 + ANX1_3
ANX3_2 ~~ ANX3_2 + ANX3_3
ANX4_2 ~~ ANX4_2 + ANX4_3
ANX6_2 ~~ ANX6_2 + ANX6_3
ANX7_2 ~~ ANX7_2 + ANX7_3
ANX8_2 ~~ ANX8_2 + ANX8_3
ANX1_3 ~~ ANX1_3
ANX3_3 ~~ ANX3_3
ANX4_3 ~~ ANX4_3
ANX6_3 ~~ ANX6_3
ANX7_3 ~~ ANX7_3
ANX8_3 ~~ ANX8_3


##### LATENT VARIABLE MEANS: #####
CC_1 ~ 0*1
CC_2 ~ 0*1
CC_3 ~ 0*1
VAL_1 ~ 0*1
VAL_2 ~ 0*1
VAL_3 ~ 0*1
ANX_1 ~ 0*1
ANX_2 ~ 0*1
ANX_3 ~ 0*1


##### LATENT VARIABLE VARIANCES & COVARIANCES: #####
CC_1 ~~ 1*CC_1
CC_2 ~~ NA*CC_2
CC_3 ~~ NA*CC_3
VAL_1 ~~ 1*VAL_1
VAL_2 ~~ NA*VAL_2
VAL_3 ~~ NA*VAL_3
ANX_1 ~~ 1*ANX_1
ANX_2 ~~ NA*ANX_2
ANX_3 ~~ NA*ANX_3


##### GROWTH MODEL: #####
iCC =~ 1*CC_1 + 1*CC_2 + 1*CC_3
sCC =~ 0*CC_1 + 1*CC_2 + 2*CC_3
iCC ~ NA*1
sCC ~ NA*1
iCC ~~ NA*iCC
sCC ~~ NA*sCC

iVAL =~ 1*VAL_1 + 1*VAL_2 + 1*VAL_3
sVAL =~ 0*VAL_1 + 1*VAL_2 + 2*VAL_3
iVAL ~ NA*1
sVAL ~ NA*1
iVAL ~~ NA*iVAL
sVAL ~~ NA*sVAL

iANX =~ 1*ANX_1 + 1*ANX_2 + 1*ANX_3
sANX =~ 0*ANX_1 + 1*ANX_2 + 2*ANX_3
iANX ~ NA*1
sANX ~ NA*1
iANX ~~ NA*iANX
sANX ~~ NA*sANX
'

Best,

Zach

[Mauricio Garnier-Villarreal]

Dear Zach

I am glad this thread has been helpful

Yes, if you have defined growth models for multiple factors, you can only paste all the syntaxes together to estimate the growth models together (be careful with the labels). I think by default lavaan will estimate the correlations between factors. 

If you want to estimate regressions between, you would want to add a new section with regressions between random intercepts/slopes

take care