You are here

Maximum Likelihood for Cross-Lagged Panel Models with Fixed Effects - Allison et al. 2017

6 posts / 0 new
Last post
ppl's picture
ppl
Offline
Joined: 01/12/2023 - 09:19
Maximum Likelihood for Cross-Lagged Panel Models with Fixed Effects - Allison et al. 2017

Hey all,

I just started recently to use OpenMx for academic reasons. I'm currently trying to reproduce the empirical example of the paper of Allison et al. 2017 in which they use a Maximum Likelihood Structural Equation Model (ML-SEM) for dynamic panel models.
The authors provide the code of the following software packages: R-lavaan, Stata, Mplus, SAS-Proc Calis. However, I'm using OpenMx in R since it seems to be more flexible for changes and adjustments than lavaan.

You can find my code and the paper in the attachment. Unfortunately I get the following error message, but I can't find my error:
Error: The job for model 'Test1' exited abnormally with the error message: fit is not finite (The continuous part of the model implied covariance (loc2) is not positive definite in data 'Test1.data' row 499.

However, even if I delete row 499, the error remains. I hope that someone here in the forum can help me. Thank you very much in advance!

AdminHunter's picture
Offline
Joined: 03/01/2013 - 11:03
Check the data and starting values

Welcome to OpenMx!

I don't see anything obviously wrong with your model specification, but there are a couple of standard things to check.

At the risk of being a bit pedantic, let me explain the error message. "fit is not finite": I tried to calculate the fit function value, but I got was infinity, negative infinity, or NA. "The continuous part of the model implied covariance is not positive definite": The model-implied covariance matrix is not acting like an actual covariance matrix. The "positive-definite" language is a generalization of the requirement that variances are always greater than zero (i.e., positive): by analogy, covariance matrices are always positive definite, but your model-implied covariance matrix appears not to be. Now, on to what you can do.

First, it looks like you're using the haven package to read in data. My experiences with haven have not been positive. Check to make sure wages is a data.frame or matrix by running is(wages). If not, you can use read.table() instead of haven, or cast it to a data.frame with ds <- as.data.frame(wages).

Second, you can check the model-implied covariance matrix at the starting values with

mxGetExpected(model1, 'covariance')
# 
round(eigen(mxGetExpected(model1, 'covariance'))$values, 3)
# [1] 473.035   1.172   0.688   0.578   0.551   0.285   0.000   0.000   0.000   0.000   0.000   0.000
#[13]   0.000   0.000   0.000   0.000   0.000  -0.138  -0.315  -1.018  -2.912  -4.927

Your starting values imply 11 zero eigenvalues and 5 negative ones. This is probably the issue. Set starting values (particularly variances) large enough that all the eigenvalues of the model-implied covariance matrix are positive. A built-in helper function for this is mxAutoStart().

smodel <- mxAutoStart(model1)
run <- mxRun(smodel)

Third, it looks like you have a free parameter named "ed" and also a variable named "ed". I don't know that this causes a problem or your problem, but it's probably not a good idea. Use a different variable name or free parameter label.

Fourth and finally, think about ways to make the model specification shorter and less tedious. This is just a style question. For example, many of the repeated paths in your script could be replaced with the following:

mxPath('ed', paste0('wks', 2:7), values=1, labels='edpath')
mxPath(from=paste0(c('wks', 'lwages', 'union'), rep(1:6, each=3)),
       to=rep(paste0('wks', 2:7), each=3),
       values=1,
       labels=rep(c('lambda', 'beta2', 'beta3'), times=6))

Cheers!

ppl's picture
ppl
Offline
Joined: 01/12/2023 - 09:19
Dear Micheal Hunter,

Dear Micheal Hunter,

thank you very much for your detailed answer. I'll check whether I can fix my problems with the help of your comments in the following days. Thanks a lot!

All the best to you!

ppl's picture
ppl
Offline
Joined: 01/12/2023 - 09:19
Any further hints? - lavaan, umx packages

Thanks again for your comment. I tried to solve the problem with the help of your instructions but unfortunately it didn't work. An updated version can be found attached ("TestOpenMx_new.R"). Among other things, I updated the starting parameters with regard to the suggestions produced by "mxAutoStart()" but, as I said, it still did not solve the problem. In order to approach to the core of the issue, I tried to estimate a short or reduced-version of my general model. But even that shows the same error message ("TestOpenMx_shortmodel").

Additionally, I converted the correct lavaan-Code, provided by the authors, by using the "metaSEM"-package into OpenMx. It produced some results but those results were completely different than the original results provided by lavaan. Do you have any ideas what could be the reason for this? I also attached this R-code ("lavaan.R").

Furthermore, I stumbled across the umx package. However, I can't find any explanations where exactly the difference to OpenMx is. Out of curiosity I tried the OpenMx script slightly modified with the umx package. However, this has also led to no positive result ("umx_new.R" - However, this might not be important right now). Maybe someone of you has experience with this and knows the difference to OpenMx?

I hope my explanations are not too excessive. Basically, I am trying to reproduce the paper pinned above with OpenMx. Since this doesn't work, I'm currently trying different possible solutions, but they don't currently lead to the right solution.

Hopefully, someone in this forum is able to help me to solve the issue. Many thanks in advance!

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Missing the residual variances on 17/23 variables

I plotted your model to take a look at it. I also had a look at the eigenvalues of the expected covariance matrix, like this:

 eigen(mxEval(solve(diag(23)-A)%&%S,model1))

which revealed a lot of negative eigenvalues. The diagram shows that there aren't residual variances on many of the variables. Changing starting values won't fix this problem.
require(umx); plot(model1)
I partly fixed it with making all variables have residual variances by replacing the few variables that did have residuals with the entire list of them. Then using mxTryHard() gets over the issue that these were all starting at zero. However, better is to explicitly tell it to start the residual variances at roughly the observed variances of the manifest variables, and to avoid giving the variables high covariances between them.

library(OpenMx)
mxOption(key='Number of Threads', value=parallel::detectCores() - 1)
 
# arrow = 1  from -> to
# arrow = 2  from <-> to
# free = FALSE -> fixed path (parameter values are ex ante fix)
# free = TRUE is a free path (parameter values are estimated)
 
# from: sources of the new paths
# values: Starting values of parameters
# connect: type of source to sink connection
# labels: name of the paths
 
# values is the starting value of a path # reasonable values for faster convergence
 
# first try to replicate model in
# "Linear dynamic panel-data estimation using maximum likelihood and structural equation modeling"
 
# yit = λyit−1 + x∗itβ + w∗iδ + αi + ξt + υit ∀i, t
# E(υit|yt−1,i , xt,i, wi, αi) = 0 ∀i, 
library(haven)
wages <- haven::read_dta("https://www3.nd.edu/~rwilliam/statafiles/wageswide.dta")
#is(wages)
 
wages[,c(2:24)]
wages <- as.data.frame(scale(wages))
#wages <- as.matrix(wages)
#is(wages)
#View(wages)
 
head(wages)
summary(wages)
## 595 households who reported a non-zero wage (measured for 7 years: 1976 - 1982)
# y = wks: number of weeks employed in each year 
# x = union (dummy): 1 if wage set by union contract
# w = lwage: ln of wage in each year
# z = ED: years of education in 1976
 
data <- mxData(wages, type = "raw")
 
#mxTypes()
 
latents <- "alpha"
 
manifests <- c("wks1", "union1", "lwage1", "wks2", "union2", "lwage2",
               "wks3", "union3", "lwage3", "wks4", "union4", "lwage4",
               "wks5", "union5", "lwage5", "wks6", "union6", "lwage6",
               "wks7", "union7", "lwage7", "ed")
 
cov_data <- cov(wages[,manifests])
 
#mxTypes()
model1 <- mxModel("Test1", type = "RAM",
                  # Manifest variables
                  manifestVars = c("wks1", "union1", "lwage1", "wks2", "union2", "lwage2",
                                   "wks3", "union3", "lwage3", "wks4", "union4", "lwage4",
                                   "wks5", "union5", "lwage5", "wks6", "union6", "lwage6",
                                   "wks7", "union7", "lwage7", "ed"),
                  # Latent variable
                  latentVars = c("alpha"),
                  #mxPath(from = "alpha", to = manifests,
                  #       free = FALSE, values = c(rep(1, 21), 0)),
                  # Correlation between alpha and all exogenous variables, except "ed"
                  mxPath(from = c("alpha"),
                         to =  c("wks1", "union1", "lwage1", "union2", "lwage2",
                                 "union3", "lwage3", "union4", "lwage4",
                                 "union5", "lwage5", "union6", "lwage6",
                                 "union7", "lwage7", "ed"),
                         arrows = 2, connect = "unique.pairs",
                         free = FALSE, values = c(rep(1, 15), 0)),
                  # Correlation between alpha and endogenous variables, fixed effects
                  mxPath(from = c("alpha"),
                         to = c("wks2", "wks3", "wks4", "wks5", "wks6", "wks7"),
                         arrows = 1, connect = "unique.pairs",
                         free = FALSE, values = c(rep(1, 6))),
                  # Correlation between exogenous and endogenous variables
                  mxPath(from = c("ed", "wks1", "lwage1", "union1"),
                         to = c("wks2"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("ed", "wks2", "lwage2", "union2"),
                         to = c("wks3"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("ed", "wks3", "lwage3", "union3"),
                         to = c("wks4"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("ed", "wks4", "lwage4", "union4"),
                         to = c("wks5"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("ed", "wks5", "lwage5", "union5"),
                         to = c("wks6"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("ed", "wks6", "lwage6", "union6"),
                         to = c("wks7"),
                         arrows = 1, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4)),
                         labels = c("edpath", "lambda", "beta2", "beta3")),
                  mxPath(from = c("wks2"),
                         to = c("union3", "union4", "union5", "union6"),
                         arrows = 2, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 4))),
                  mxPath(from = c("wks3"),
                         to = c("union4", "union5", "union6"),
                         arrows = 2, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 3))),
                  mxPath(from = c("wks4"),
                         to = c("union5", "union6"),
                         arrows = 2, connect = "unique.pairs",
                         free = TRUE, values = c(rep(1, 2))),
                  mxPath(from = c("wks5"),
                         to = c("union6"),
                         arrows = 2, connect = "unique.pairs",
                         free = TRUE, values = 1),
                  # Residual variances
                  mxPath(from = manifests,
                         arrows = 2, connect = "single",
                         free = TRUE, values = c(rep(10, 5)),
                         labels = rep("residual", 6)),
                  mxPath(from = "one", to = manifests,
                         arrows = 1, free = TRUE, values = 1),
                  #Data
                  mxData(wages, type = "raw", numObs = 595))
model1$matrices$A
round(eigen(mxGetExpected(model1, "covariance"))$values, 3)
smodel1 <- mxAutoStart(model1)
mxGetExpected(smodel1, "covariance")
fitmodel <- mxTryHard(model1)
summary(fitmodel)

HTH!

File attachments: 
ppl's picture
ppl
Offline
Joined: 01/12/2023 - 09:19
Dear Michael Neale,

Dear Michael Neale,

thank you very much for your helpful answer. It helped a lot to get rid of the error messages!

All the best to you!