Correcting non-positive definite matrices and misspecification related to starting values

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.

Misspecification

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.

Non-positive definite matrices

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 - $\lambda$ I) = 0. This equatin will result in a polynomial where the values of $\lambda$ 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.

lf-araujo/non-pos def.md

Correcting non-positive definite matrices and misspecification related to starting values

Misspecification

Non-positive definite matrices