Longitudinal sem model - constrain initial time point variance to be equal to subsequent time point variances?

[Erik Ruzek]

Hi all,

I am trying to replicate models in Nestler & Humberg's (2023) article in the SEM journal. In it, they show how to estimate a series of longitudinal models in both mixed effects and SEM formulations in R (nlme and lavaan).

The model of interest is the first one, the random intercept - autoregressive panel model. The data can be downloaded at osf along with the R code. First is the mixed version in Stata:

mixed e i.time || id: , residuals(ar 1, t(time))

Next is the trickier sem version, and this is where I am having trouble. The issue is that in mixed, a stationarity constraint is imposed on the variances of the timepoint specific residuals. I cannot figure out how to do this in sem. My code is below:

reshape wide e time1vsrest, i(id) j(time)

sem (RI -> e1@1 e2@1 e3@1 e4@1 e5@1) /// random intercept
    (E1 -> e1@1) (E2 -> e2@1) (E3 -> e3@1) /// person-centered variables
    (E4 -> e4@1) (E5 -> e5@1) /// person-centered variables
    (E2 <- E1@c1) (E3 <- E2@c1) (E4 <- E3@c1) /// autoregressions
    (E5 <- E4@c1) , ///
    covstructure(e._OEn, zero) ///
    covstructure(_OEx, diagonal) ///
    covstructure(e._LEn, identity) ///
    covariance(RI*E1@0) ///
    noivstart technique(nr bhhh) ///
    nocnsreport nodescribe

The sem code above was modified from some code shared by Joseph Coveney in a 2021 thread.

Again, my main question is how do I fix the residual variance of E1 to be the same as E2-E5?

[Joseph Coveney]

. . . my main question is how do I fix the residual variance of E1 to be the same as E2-E5?

E1 is exogenous in your model and so it doesn't have a residual variance. You need to make it endogenous.

version 18.0

clear *

quietly import delimited IllustrativeData.txt, delimiter(space)

mixed e ibn.time, noconstant || id: , residuals(ar 1, t(time)) nolrtest nolog

quietly reshape wide e time1vsrest, i(id) j(time)

generate byte k = 1
sem ///
    (e1@1 e2@1 e3@1 e4@1 e5@1 <- RI) ///
    (e1@1 <- E1) (e2@1 <- E2) (e3@1 <- E3) (e4@1 <- E4) (e5@1 <- E5) ///
    (E1@0 <- k) (E2 <- E1@c1) (E3 <- E2@c1) (E4 <- E3@c1) (E5 <- E4@c1), ///
        covstructure(e._OEn, zero) covstructure(e._LEn, identity) ///
            noivstart ///
                nocnsreport nodescribe nolog
assert e(coverged)
mata: rowsum(st_matrix("e(ilog)") :!= 0) // Converged in eight iterations

exit

The tack shown above addresses your main question although you'll have to verify for yourself that the SEM is what you're after.

The link to the server directory containing the dataset and R code I think is here and that directly to the dataset itself is here.

[Joseph Coveney]

Actually, you don't need to explicitly generate the constant. You can render E1 endogenous with the implied constant, as follows:

sem ///
    (e1@1 e2@1 e3@1 e4@1 e5@1 <- RI) ///
    (e1@1 <- E1) (e2@1 <- E2) (e3@1 <- E3) (e4@1 <- E4) (e5@1 <- E5) ///
    (E1@0 <- ) (E2 <- E1@c1) (E3 <- E2@c1) (E4 <- E3@c1) (E5 <- E4@c1), ///
        covstructure(e._OEn, zero) covstructure(e._LEn, identity) ///
            noivstart ///
                nocnsreport nodescribe nolog

[Erik Ruzek]

Thank you so much, Joseph Coveney! When I look at the lavaan code in the osf directory for the paper, they impose some constraints on the estimation of the residual variances to impose stationarity. It is not just constraining the residual variances to be the same. The variances are further augmented by the autoregressive parameter, shown below, where v is the constraint on the E1 variance and vv are the constrained equal variances on E2-E5, and a is the autoregressive parameter label:

v == vv/(1-a*a)

I believe something like this is possible in Stata through the constraints syntax, but I will have to investigate more. If you have any ideas, please do suggest them!

[Joseph Coveney]

I believe something like this is possible in Stata through the constraints syntax . . .

Unfortunately, no. Stata's constraint cannot accept nonlinear constraints.

There is a workaround involving minbound that will allow you to impose the necessary nonlinear constraint and obtain point estimates that are identical to those with R's lavaan (and Stata's mixed above in #2). The workaround is illustrated in the code below:

version 18.0

clear *

quietly import delimited IllustrativeData.txt, delimiter(space)

quietly reshape wide e time1vsrest, i(id) j(time)

program define ri_apm, rclass
    version 18.0
    args a

    constraint define 1 [/]var(e.E2) = [/]var(e.E1) * (1 - `a' * `a')
    constraint define 2 [/]var(e.E3) = [/]var(e.E2)
    constraint define 3 [/]var(e.E4) = [/]var(e.E2)
    constraint define 4 [/]var(e.E5) = [/]var(e.E2)
    constraint define 5 [E2]E1 = `a'

    sem ///
        (e1@1 e2@1 e3@1 e4@1 e5@1 <- RI) ///
        (e1@1 <- E1) (e2@1 <- E2) (e3@1 <- E3) (e4@1 <- E4) (e5@1 <- E5) ///
        (E1@0 <- ) (E2 <- E1@rho) (E3 <- E2@rho) (E4 <- E3@rho) (E5 <- E4@rho),  ///
            constraints(1/5) covstructure(e._OEn, zero) ///
                noivstart

    return scalar fx = e(ll) * e(ll)
end

minbound ri_apm, range(0.2 0.3) trace

sem , nocnsreport

exit

Because of the constraint on the latent factors' serial regressions, the standard errors in the fitted model are a bit optimistic, and I suppose that you'll need to resort to bootstrapping for more realistic standard errors.

Allowing nonlinear constraints (and ad hoc parameters) for sem would be worth asking for in the Stata 19 wishlist thread inasmuch as MPlus allows them.

[Erik Ruzek]

Joseph Coveney, I am speechless. Thank you so very much for following up with this detailed explanation and for always being so helpful and thorough on this forum. I have learned a ton both about Stata coding and mixed effects and SEM from your posts through the years.