# Ordinal covariance is not positive definite in data

I would like to ask for help with an error:

Details about my openMX run:

OpenMx version: 2.7.12 [GIT v2.7.12-dirty]

R version: R version 3.4.1 (2017-06-30)

Platform: x86_64-w64-mingw32

Default optimiser: CSOLNP

I am trying to run a univariate model with three covariates-one categorical and two continuous. Unfortunately, when I get to the phase of model running (mxRun), no matter what I do, I receive the following error:

"Running ACE_Hetero with 18 parameters

Error: The job for model 'ACE_Hetero' exited abnormally with the error message: fit is not finite (Ordinal covariance is not positive definite in data 'DZF.data' row 112 (loc2))

In addition: Warning message:

In model 'ACE_Hetero' Optimizer returned a non-zero status code 10. Starting values are not feasible. Consider mxTryHard()".

In attempts to resolve the first error, I tried removing observations with missing values, removing one of the covariates from the model, replacing the starting values, playing with the "connect" command when specifying the covariates covariance matrix... and nothing helped. It is always the same error.

I an attaching the syntax. I would really appreciate your input regarding what is wrong.

Thank you!!

#define start values

meanFAC1 <-mean(na.omit(WDtotal$POSFAC_1_Res))

meanFAC2 <-mean(na.omit(WDtotal$POSFAC_2_Res))

StartMean <-mean(c(meanFAC1,meanFAC2))

StartMean.Age <-mean(WDtotal$Age6Y)

StartMean.Relig <-mean(WDtotal$MomReligion, na.rm=T)

StartMean.Sib <- mean(WDtotal$sibNum,na.rm=T)

```
```require("OpenMx")

# Select Variables for Analysis

selVars <- c('POSFAC_1_Res','POSFAC_2_Res','LabHome6', 'MomReligion','sibNum')

aceVars <- c("A1","C1","E1","A2","C2","E2")

# Select Data for Analysis

mzM <- subset(WDtotal, zygosity==1 & sex_a=='sex1',selVars)

mzF <- subset(WDtotal, zygosity==1 & sex_a=='sex2',selVars)

dzM <- subset(WDtotal, zygosity==2 & sex_a=='sex1',selVars)

dzF <- subset(WDtotal, zygosity==2 & sex_a=='sex2',selVars)

# Path objects for Multiple Groups

manifestVars=selVars

latentVars=aceVars

# specify paths for the variances and means of the latent variables

latVariances <- mxPath( from=aceVars, arrows=2,free=FALSE, values=1 )

varCovar <- mxPath( from= c('LabHome6','MomReligion','sibNum'), arrows=2,

free=(c(F,T,T)), values=c(1,1,1),

label=c('placeVar','religion','siblings'),

connect='unique.bivariate' )

# means of latent variables

latMeans <- mxPath( from="one", to=aceVars, arrows=1,free=FALSE, values=0 )

#specify paths for the means of the observed variables. Because the scores are residual scores from sex and age, I constrain the means to be the same for boys and girls.

obsMeansM <- mxPath( from="one", to=selVars, arrows=1,

free=c(T,T,F,T,T), values=c(StartMean, StartMean,0,StartMean.Relig,StartMean.Sib),

labels=c("meanFAC","meanFAC","meanPlace","meanReligion","meanSiblings") )

obsMeansF <- mxPath( from="one", to=selVars, arrows=1,

free=c(T,T,F,T,T), values=c(StartMean, StartMean,0,StartMean.Relig,StartMean.Sib),

labels=c("meanFAC","meanFAC","meanPlace","meanReligion","meanSiblings") )

# path coefficients for twin 1.

pathAceTM1 <- mxPath( from=c("A1","C1","E1", "LabHome6","MomReligion","sibNum"),

to="POSFAC_1_Res", arrows=1, free=c(T,T,T,T,T,T),

values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartMean.Relig,StartMean.Sib),

label=c("aM","cM","eM","place","Religion","Siblings") )

pathAceTF1 <- mxPath( from=c("A1","C1","E1", "LabHome6","MomReligion","sibNum"),

to="POSFAC_1_Res", arrows=1, free=c(T,T,T,T,T,T),

values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartMean.Relig,StartMean.Sib),

label=c("aF","cF","eF","place","Religion","Siblings") )

# path coefficients for twin 2

pathAceTM2 <- mxPath( from=c("A2","C2","E2", "LabHome6","MomReligion","sibNum"),

to="POSFAC_2_Res", arrows=1, free=c(T,T,T,T,T,T),

values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartMean.Relig,StartMean.Sib),

label=c("aM","cM","eM","place","Religion","Siblings") )

