Maximum Likelihood for Cross-Lagged Panel Models with Fixed Effects - Allison et al. 2017
ppl
Joined: 01/12/2023
| Attachment | Size |
|---|---|
| OpenMX_code | 6.91 KB |
| Allison et al. 2017 - ML for Cross-Lagged Panel Models with Fixed Effects | 776.15 KB |
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!
Check the data and starting values
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
havenpackage to read in data. My experiences withhavenhave not been positive. Check to make surewagesis adata.frameormatrixby runningis(wages). If not, you can useread.table()instead ofhaven, or cast it to adata.framewithds <- 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!
Log in or register to post comments
In reply to Check the data and starting values by AdminHunter
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!
Log in or register to post comments
Any further hints? - lavaan, umx packages
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!
Log in or register to post comments
In reply to Any further hints? - lavaan, umx packages by ppl
Missing the residual variances on 17/23 variables
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!
Log in or register to post comments
In reply to Missing the residual variances on 17/23 variables by AdminNeale
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!
Log in or register to post comments