Login

News

Ordinal covariance is not positive definite in data

Posted on
No user picture. lior abramson Joined: 07/21/2017
Forums
OpenMx Error Messages
Hi,
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)

Replied on Wed, 09/27/2017 - 11:46
Picture of user. AdminRobK Joined: 01/24/2014

I don't see any obvious problem with the script, but it's hard to troubleshoot without the data or even the start values you're using. Could you provide the numeric values of 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
Replied on Fri, 09/29/2017 - 05:53
No user picture. lior abramson Joined: 07/21/2017

Hi,
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
Replied on Fri, 09/29/2017 - 11:36
Picture of user. AdminRobK Joined: 01/24/2014

I guess the start values could somehow still be really bad, since Nelder-Mead couldn't even get off the ground. You can get the model-expected covariance matrices at the start values for each of your groups via (for example)

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,

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?

, 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
Replied on Sun, 10/01/2017 - 21:57
No user picture. lior abramson Joined: 07/21/2017

So, I used the mxGetExpected and this is what I found out:

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
Replied on Mon, 10/02/2017 - 10:23
Picture of user. AdminRobK Joined: 01/24/2014

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?

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 to itself, 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
  Login

News