pathAceTF2 <- mxPath( from=c("A2","C2","E2", "LabHome6","MomReligion","sibNum"),

to="POSFAC_2_Res", arrows=1, free=c(T,T,T,T,T,T),

values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartMean.Relig,StartMean.Sib),

label=c("aF","cF","eF","place","Religion","Siblings") )

# covariance between C1 & C2

covC1C2 <- mxPath( from="C1", to="C2", arrows=2,free=FALSE, values=1 )

# covariance between A1 & A2 in MZ's

covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,free=FALSE, values=1 )

# covariance between A1 & A2 in DZ's

covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,free=FALSE, values=.5 )

#call the data.frame with the MZ raw data, mzData, and the DZ raw data

dataMZM <- mxData( observed=mzM, type="raw" )

dataMZF <- mxData( observed=mzF, type="raw" )

dataDZM <- mxData( observed=dzM, type="raw" )

dataDZF <- mxData( observed=dzF, type="raw" )

# Combine Groups

#arrange together all the paths that are the same for MZ and DZ

pathsM <- list(latVariances, latMeans, obsMeansM,

pathAceTM1, pathAceTM2, covC1C2, varCovar )

pathsF <- list( latVariances, latMeans, obsMeansF,

pathAceTF1, pathAceTF2, covC1C2, varCovar )

#Build a model seperately for MZ and DZ and for each sex

modelMZM <- mxModel(model="MZM", type="RAM", manifestVars=selVars,

latentVars=aceVars, pathsM, threshold, covA1A2_MZ, dataMZM )

modelDZM <- mxModel(model="DZM", type="RAM", manifestVars=selVars,

latentVars=aceVars, pathsM, threshold, covA1A2_DZ, dataDZM )

modelMZF <- mxModel(model="MZF", type="RAM", manifestVars=selVars,

latentVars=aceVars, pathsF, threshold, covA1A2_MZ, dataMZF )

modelDZF <- mxModel(model="DZF", type="RAM", manifestVars=selVars,

latentVars=aceVars, pathsF, threshold, covA1A2_DZ, dataDZF )

obj <- mxFitFunctionMultigroup(c("MZM","DZM","MZF","DZF"))

modelACE.Hetero <- mxModel(model="ACE_Hetero", modelMZM, modelDZM,modelMZF, modelDZF, obj, CI )

fitACE.Hetero <- mxRun(modelACE.Hetero)

sumACE.Hetero <- summary(fitACE.Hetero)

## some suggestions

`meanFAC1`

etc.? Also, is there anything unusual about row 112 of the female DZ dataset?I have a few generic suggestions. One is to go with the advice in the warning message and use

`mxTryHardOrdinal()`

instead of`mxRun()`

. Another idea is to run your model through`mxAutoStart()`

before`mxRun()`

. Something else you could try is to switch optimizers (e.g.,`mxOption(NULL,"Default optimizer","SLSQP")`

switches to SLSQP), but I doubt that will make a difference.If none of those things help, then I'd be curious to see what happens if you do

plan <- omxDefaultComputePlan()

plan$steps$GD <- mxComputeNelderMead()

and include

`plan`

in your container model,`modelACE.Hetero`

. You'd then be using my Nelder-Mead implementation as the optimizer, and Nelder-Mead can still work if the start values are infeasible (as long as at least one vertex of the initial simplex is feasible).Other possibilities would be to use weighted least-squares instead of maximum-likelihood, or to pass a non-default value for argument

`jointConditionOn`

to`mxFitFunctionML()`

, but I'm not sure offhand how you would do either of those things with MxModels that have`type="RAM"`

.Log in or register to post comments

In reply to some suggestions by AdminRobK

## thanks and replies

Thank you for all the offers.

1)Regarding the start values-

The start values of the means are the mean of each variable, which are: 4.028427e-17 for the observed factor ("StartMean"), 0 for the categorical covariate, and 3.015 (StartMean.Relig) and 1.78 (StartMean.Sib) for the two continuous covariates.

The start values of the ACE variances ("StartVar.ACE") are the variance of the observed variable divided by 3, which is .52. the start values of the covariates are 1 for the categorical covariate and the variance of the continuous covariates-approximately 1 and 2.5.

2) Regarding the start values: I don't think there is anything unusual in row 112. What I see is:

POSFAC_1_Res POSFAC_2_Res LabHome6 MomReligion sibNum

305 -1.970493 0.4489254 LabHome61 2 3

It is very similar to what I see in other rows of the data... Moreover, following your question I tried to remove this row from the data, and it just told me that there is a problem in another row (I did this procedure a few more times, each time I got problem with a new row).

3) I tried your suggestions, including adding the Nelder-Med implementation. I got, aside from the previous error, the following message:

