In OpenMx, as model complexity grows, the higher the chance of not setting starting values that will help the optimizer reach a solution. This is particularly true for longitudinal models. This post aims at defining a set of rules that can help you define good starting points and also present what is current available. Starting values can lead to time-consuming, annoying errors and thinking about them during model creation can save you time.
Here are two initial steps for variances and covariances:
- If you know from previous work the variances for certain latent variables and or error measurements, start with these and with covariances that are half of these values;
- If you don’t know the variances, start with 1 for all variances and 0.5 for covariances.
If your model does not have a longitudinal aspect to it, but still runs
into non-pos def issues, then there are two tools that can help, as
these cases tend to be related to zeroed starting points (or values too
close to zero). In these cases one can use mxAutostart() or the
umx() internal function xmuValues(onlyTouchZeros=T), which should be
enough. Both requires data attached to the mx object.
Let us start with a simple misspecification case, which will help us
introduce the use of mxGetExpected to identify problems with the
model.
library(umx)
options(mxByrow = T)
manifests <- c("Y2", "X2", "Y3", "X3", "Y4", "X4", "X1", "Y1", "PGSx", "PGSy")
latents <- c("xir", "yir", "innoX1", "innoX2", "innoX3", "innoX4", "innoY1", "innoY2", "innoY3", "innoY4")
allVariables <- c(manifests, latents)
Al <- c(
NA,NA ,NA ,NA ,NA,NA,NA,"Y1_to_X2","PGSx_to_X2" ,NA,NA,NA ,"innoX1_to_X2",NA,NA ,NA ,NA ,NA ,NA,NA,
NA,NA ,NA ,"Y2_to_X3",NA,NA ,NA ,NA ,"PGSx_to_X3" ,NA,NA,NA ,NA,"innoX2_to_X3" ,NA,NA ,NA ,NA ,NA,NA,
NA,NA ,NA ,NA ,"Y3_to_X4" ,NA ,NA ,NA,"PGSx_to_X4" ,NA ,NA ,NA ,NA,NA,"innoX3_to_X4",NA,NA ,NA,NA ,NA,
NA,NA ,NA ,NA,NA,NA,"X1_to_Y2" ,NA ,NA ,"PGSy_to_Y2",NA,NA ,NA,NA,NA ,NA,"innoY1_to_Y2",NA ,NA ,NA,
"X2_to_Y3" ,NA,NA ,NA ,NA,NA,NA,NA,NA ,"PGSy_to_Y3" ,NA ,NA,NA,NA,NA ,NA,NA,"innoY2_to_Y3",NA,NA,
NA,"X3_to_Y4",NA,NA ,NA,NA,NA,NA,NA ,"PGSy_to_Y4" ,NA ,NA ,NA,NA,NA ,NA,NA ,NA ,"innoY3_to_Y4",NA,
NA,NA,NA,NA ,NA ,NA ,NA,NA ,"PGSx_to_X1" ,NA ,NA ,NA,NA,NA,NA,NA,NA,NA ,NA ,NA,
NA,NA,NA,NA,NA ,NA,NA,NA,NA ,"PGSy_to_Y1" ,NA,NA ,NA,NA,NA ,NA,NA,NA ,NA,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA,
NA,NA,NA,NA ,NA ,NA,NA,NA,NA,NA ,NA ,NA ,NA,NA,NA,NA ,NA,NA,NA ,NA
)
Sl <- c(
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,"PGSx_with_PGSy", NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,"PGSx_with_PGSy",NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,"xir_with_yir",NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, "xir_with_yir",NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
NA, NA,NA,NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA,NA,NA
)
ldX <- mxModel("long2", type = "RAM", manifestVars = manifests, latentVars = latents,
umxMatrix(
name = "A", type = "Full", nrow = length(allVariables), ncol = length(allVariables),
values = c(
# X2 X3 X4 Y2 Y3 Y4 X1 Y1 PGSx PGSy xir yir innoX1 innoX2 innoX3 innoX4 innoY1 innoY2 innoY3 innoY4
0,0, 0, 0, 0, 0, 0, 0,0,0,1,0, 0, 1, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,1,0, 0, 0, 1, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,1,0, 0, 0, 0, 1, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,1, 0, 0, 0, 0, 0, 1, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,1, 0, 0, 0, 0, 0, 0, 1, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,1, 0, 0, 0, 0, 0, 0, 0, 1,
0,0, 0, 0, 0, 0, 0, 0,0,0,1,0, 1, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,1, 0, 0, 0, 0, 1, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0,
0,0, 0, 0, 0, 0, 0, 0,0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0
), labels = Al, free = (!is.na(Al))
),
umxMatrix(name = "S", type = "Full", nrow = length(allVariables), ncol = length(allVariables), values = c(
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1
), labels = Sl, free = (!is.na(Sl))),
mxMatrix(name = "M", "Full", nrow = 1, ncol = length(allVariables), values = c(
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
), free = c(
T, T, T, T,T, T, T, T, F, F, F, F, F, F, F, F, F, F, F, F
), labels = c(
"one_to_X2", "one_to_X3", "one_to_X4", "one_to_Y2", "one_to_Y3", "one_to_Y4", "one_to_X1", "one_to_Y1", "one_to_PGSx", "one_to_PGSy", "one_to_xir", "one_to_yir", "one_to_innoX1", "one_to_innoX2",
"one_to_innoX3", "one_to_innoX4", "one_to_innoY1", "one_to_innoY2", "one_to_innoY3", "one_to_innoY4"
), ),
mxExpectationRAM("A", "S", "F", "M", dimnames = c(manifests, latents)),
mxFitFunctionML())
mldX <- mxGenerateData(ldX, nrows = 1000, returnModel = T, empirical = T)
result <- mxTryHard(mldX) ## Beginning initial fit attemptBeginning fit attempt 1 of at maximum 10 extra triesBeginning fit attempt 2 of at maximum 10 extra triesBeginning fit attempt 3 of at maximum 10 extra triesBeginning fit attempt 4 of at maximum 10 extra triesBeginning fit attempt 5 of at maximum 10 extra triesBeginning fit attempt 6 of at maximum 10 extra triesBeginning fit attempt 7 of at maximum 10 extra triesBeginning fit attempt 8 of at maximum 10 extra triesBeginning fit attempt 9 of at maximum 10 extra triesBeginning fit attempt 10 of at maximum 10 extra tries
##
## All fit attempts resulted in errors - check starting values or model specification
Which results in fit is not finite (The continuous part of the model implied covariance (loc2) is not positive definite. Can you spot the
problem with starting values?
The error is in not setting covariances to the observed PRSs. How did we pick this up? Simply by running:
mxGetExpected(ldX,"covariance")## Y2 X2 Y3 X3 Y4 X4 X1 Y1 PGSx PGSy
## Y2 1 0 0 0 0 0 0 0 0 0
## X2 0 1 0 0 0 0 0 0 0 0
## Y3 0 0 1 0 0 0 0 0 0 0
## X3 0 0 0 1 0 0 0 0 0 0
## Y4 0 0 0 0 1 0 0 0 0 0
## X4 0 0 0 0 0 1 0 0 0 0
## X1 0 0 0 0 0 0 1 0 0 0
## Y1 0 0 0 0 0 0 0 1 0 0
## PGSx 0 0 0 0 0 0 0 0 0 0
## PGSy 0 0 0 0 0 0 0 0 0 0
Or visually:
mxGetExpected(ldX,"covariance") |> cov2cor() |> heatmap(Rowv = NA, Colv = NA)## Warning in cov2cor(mxGetExpected(ldX, "covariance")): diag(.) had 0 or NA
## entries; non-finite result is doubtful
Which results in:
> mxGetExpected(ldX,"covariance")
Y2 X2 Y3 X3 Y4 X4 X1 Y1 PGSx PGSy
Y2 1 0 0 0 0 0 0 0 0 0
X2 0 1 0 0 0 0 0 0 0 0
Y3 0 0 1 0 0 0 0 0 0 0
X3 0 0 0 1 0 0 0 0 0 0
Y4 0 0 0 0 1 0 0 0 0 0
X4 0 0 0 0 0 1 0 0 0 0
X1 0 0 0 0 0 0 1 0 0 0
Y1 0 0 0 0 0 0 0 1 0 0
PGSx 0 0 0 0 0 0 0 0 0 0
PGSy 0 0 0 0 0 0 0 0 0 0
This issue can then be fixed by altering the values in the correct places (note the value, 0.5):
ldX$S$values[9,9]<-0.5; ldX$S$values[10,10]<-0.5
mxGetExpected(ldX,"covariance")## Y2 X2 Y3 X3 Y4 X4 X1 Y1 PGSx PGSy
## Y2 1 0 0 0 0 0 0 0 0.0 0.0
## X2 0 1 0 0 0 0 0 0 0.0 0.0
## Y3 0 0 1 0 0 0 0 0 0.0 0.0
## X3 0 0 0 1 0 0 0 0 0.0 0.0
## Y4 0 0 0 0 1 0 0 0 0.0 0.0
## X4 0 0 0 0 0 1 0 0 0.0 0.0
## X1 0 0 0 0 0 0 1 0 0.0 0.0
## Y1 0 0 0 0 0 0 0 1 0.0 0.0
## PGSx 0 0 0 0 0 0 0 0 0.5 0.0
## PGSy 0 0 0 0 0 0 0 0 0.0 0.5
mldX <- mxGenerateData(ldX, nrows = 1000, returnModel = T, empirical = T)
out <- mxRun(mldX)## Running long2 with 30 parameters
These edits solves the issue and hopefully this section helped show how mxGetExpected can be used to debug the model.
In longitudinal models where at least two phenotypes or their variances are allowed to correlate one can run into non-positive definite matrices. The general idea to avoid this is to set correlation paths starting values that are smaller than the variances of the observed/latent variables. Consider the rather large model below, where X influences causally Y over time points (t_one-1, t_one, and temp). However, observe that the effect is not being carried over by the observed variables, but by latent variables that incorporates the variances coming from the observed. It is perhaps beyond the point, but this model was set up this way so that the individual variances were separated from the variances from the random intercepts xir and yir.
Now the latent variables p and q, have themselves variances associated to them: u and v. U and v are allowed to correlate within wave. Check the model below:
library(umx)
temp = 3; t_one = 2; px = 0.1; py = 0.1; p2q= 0.2; q2p=0.2; immedx = 0.1; immedy = 0.1; arp = 0.7; arq = 0.7; summary = T; id = F; plot = F
model <- mxModel("Hamaker", type = "RAM",
manifestVars = c("PSx","PSy",paste0("X",temp), paste0("X",t_one-1),paste0("X", t_one),
paste0("Y",temp), paste0("Y",t_one-1), paste0("Y",t_one)),
latentVars = c("xir", "yir",paste0("p", temp), paste0("p", t_one-1),paste0("p", t_one),
paste0("q",temp), paste0("q", t_one-1), paste0("q",t_one),
paste0("u",temp), paste0("u",t_one),
paste0("v",temp), paste0("v",t_one)) ,
# Immediate
umxPath(paste0("X",temp), to = paste0("Y",temp), value = immedx) ,
umxPath(paste0("Y",temp), to = paste0("X",temp), value = immedy) ,
umxPath(paste0("X",t_one), to = paste0("Y",t_one), value = immedx) ,
umxPath(paste0("Y",t_one), to = paste0("X",t_one), value = immedy) ,
umxPath(paste0("X",t_one-1), to = paste0("Y",t_one-1), value = immedx) ,
umxPath(paste0("Y",t_one-1), to = paste0("X",t_one-1), value = immedy) ,
# Instruments
umxPath("PSx", to = c(paste0("X",temp),paste0("X",t_one),paste0("X",t_one-1)), value=px),
umxPath("PSy", to = c(paste0("Y",temp),paste0("Y",t_one),paste0("Y",t_one-1)), value=py),
# RIs
umxPath("xir", to = c(paste0("X",temp),paste0("X",t_one),paste0("X",t_one-1)), fixedAt=1),
umxPath("yir", to = c(paste0("Y",temp),paste0("Y",t_one),paste0("Y",t_one-1)), fixedAt=1),
# Causal
umxPath(paste0("p",t_one-1), to = paste0("q",t_one), value=p2q),
umxPath(paste0("q",t_one-1), to = paste0("p",t_one), value=q2p),
umxPath(paste0("p",t_one), to = paste0("q",temp), value=p2q),
umxPath(paste0("q",t_one), to = paste0("p",temp), value=q2p),
# AR
umxPath(paste0("p",t_one-1), to = paste0("p",t_one), value=arp),
umxPath(paste0("q",t_one-1), to = paste0("q",t_one), value=arq),
umxPath(paste0("p",t_one), to = paste0("p",temp), value=arp),
umxPath(paste0("q",t_one), to = paste0("q",temp), value=arq),
# there is no u or v in this first wave
umxPath(var = c(paste0("p",t_one-1), paste0("q",t_one-1)), value=0.2),
umxPath(paste0("p",t_one-1), with = paste0("q",t_one-1)) ,
# p,q,u,v
umxPath(paste0("p",temp), to = paste0("X",temp), fixedAt=1),
umxPath(paste0("p",t_one), to = paste0("X",t_one), fixedAt=1),
umxPath(paste0("p",t_one-1), to = paste0("X",t_one-1), fixedAt=1),
umxPath(paste0("q",temp), to = paste0("Y",temp), fixedAt=1),
umxPath(paste0("q",t_one), to = paste0("Y",t_one), fixedAt=1),
umxPath(paste0("q",t_one-1), to = paste0("Y",t_one-1), fixedAt=1),
umxPath(paste0("u",temp), to = paste0("p",temp), fixedAt=1),
umxPath(paste0("v",temp), to = paste0("q",temp), fixedAt=1),
umxPath(paste0("u",temp), with = paste0("v",temp), value=0.2),
umxPath(paste0("u",t_one), to = paste0("p",t_one), fixedAt=1),
umxPath(paste0("v",t_one), to = paste0("q",t_one), fixedAt=1),
umxPath(paste0("u",t_one), with = paste0("v",t_one), value=0.2),
umxPath(var = c(paste0("u",t_one), paste0("v",t_one)), value=0.2),
umxPath(var = c(paste0("u",temp), paste0("v",temp)), value=0.2),
umxPath(var = c("xir", "yir")),
umxPath(means = c("PSx","PSy",paste0("X",temp), paste0("X", t_one-1),paste0("X", t_one),
paste0("Y",temp), paste0("Y",t_one-1),paste0("Y",t_one))),
umxPath("PSx", with ="PSx", value=0.2),
umxPath("PSy", with ="PSy", value=0.2),
umxPath("PSx", with = "PSy") ,
umxPath("xir", with = "yir")
)
# mxCheckIdentification(model)$status # TRUE
# plot(model, means = F,spline=F) # If you want to visualize the model
run <- mxGenerateData(model, nrows = 1000, empirical = T, returnModel =T)
out <- mxRun(run)The above exits with the error that the model implied covariance is not positive definite. The steps for debugging the error includes:
- Extracting the expected covariance matrix with
badExpCov <- mxGetExpected(model, "covariance")> badExpCov
PSx PSy X3 X1 X2 Y3 Y1 Y2
PSx 0.200000000 0.000000000 0.020202020 0.020202020 0.020202020 0.002020202 0.002020202 0.002020202
PSy 0.000000000 0.200000000 0.002020202 0.002020202 0.002020202 0.020202020 0.020202020 0.020202020
X3 0.020202020 0.002020202 0.535139884 0.122722171 0.324613815 0.523156413 0.079746965 0.302299765
X1 0.020202020 0.002020202 0.122722171 0.208162432 0.154494439 0.079746965 0.041220284 0.070196919
X2 0.020202020 0.002020202 0.324613815 0.154494439 0.369635751 0.302299765 0.070196919 0.326660545
Y3 0.002020202 0.020202020 0.523156413 0.079746965 0.302299765 0.535139884 0.122722171 0.324613815
Y1 0.002020202 0.020202020 0.079746965 0.041220284 0.070196919 0.122722171 0.208162432 0.154494439
Y2 0.002020202 0.020202020 0.302299765 0.070196919 0.326660545 0.324613815 0.154494439 0.369635751
- Check the eigenvalues of it if nothing obvious (like zero or negative values on diagonal)
eigen(badExpCov)
eigen() decomposition
$values
[1] 1.596298e+00 2.664973e-01 2.404901e-01 1.987381e-01 1.814107e-01 1.424416e-01 1.763735e-17
[8] -1.290923e-16
$vectors
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] -0.01783117 -0.1248507 -0.4071083 0.68457821 -0.57578040 0.1242861 -0.05149735 0.00936744
[2,] -0.01783117 -0.1248507 0.4071083 0.68457821 0.57578040 0.1242861 0.05149735 -0.00936744
[3,] -0.54765994 0.4157610 -0.1506645 0.08700418 0.05418694 -0.1401487 0.63411315 0.26882059
[4,] -0.15269004 -0.3968405 -0.4916283 0.02611207 0.37377878 -0.5643771 -0.24941428 0.23745213
[5,] -0.42004333 -0.3925210 -0.2643191 -0.15199047 0.16071755 0.3826092 0.18177199 -0.60931456
[6,] -0.54765994 0.4157610 0.1506645 0.08700418 -0.05418694 -0.1401487 -0.63411315 -0.26882059
[7,] -0.15269004 -0.3968405 0.4916283 0.02611207 -0.37377878 -0.5643771 0.24941428 -0.23745213
[8,] -0.42004333 -0.3925210 0.2643191 -0.15199047 -0.16071755 0.3826092 -0.18177199 0.60931456
Eigenvalues are a special set of numbers that are associated with square
matrices. They are known as characteristic roots, characteristic values,
proper values, latent roots, and principal components. Eigenvalues
are commonly used in factor analysis and genomic population
stratification. The eigenvalue of a square matrix D is det(D -
I) = 0. This equatin will result in a polynomial where the
values of
are the eigenvalues.
Examining the eigenvalues can help identify the paramaters that are
causing problems. In the output above, the last eigenvalue is negative,
so possibly the eigenvectors will suggest which variables are involved.
Looks like biggest numbers in column 8 of eigenvectors are 5th & 8th
rows. Clearly, this model was not written following the recommendations
at the beginning of this text. Variance for v and u are identical to the
covariance between them (0.2). There are two solutions possible, both
works. Either one raises the starting value for the variance or one
reduces the value for the correlation. In general blowing up values for
variances won’t cause great problems, especially in longitudinal models.
My inclination is always to raise the variance to 1. In this case we did
the opposite, we reduced the covariance. These rows in the eigen vectors
are those associated with X2 and Y2, therefore the error should
disappear if you give more variable-specific variance to X2 and Y2. This
can be done by tweaking the script line umxPath(paste0("u",t_one), with = paste0("v",t_one), value=0). Finally, re-run the code.
You can inspect if you edited the correct parts by running:
badExpCov - mxGetExpected(run, 'covariance')And by using mxRun() you can see that the error is gone.