You are here

How to specify a random intercept cross-lagged panel model in OpenMx?

9 posts / 0 new
Last post
cjvanlissa's picture
Joined: 04/10/2019 - 12:43
How to specify a random intercept cross-lagged panel model in OpenMx?

Dear readers,

I am trying to specify an RI-CLPM in OpenMx. However, if I constrain the variance of the observed indicators to zero, the model implied covariance matrix is not positive definite - but these variances are supposed to be zero, as they are partialized into a time-invariant component and a time-variant component. E.g., the variance x1 <-> x1 should be split up into RIx <-> RIx and cx1 <-> cx1.

Can anyone help me correctly specify this model in OpenMx?

df <- data.frame(matrix(rnorm(500*8), ncol = 8))
names(df) <- paste0(rep(c("x", "y"), each = 4), 1:4)
xobs <- grep("^x", names(df), value = TRUE)
yobs <- grep("^y", names(df), value = TRUE)
xvars <- paste0("c", grep("^x", names(df), value = TRUE))
yvars <- paste0("c", grep("^y", names(df), value = TRUE))
xminone <- xvars[-1]
xminend <- xvars[-length(xvars)]
yminone <- yvars[-1]
yminend <- yvars[-length(yvars)]
              manifestVars=c(xobs, yobs),
              latentVars = c("RIx", "RIy", xvars, yvars),
              mxData(observed=df, type="raw"),
              mxPath(from="RIx", to=xobs, arrows=1, free=FALSE, values=1),
              mxPath(from="RIy", to=yobs, arrows=1, free=FALSE, values=1),
              mxPath(from=xvars, to=xobs, arrows=1, free=FALSE, values=1),
              mxPath(from=yvars, to=yobs, arrows=1, free=FALSE, values=1),
              # RI cors
              mxPath(from="RIy", to="RIx", arrows=2, free=TRUE, values=0),
              mxPath(from=c("RIy", "RIx"), arrows=2, free=TRUE, values=0),
              mxPath(from=xminend, to=xminone, arrows=1, free=TRUE, values=c(.1), labels="x_x"),
              mxPath(from=yminend, to=yminone, arrows=1, free=TRUE, values=c(.1), labels="y_y"),
              mxPath(from=xminend, to=yminone, arrows=1, free=TRUE, values=c(.1), labels="y_x"),
              mxPath(from=yminend, to=xminone, arrows=1, free=TRUE, values=c(.1), labels="x_y"),
              # Vars
              mxPath(from=xvars[1], arrows=2, free=TRUE, values=c(.1), labels = "vx1"),
              mxPath(from=yvars[1], arrows=2, free=TRUE, values=c(.1), labels = "vy1"),
              mxPath(from=xminone, arrows=2, free=TRUE, values=c(.1), labels = "vx"),
              mxPath(from=yminone, arrows=2, free=TRUE, values=c(.1), labels = "vy"),
              #mxPath(from=c(xobs, yobs), arrows=2, free=TRUE, values=c(.1)),
              mxPath(from=c(xobs, yobs), arrows=2, free=FALSE, values=0),
              # Covars
              mxPath(from=xvars[1], to = yvars[1], arrows=2, free=TRUE, values=0, labels = "rxy1"),
              mxPath(from=xminone, to = yminone, arrows=2, free=TRUE, values=0, labels = "rxy"),
              mxPath(from = 'one', to = c(xvars, yvars))
res <- mxRun(riclpm)
mhunter's picture
Joined: 07/31/2009 - 15:26
Forgot variances of RIx and RIy?

I ran


It looked like you set the variances of RIx and RIy to zero. So, I added

       mxPath(c('RIx', 'RIy'), arrows=2, values=1.5),

into your model and it ran, only giving a warning about the Hessian being non-positive definite (probably because the data were nonsense).

Does that help?

cjvanlissa's picture
Joined: 04/10/2019 - 12:43
Dear mhunter, thank you! My

Dear mhunter, thank you! My code already contained the line:

mxPath(from=c("RIy", "RIx"), arrows=2, free=TRUE, values=0)

If I set this to

mxPath(from=c("RIy", "RIx"), arrows=2, free=TRUE, values=1.5)

I indeed also get the error about the Hessian (also with the real data). Also, there are no estimates - instead, the starting values are returned:

