Expected covariance matrix is not positive-definite in data row... at iteration...

22 replies [Last post]
dadrivr's picture
Offline
Joined: 01/19/2010

Hi guys,

I'm trying to fit a saturated model where the variable, 'manifests', includes all of the variables in the model. The non-saturated structural model runs fine, but I get an error when I fit the saturated model:

Error: The job for model 'Saturated Model2' exited abnormally with the error message: Expected covariance matrix is not positive-definite in data row 79 at iteration 0.

I have tried fiddling with the starting values, but to no avail. By the way, I know that all=TRUE is being deprecated, and I have written the same code in a longer form without all=TRUE, but this was shorter, so I thought it'd be easier to post. Below is the model:

saturated_model <- mxModel(paste("Saturated Model"),
type="RAM",
manifestVars=manifests,

#Variances and Covariances
mxPath(
from=manifests,
arrows=2,
free=TRUE,
all=TRUE,
values=.5,

#Means
mxPath(
from="one",
to=manifests,
arrows=1,
free=TRUE,

mxData(observed=modeldata, type="raw")
)
saturatedModel_run <- mxRun(saturated_model)

Any help would be greatly appreciated. Thanks!

Ryne's picture
Offline
Joined: 07/31/2009
So this is a starting values

So this is a starting values problem, and one caused by your use of the "all" argument. You're creating all of your variances and covariances with the same function call, so whatever starting value you assign to the variances gets assigned to the covariances as well. That means that you've assigned all variables to be perfectly correlated.

I'll also point out that the "connect" argument is much more flexible than the "all" argument. connect="all.pairs" does what "all=TRUE" does, creating all possible paths for a set of variables. For two-headed arrows, that's actually a pain, as it creates a with b AND b with a. You have four options for connect, which will serve your needs well.

all.pairs does what all=TRUE does.
unique.pairs will do all variables with all variables EXCEPT repeat the redundant paths. a with b will be created, but b with a won't.
all.bivariate and unique.bivariate do the same thing, but leave out the variance terms (a with a is dropped).

So what do you do? Use two mxPath statements. In the first, leave out the all and connect arguments and create the variance terms. In the second, use connect=unique.bivariate to make all of the covariances. Make sure the start values for the covariances are lower than the variances and you're all set.

yenni's picture
Offline
Joined: 10/08/2013
Expected covariance matrix is not positive-definite

Hi guys,

I'm trying to fit a saturated model. The non-saturated structural model runs fine, but I get an error when I fit the saturated model: Error: The job for model 'Saturated Model' exited abnormally with the error message: Expected covariance matrix is not positive-definite in data row 1041 at iteration 0(minor iteration 1).

I have tried fiddling with the starting values and did several suggestions from this forum, but it still failed,so I thought it'd be easier to post.
I post you the data set and my saturated file. I hope someone can help me. Thanks.

AttachmentSize
saturatedmodelYA_1.txt 55.88 KB
han7com.dat 328.83 KB
Ryne's picture
Offline
Joined: 07/31/2009
This isn't a saturated model.

This isn't a saturated model. This is a common factor model with no residual variance terms. As all 50-something manifest variables are linearly dependent on the 9 or so latent variables, your model is not positive definite. Add residual variance terms for the manifest variables (the diagonal of the S matrix) and the model will be identified.

Why are you using this as a saturated model?

yenni's picture
Offline
Joined: 10/08/2013
my saturated model

Thanks for your response. I'm trying on it now.
Actually, I'm trying to compare OpenMx results and LISREL results according to the paper of Statistica Neerlandica(Toharudin, Oud, Billiet 2007).
The structural model runs fine, although some parameters was not same with LISREL results.
My problem is I can not calculate chi square, so I try it to calculate Saturated Likelihood.

I have also tried with syntax below to calculate it, but it keeps running forever.
###compute Likelihood of saturated model and compare with likelihood of postulated model
Sat.Model <- function(data1){ # data1 are the data WITHOUT the time intervals at the end
SatModel <- mxModel("SatModel",
mxMatrix(type="Full", nrow=1, ncol=ncol(data1), free=TRUE, values=sapply(data1, mean, na.rm=T), name="expMean"),
mxMatrix(type="Symm", nrow=ncol(data1), ncol=ncol(data1), free=TRUE, values=var(data1, na.rm=T), name="expCov"),
mxData(observed=data1, type="raw"),
mxFIMLObjective(covariance="expCov", means="expMean", dimnames=names(data1))
)
SatModelFit<-mxRun(SatModel)
summary_SatFit<-summary(SatModelFit)
return(c(summary_SatFit$Minus2LogLikelihood, summary_SatFit$degreesOfFreedom))
}

chi_square <- function(data1, summary_fit)
{
temp <- Sat.Model(data1)
LLsat <- temp[1]
df_sat <- temp[2]

summary_fit$SaturatedLikelihood <- temp[1]
summary_fit$saturatedDoF <- temp[2]
summary_fit$Chi <- summary_fit$Minus2LogLikelihood - LLsat
chi_df <- summary_fit$degreesOfFreedom - df_sat
summary_fit$p <- (1-pchisq(summary_fit$Chi, chi_df))
summary_fit$RMSEA <- sqrt(((summary_fit$Chi/chi_df)- 1)/summary_fit$numObs )
return(summary_fit)
}

at the same time I'm also running your suggestion, but it keeps still running. I'm afraid that the result will be same with the previous one.

mhunter's picture
Offline
Joined: 07/31/2009
Try something like the

Try something like the following:

###compute Likelihood of saturated model and compare with likelihood of postulated model
Sat.Model <- function(data1){	# data1 are the data WITHOUT the time intervals at the end
	SatModel <- mxModel("SatModel",
		mxMatrix(type="Full", nrow=1, ncol=ncol(data1), free=TRUE, values=sapply(data1, mean, na.rm=T), name="expMean"),
		mxMatrix(type="Symm", nrow=ncol(data1), ncol=ncol(data1), free=TRUE, values=var(data1, na.rm=T), name="expCov"),
		mxData(observed=data1, type="raw"),
		mxFIMLObjective(covariance="expCov", means="expMean", dimnames=names(data1))
	)
	return(SatModel)
}
mySaturatedModel <- Sat.Model(myDataToFit)
mySaturatedFit <- mxRun(mySaturatedModel)
 
myComparisonModel <- mxModel(
#...
)
 
myComparisonFit <- mxRun(myComparisonModel)
 
summary(myComparisonFit, SaturatedLikelihood=mySaturatedFit)

If that does not work, then let us know where in the process it failed and we should be able to fix it!

Ryne's picture
Offline
Joined: 07/31/2009
This is a starting values

This is a starting values problem. No matter what constant value you pick for the single "variances and covariance" path, your expected covariance matrix will not be positive definite because all variables will be perfectly correlated. Instead, make one set of paths (all=TRUE) for the covariances with one constant, then make a second set of paths for just variances (all=FALSE) and give that set of paths a different, higher constant. Whenever you get this error at iteration 0 (the OpenMx backend is in C, which starts counting from zero), the problem exists in your initial starting values; higher iterations means that the optimization led to a bad expected covariance.

#Covariances
mxPath(
from=manifests,
arrows=2,
free=TRUE,
all=TRUE,
values=.5),

#Variances
mxPath(
from=manifests,
arrows=2,
free=TRUE,
all=FALSE,
values=1),

brauer's picture
Offline
Joined: 01/28/2012
Starting values and false parameter estimates

Hi,

I have to admit that I don't understand your answer. I used the script below to test a fully saturated model with four observed variables. I gave the paths a starting value of .1, the covariance a starting value of .5, and the variances a starting value of 1. The model runs but it produces an output that is different form the one I get with EQS or AMOS (different parameter estimates). Where is the problem? I tried to take out the mean structure or to rescale the variables so that they have more similar variances but this didn't help either. Any ideas?

-- Markus

data_apgar <- read.spss("data_apgar.sav", to.data.frame = T)
manifests <- c ("apgar", "gestat", "wgtgain", "smokes")
apgar_model1 <- mxModel ("Apgar Model 1",
type="RAM",
manifestVars = manifests,
mxPath(from="smokes", to="gestat", values=.1, label="b"),
mxPath(from="smokes", to="apgar", values=.1, label="c"),
mxPath(from="wgtgain", to="gestat", values=.1, label="d"),
mxPath(from="wgtgain", to="apgar", values=.1, label="e"),
mxPath(from="gestat", to="apgar", values=.1, label="f"),
mxPath(from=manifests, arrows=2, free=TRUE, values=1),
mxPath(from="smokes", to="wgtgain", arrows=2, values=.5, label="a"),
mxPath(from="one", to=manifests, arrows=1, free=TRUE, values=1),
mxData(data_apgar, type="raw")
)
apgar1out <- mxRun(apgar_model1)
apgar1out@output
summary(apgar1out)

tbates's picture
Offline
Joined: 07/31/2009
data_apgar

Any chance you can post the dataset and the AMOS/EQS results? Should be easy to solve then.

brauer's picture
Offline
Joined: 01/28/2012
Expected covariance matrix is not positive definite

Sure, here is the dataset and the Amos output. Thanks in advance for your help.

AttachmentSize
data_apgar.sav 1.43 KB
apgar1_out new.doc 81.5 KB
tbates's picture
Offline
Joined: 07/31/2009
same?

Hi there,
I ran the model here under OpenMx, and I get the same result as the one you posted from M+ to 3 decimals.

AttachmentSize
samizdat.pdf 286.58 KB
brauer's picture
Offline
Joined: 01/28/2012
graphic output

Hi again,

How did you generate the graphic representation of the model in the file you sent me ("samizdat.pdf")? Is there a way to get R to draw the graph for me? With the path coefficients? Thank you!

rgore's picture
Offline
Joined: 01/27/2011
My hunch is ...

While I don't know for sure I'd guess the "Show RAM Model in Graphviz Format" functionality of OpenMx was used. See here for how to use it, then render the file in Graphviz.

brauer's picture
Offline
Joined: 01/28/2012
Expected covariance matrix is not positive definite

Hi again,

Your help is greatly appreciated. I am not sure I know how to read the output. I pasted the output in a word document (see attached doc). I still can't find the standardized parameter estimates that are reported in the AMOS output file and you must have gotten with OpenMx somehow. For example, the path from "smokes" to "wgtgain" is -.461. Where did you find this estimate in the output? I searched the document for "-.46" (and -0.46) but didn't find it.

Then there are a number of things that are not in the output but that I probably have to ask for explicitly with the OpenMx script: The significance level of the paths, the standardized paths, the correlations between exogenous variables, the sample covariance matrix, the model-implied covariance matrix, the residual covariance matrix, the squared multiple correlations (how much variance my model explains in each of the endogenous variables). If you can tell me what I have to do to get the relevant output, that would be great. If not can you direct me to the appropriate documentation? Thanks a million in advance.

AttachmentSize
Output for apgar1out.doc 54.5 KB
tbates's picture
Offline
Joined: 07/31/2009
summary and more

Yes, you're right: the OpenMx summary gives a range of data, but you can, and should "ask for more".

Significance levels of paths are best tested with CIs (which will be included in the summary when requested in the model).
Standardization, I attach a script for (written by Ryne originally?)

Many things, though, are standard R: so sample size and count of manifests are just
dim(yourData)

Correlations and covariances, again you have R
cor(yourData)
cov(yourData)

etc.

AttachmentSize
standardizeRAM.R 3.08 KB
brauer's picture
Offline
Joined: 01/28/2012
Significance levels, standardized paths

Hi, thank you for your help. I looked at the description of MxPath and MxModel but couldn't figure out how to ask for confidence intervals in the model. Can you tell me how?

Standardization: Is there really no easier way than to write a 75-line script? My students are going to flip out if I show them the script you sent me! Getting standardized estimates seems such a frequent thing to do that I can't believe that there is an easier way to do it.

I still didn't entirely understand how to get the model-implied covariance matrix (I assume it's EC1 <- mxEval(S, apgar1out)) and the residual covariance matrix.

tbates's picture
Offline
Joined: 07/31/2009
CIs, and long functions

Hi brauer,

To compute confidence intervals, you add mxCI() statements to the model, then call mxRun with intervals=TRUE.

mxCI statements are simply a list of labels, or matrices, or matrix cell addresses to caluculate Cis on.

So, from memory, you could say "mxCI( c('A') )" to get CIs on all the asymmetric paths, or mxCI(c('a','b','c')) to get CIs on your labelled paths.

The standardization function is lengthy. Because RAM is one arbitrary type of model that OpenMx can solve, it does not come packaged with all the functions that a person focussed on RAM-specified path models would use and want.

The idea is that you can build a small library of helper functions to automate the tasks you need repeatedly, hopefully sharing those here, like the standardization function.

Each group would then package up the (relatively few) functions they need, and source() those at the top of each script.

So if you drop the stand function and any others you make up for teaching into a file, that is instantly available to all your students as if it was built-in, but customisable to your exact needs and style.

I agree that a help file from the "mothership" as it were is a good idea, and I think that will emerge as the user base grows and settles on some standard shared requirements. Your questions help with that, so thank you!

Hope this helps,
tim

dadrivr's picture
Offline
Joined: 01/19/2010
Great, that solved the

Great, that solved the problem. Do you know if there will be an alternative to "all=TRUE" for quick model specification in the next version of OpenMx? How should I code these models in the interim (i.e., omitting 'all=TRUE')? Do I need to code each set of paths separately? Thanks!

tbates's picture
Offline
Joined: 07/31/2009
connect="unique pairs"

"all=" has been replaced with "connect="

This is MUCH easier, more powerful, and easier to understand.

You can now connect from and to in the following ways:

"all.pairs": All possible combinations amongst from and to
"all.bivariate": As above, but excluding combinations from an parameter back to itself (covariances)

"unique.pairs": All possible combinations excluding mirror images (reversed from and to)
"unique.bivariate". "from" != "to" && no mirror images

"single": Just pairs up the from and to elements, as item 1 of "from" to item 1 of "to", etc.

prast's picture
Offline
Joined: 08/18/2010
connect='all.pairs'

Hi, I have code which used to run with the all=TRUE argument but now, somehow, the exact same model does not run anymore.
I changed the all=TRUE to connect='all.pairs' and arrows=2 to arrows=1 (see old values behind the comment #). I would think that this should be the equivalent:

mxPath(
from=c("L1","S1","L2","S2"), # bivariate GCM with correlated intercepts and slopes
arrows=1,
#arrows=2,
connect='all.pairs',
#all=TRUE,
free=TRUE,
values=c(SigmaDep),
labels=c("varL1", "cL1S1", "cL1L2", "cL1S2",
"cL1S1", "varS1", "cS1L2", "cS1S2",
"cL1L2", "cS1L2", "varL2", "cL2S2",
"cL1S2", "cS1S2", "cL2S2", "varS2"),
lbound=c( 0, NA, NA, NA,
NA, 0, NA, NA,
NA, NA, 0, NA,
NA, NA, NA, 0)
),

After running the model I get 'The job for model 'bivgcm' exited abnormally with the error message: Expected covariance matrix is not positive-definite in data row 122 at major iteration 0 (minor iteration 1).'
Where am I going wrong?
Somehow the labels and values now showp up in the A matrix instead of in the S matrix. Is this intended OpenMx 1.2.?

mspiegel's picture
Offline
Joined: 07/31/2009
We are deprecating the

We are deprecating the ability of mxPath() to generate redundant specifications of the same path. This is why mxPath(c('a','b'), connect = 'all.pairs', arrows = 2) is not permitted, because the path from 'b' to 'a' is a redundant specification of the path from 'a' to 'b'. You can rewrite your path specification using connect = "unique.pairs" and arrows = 2:

mxPath(from = c("L1","S1","L2","S2"), 
    arrows = 2,
    connect = "unique.pairs",
    free = TRUE,
    labels = c("varL1", "cL1S1", "cL1L2", "cL1S2",
        "varS1", "cS1L2", "cS1S2", "varL2", "cL2S2",
        "varS2"))

You'll need to determine how to re-write your 'values' and 'lbound' arguments. Another alternative is to split this into two constructions, using "unique.bivariate" and "single". Sometimes it is easier to specify the 'values' and other arguments with this alternative construction:
mxPath(from = c("L1","S1","L2","S2"),
arrows = 2, connect = "unique.bivariate",
free = TRUE, labels = c("cL1S1", "cL1L2", "cL1S2",
"cS1L2", "cS1S2", "cL2S2"))
 
mxPath(from = c("L1","S1","L2","S2"),
arrows = 2, connect = "single",
free = TRUE, labels = c("varL1", "varS1", "varL2", "varS2"))

prast's picture
Offline
Joined: 08/18/2010
Thank you for your swift

Thank you for your swift reply!
The second approach did the trick.

tbates's picture
Offline
Joined: 07/31/2009
connect = "unique.pairs"

The problem (and the change from S (ymmetric) to A (symmetric) is because you have changed your paths from 2 to 1 arrow. You don't need (or want) to do that.

To replicate what all did by default but without the illegal mirror-pairs, just say:

 connect = "unique.pairs"

To link each source to everything except itself, and exclude equivalent paths, say:
 connect = "unique.bivariate"

The full option list is:

 connect = c("single", "all.pairs", "unique.pairs",  "all.bivariate", "unique.bivariate"),