Three-level longitudinal SEM--is it possible?

[Courtney McDermott]

I'm working with three level data (two level with repeated measurements): Emotional Self Awareness (ESA) measures were taken at three time points, from students nested within schools. Five of the schools were given an intervention, five were control. I've already analysed the data using a mixed model in the lme4 package, but I'm trying to get SEM to match those results. The variables in the dataset include:

ESA1

ESA2

ESA3 (three time points that measured emotional self-awareness)

trtGrp (treatment group)

School (cluster group)

idChild (child ID)

For the mixed model, I ran a random intercept model, but ONLY at the school level (and fixed slope).  I defined a compound symmetry structure (at the student level).  The outcome, ESA, is regressed on Time, treatment group, and the interaction of Time*treatment group. The interactions are the coefficients of interest, and these are what I'm having trouble matching to an SEM.

My mixed model syntax (data is in long format) is as follows:

test <- lme(ESA ~ Time + trtGrp + trtGrp*Time, 

            random = list(School = ~1, idChild = ~1), 

            weights =varIdent(School),

            corr = corCompSymm(form= ~Time),

            data=DataESA, na.action="na.omit",

            method = "REML")

Using lavaan, I'm able to get the variance/covariance at both the child and school levels in an empty model to match the mixed model using this syntax:

## empty model, with correct variance/covariance

# data is in wide format

model = '  

level: 1

#compound symmetry var/covar structure

ESA1 ~~ v*ESA1

ESA2 ~~ v*ESA2

ESA3 ~~ v*ESA3

ESA1 ~~ c*ESA2 + c*ESA3

ESA2 ~~ c*ESA3

level: 2

#unstructured covar

ESA1 ~~ ESA2 + ESA3

ESA2 ~~ ESA3

'

fit = sem(model, data = DataESA_WIDE, cluster = "School")

summary(fit)

However, the problem comes when I add the predictors/regressions. From articles/tutorials I've seen online, it's easiest to use a latent intercept and slope for repeated-measures data. I have two slopes since the slope from T1 to T2 to T3 are not linear, so f1=the slope from t1 to t2, and f2 is the slope from t1 to t3 (essentially dummy coding). When I run a single-level model (model 2), the regression coefficients and SEs are identical to the mixed model, but the variance and covariance are inflated since it does not take into account the school clustering. When I attempt to use multilevel SEM (model 3), I can get the correct var/covar, but incorrect regression coefficients.

## Gives correct coefficients and SEs for regressions

#problem is that only one level, doesn't take into account school clustering

model2 = ' 

#created latent intercept (b1), latent slope for time 1 to time 2,

# and latent slope for time 1 to time 3 

b1 =~ 1*ESA1 + 1*ESA2 + 1*ESA1

f1 =~ 0*ESA1 + 1*ESA2 

f2 =~ 0*ESA1 + 1*ESA3 

#all fixed effects, since only random effect is intercept by cluster (school)

b1 ~ 1

b1 ~~ 0*b1

f1 ~~ 0*f1

f2 ~~ 0*f2

f1 ~~ 0*f2

b1 ~~ 0*f1

b1 ~~ 0*f2

#regressed on Treatment Group

b1 ~ trtGrp

f1 ~ trtGrp 

f2 ~ trtGrp

#compound symmetry

ESA1 ~~ v*ESA1

ESA2 ~~ v*ESA2

ESA3 ~~ v*ESA3

ESA1 ~~ c*ESA2 + c*ESA3

ESA2 ~~ c*ESA3

'

fit1 = sem(model1, data = DataESA_WIDE)

summary(fit1)

## Multilevel model, but incorrect coefficients/SEs; correct variance and covariance

model3 = ' 

level: 1

#latent intercept and slopes

b1 =~ 1*ESA1 + 1*ESA2 + 1*ESA3 

f1 =~ 0*ESA1 + 1*ESA2

f2 =~ 0*ESA1 + 1*ESA3

f1 ~ 1

f2 ~ 1

b1 ~ 1

#fixed effects at level 1 (no random intercept or slope)

f1 ~~ 0*f1

f2 ~~ 0*f2

f1 ~~ 0*f2

b1 ~~ 0*b1

b1 ~~ 0*f1 + 0*f2

#regressed on Treatment Group

b1 ~ trtGrp

f1 ~ trtGrp

f2 ~ trtGrp

#compound symmetry

ESA1 ~~ v*ESA1

ESA2 ~~ v*ESA2

ESA3 ~~ v*ESA3

ESA1 ~~ c*ESA2 + c*ESA3

ESA2 ~~ c*ESA3

level:2

#unstructured variance at level 2

ESA1 ~~ ESA1

ESA2 ~~ ESA2

ESA3 ~~ ESA3

ESA1 ~~ ESA2 + ESA3

ESA2 ~~ ESA3

'

fit1 = sem(model1, data = DataESA_WIDE, cluster = "School")

summary(fit1)

I'm new to structural equation modeling, and still fairly new to R, so I'm unable to figure out what I'm doing wrong. If anyone knows how to fix my code so that it reflects a random-intercept (ONLY for schools) model for repeated measures, with clustered data, I would really appreciate it. If it is not possible, please also let me know that too!