Warning message:
In model 'riclpm' Optimizer returned a non-zero status code 5. The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED). 
> summary(res)
Summary of riclpm 
The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED). 
free parameters:
              name matrix row col Estimate  Std.Error A
1              x_x      A cx2 cx1      0.1         NA !
2              y_x      A cy2 cx1      0.1         NA !
3              x_y      A cx2 cy1      0.1         NA !
4              y_y      A cy2 cy1      0.1         NA !
5    riclpm.S[9,9]      S RIx RIx      1.5         NA !
6   riclpm.S[9,10]      S RIx RIy      0.1         NA !
7  riclpm.S[10,10]      S RIy RIy      1.5         NA !
8              vx1      S cx1 cx1      0.1         NA !
9               vx      S cx2 cx2      0.1         NA !
10            rxy1      S cx1 cy1      0.1         NA !
11             vy1      S cy1 cy1      0.1         NA !
12             rxy      S cx2 cy2      0.1         NA !
13              vy      S cy2 cy2      0.1         NA !
14  riclpm.M[1,11]      M   1 cx1      0.1 8601467937 !
15  riclpm.M[1,12]      M   1 cx2      0.1 8866184397 !
16  riclpm.M[1,13]      M   1 cx3      0.1 9126759905 !
17  riclpm.M[1,14]      M   1 cx4      0.1 7286763147 !
18  riclpm.M[1,15]      M   1 cy1      0.1 8601467226 !
19  riclpm.M[1,16]      M   1 cy2      0.1 8866192940 !
20  riclpm.M[1,17]      M   1 cy3      0.1 9126760586 !
21  riclpm.M[1,18]      M   1 cy4      0.1 7286752751 !
cjvanlissa's picture
Joined: 04/10/2019 - 12:43
Interestingly, if I re-run

Interestingly, if I re-run that model with mxTryHard, the results are identical to the lavaan solution. So... I guess this points to a problem of bad starting values? I still don't understand why the error about the Hessian appears, nor why the model doesn't converge with mxRun.

AdminRobK's picture
Joined: 01/24/2014 - 12:15
Starting values are not feasible

If I use mxOption(NULL,"Default optimizer","NPSOL") to switch optimizers to NPSOL, and then run the model, NPSOL returns status code 10 ("Starting values are not feasible"). Looking at the gradient with res$output$gradient reveals that some of the gradient elements are NA at the start values.

I'm not sure why there's no warning about the infeasible start when running the model with SLSQP, which is the on-load default optimizer. I think there ought to be such a warning. It's a useful diagnostic.

cjvanlissa's picture
Joined: 04/10/2019 - 12:43
Thank you, AdminRobK! I

Thank you, AdminRobK! I understand what you're saying, but I do not know what I should change to avoid the gradients being NA at the start values. Could you point me in the right direction please?

AdminRobK's picture
Joined: 01/24/2014 - 12:15
start values

I don't have any specific advice, beyond adjusting the start values. It looks like the set of start values you chose initialized optimization near the boundary of the feasible space, so that numerically differentiating the fitfunction led to non-finite values of some of the gradient elements.

As far as finding good start values is concerned, mhunter has already in this thread mentioned mxTryHard() and mxAutoStart().

mhunter's picture
Joined: 07/31/2009 - 15:26
Starting values matter

Ahh! I missed that you already had the variances for RIx and RIy. The problem seems to have been bad starting values.

As a general rule, starting values of zero are not a good idea for many parameters in many models. The initial problem was that you started a covariance matrix at all zeros. OpenMx probably (depending on mxOptions) moved these starting values to all 0.1. However, a covariance matrix of all 0.1 is also nonsense. So, you gave OpenMx nonsense implausible starting values and OpenMx did not work.

When I changed the starting value of the variances of RIx and RIy to 1.5, OpenMx did not find a reasonable solution, but at least it was able to start. mxTryHard() tries several "random" starting values, and at lease one of those worked to get a reasonable solution.

OpenMx requires users to think about plausible starting values for their parameters. lavaan generally does not. One option, if you're really struggling to think of plausible starting values is to use mxAutoStart() which uses unweighted (or optionally weighted) least squares to find starting values on many models.

ss <- mxAutoStart(riclpm)
res <- mxRun(ss)
cjvanlissa's picture
Joined: 04/10/2019 - 12:43
Ahh, of course. I somehow

Ahh, of course. I somehow misremembered that mxRun wrapped mxAutoStart. Thank you!