"The job for model 'ACE_Hetero' exited abnormally with the error message: 1:initial simplex is not feasible; specify it differently, try different start values, or use mxTryHard()" ......"In model 'ACE_Hetero' Optimizer returned a non-zero status code 3. The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution. "

4) Now I tried to run the syntax only with one covariate (each time with a different covariate) and the syntax worked! However, whenever I ran it with more than one covariate, the error message appeared again.

Maybe the problem is in the relation between the covariates?

Thank you so much

Log in or register to post comments

In reply to thanks and replies by lior abramson

## more ideas

mxGetExpected(modelACE.Hetero, "covariance", "MZM")

. Check to see if there is any obvious sign of non-positive-definiteness, such as a non-positive variance, or a covariance equal to one of its corresponding variances. You can further use

`eigen()`

to check if the model-expected covariance matrices are positive-definite at the start values, which will be the case for each such matrix if all of its eigenvalues are positive.However, this,

, makes me suspect that one of the covariates is perfectly correlated with another covariate within one or more of your groups, or that, in one or more of your groups, the joint covariance matrix of the three covariates is non-positive-definite (or even that the joint covariance matrix of phenotypes and two of the three covariates is non-positive-definite). An approximate, quick-and-dirty check of the latter might be to do a simple OLS regression of one covariate onto the other two, and see if the result is pathological in some way, or to see whether the covariates' Spearman correlation matrix is positive-definite.

BTW, it might simplify troubleshooting if you treat your ordered-categorical covariate as though it were continuous for the time being. In particular, I'd be curious to see if that changes the error messages you see.

Edit: I could be missing something, but I don't see that you set the variances of the covariates to be free and positive anywhere in the script. I think you need to redefine

`varCovar`

differently. If that's the issue, then it would explain why nothing I suggested so far has helped.Double edit: Yes, you could be encountering a known issue.

Log in or register to post comments

In reply to more ideas by AdminRobK

## further checks and hopefully final answer

1) The model actually works when I use two out of three covariates (and any combination, i.e., 1-2,1-3,2-3, works OK). In these cases, the expected covariance matrix looks, for example, like this:

$DZF

POSFAC_1_Res POSFAC_2_Res MomReligion sibNum

POSFAC_1_Res 8.051072 7.650536 1 2.5

POSFAC_2_Res 7.650536 8.051072 1 2.5

MomReligion 1.000000 1.000000 1 0.0

sibNum 2.500000 2.500000 0 1.0

However, when I enter all the three covariates the expected covariate matrix looks like this:

$DZF

POSFAC_1_Res POSFAC_2_Res LabHome6 MomReligion sibNum

POSFAC_1_Res 12.80107 12.40054 3.5 3.5 2

POSFAC_2_Res 12.40054 12.80107 3.5 3.5 2

LabHome6 3.50 3.50 0.0 1.0 1

MomReligion 3.50 3.50 1.0 0.0 1

sibNum 2.00 2.00 1.0 1.0 0

I am guessing the problem is that the covariance of each covariate with itself is somehow becoming zero. Is that an "obvious sign of non-positive-definiteness"?

two of the covariates are highly correlated with one another (r=.66) but they are not perfectly correlated. The OlS regression and spearman correlation does not indicate anything weird, I think.

In light of this, I am guessing that the only solution would be to drop one of the two covariates which are highly correlated from the model. Would you agree that this is the best solution?

3) regarding the definition of the covariates' variances: I was sure that by this command:

"varCovar <- mxPath( from= c('LabHome6','MomReligion','sibNum'), arrows=2,

free=(c(T,T,T)), values=c(1,1,1),

label=c('placeVar','religion','siblings'),

connect='unique.bivariate' )"

I am setting the variances to be free ("free=c(T,T,T)") and positive ("values=c(1,1,1). Am I missing something regarding the definitions?

Thank you for your time and patience

Log in or register to post comments

In reply to further checks and hopefully final answer by lior abramson

## and one more thing

Log in or register to post comments

In reply to further checks and hopefully final answer by lior abramson

## 'connect' argument

Per the documentation for

`mxPath()`

,`connect='unique.bivariate'`

only creates paths between the three variables (i.e., covariance paths). Unless I am overlooking it, I don't see that your script creates any double-headed paths going from each covariate toitself, to create paths for the variances of the covariates. I think you need to use`connect='unique.pairs'`

instead.Log in or register to post comments

In reply to 'connect' argument by AdminRobK

## It worked!

After changing this I could run the mxAutoStart command, and by using mxTryHardOrdinal the model was able to run. Thanks!

Log in or register to post comments

In reply to It worked! by lior abramson

## OK! Glad to hear it's

Log in or register to post comments