Hi all,

I’m trying to run a saturated model with restrictions in a multivarite analysis with 6 continuous variables.

In this saturated model with restrictions (see below) I had the following problem (that I did not have in the saturated model without restrictions):

Error: The job for model 'Con' exited abnormally with the error message: Expected covariance matrix is not positive-definite in data row 431 at major iteration 0 (minor iteration 1).

However, when I use the same script but instead of 6 with 5 variables or less, the script runs perfectly. I don’t know what it’s happing… I think it’s something related to the starting values, but I don’t know how I can solve this.

Any help would be very helpful.

Many thanks in advance.

#######################################################################

# Restrictions: means and variances equated across birth-order & zygosity groups;

# One set of Cross-trait cor, symmetric Cross-trait Cross-twin correlations (MZ and DZ)

# Create start values

# To avoid starting the optimization at the solution add some random noise

jiggle <-rnorm((nv*(nv+1)/2), mean = 0, sd=.1) #calculo jiggle necesario segun estos tamaños

Stmean <-colMeans(Data,na.rm=TRUE)

Stsd <-sd(Data,na.rm=TRUE)+jiggle[1:nv]

StWithinperson <-vechs(cor(Data,use="complete")) +jiggle[1:(nv*(nv-1)/2)]

StBetweenMZ <-vech(cov(mzData[,1:nv],mzData[,(nv+1):(2*nv)],use="complete")) +jiggle

StBetweenDZ <-vech(cov(dzData[,1:nv],dzData[,(nv+1):(2*nv)],use="complete")) +jiggle

# Create Labels for Column and Diagonal Matrices

mLabs <- paste("m",1:nv,sep="")

sdLabs <- paste("sd",1:nv,sep="")

# Create Lables for a Correlation Matrix

rphLabs <- paste("r",1:ncor,sep="")

# Create Labels for Lower Triangular Matrices

MZbLabs <- paste("mz", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

DZbLabs <- paste("dz", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="")

# Specify Matrices to hold parameter estimates for MZ and DZ data

# elements for one overall Means

meansG <-mxMatrix("Full", 1, ntv, free = TRUE, values = Stmean, labels=c(mLabs,mLabs), name="ExpMeans")

# elements for the SD

sdZ <-mxMatrix("Zero", nv, nv, free=F, name="padding")

sdG <-mxMatrix("Diag", nv, nv, free = TRUE, values = Stsd, labels=c(sdLabs), name="SD")

sdT <-mxAlgebra(rbind(cbind(SD,padding), cbind(padding, SD)), name="SDtwin")

# elements for the correlations

Rph <-mxMatrix("Stand", nv, nv, free = TRUE, values = StWithinperson, labels=rphLabs, name="within")

MZb <-mxMatrix("Symm", nv, nv, free = TRUE, values = StBetweenMZ, labels=MZbLabs, name="BetweenMZ")

DZb <-mxMatrix("Symm", nv, nv, free = TRUE, values = StBetweenDZ, labels=DZbLabs, name="BetweenDZ")

corMZ <-mxAlgebra(rbind(cbind(within,BetweenMZ), cbind(BetweenMZ, within)), name="RMZ")

corDZ <-mxAlgebra(rbind(cbind(within,BetweenDZ), cbind(BetweenDZ, within)), name="RDZ")

# generate expected covariance matrices

covMZ <-mxAlgebra(SDtwin %*% RMZ %*% t(SDtwin), name="ExpCovMZ")

covDZ <-mxAlgebra(SDtwin %*% RDZ %*% t(SDtwin), name="ExpCovDZ")

# Data objects for Multiple Groups

dataMZ <- mxData( observed=mzData, type="raw" )

dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups

objMZ <- mxFIMLObjective( covariance="ExpCovMZ", means="ExpMeans", dimnames=selVars )

objDZ <- mxFIMLObjective( covariance="ExpCovDZ", means="ExpMeans", dimnames=selVars )

# Combine Groups

pars <- list(meansG, sdZ, sdG, sdT, Rph)

modelMZ <- mxModel(pars, MZb, corMZ, covMZ, dataMZ, objMZ, name="MZ" )

modelDZ <- mxModel(pars, DZb, corDZ, covDZ, dataDZ, objDZ, name="DZ" )

minus2ll <- mxAlgebra(expression=MZ.objective + DZ.objective, name="m2LL" )

obj <- mxAlgebraObjective( "m2LL" )

ConCorModel <- mxModel("Con", pars, modelMZ, modelDZ, minus2ll, obj )

ConCorFit <- mxRun(ConCorModel)

ConCorSumm <- summary(ConCorFit)

Hi Clara

It's hard to figure out what is wrong with this one without the data. I recommend that you mxRun with unsafe-T and look at the expected covariance matrices (ConCorFit$MZ.ExpCovMZ and the DZ equivalent). You could also use mxEval() to examine the matrices for a specific data row (e.g. 431) in case it is something to do with the definition variables.

Hi Neale,

Thanks a lot for your help.

With the covariance matrices I extract the eigenvectors and the eigenvalues. See below the “numbers” for the MZ and DZ groups. I think because there are some negative values, that’s why the program says “Expected covariance matrix is not positive-definite”. I tried to add an “arbitrary” number to the diagonal of the initial MZ covariance matrix in order to try to define it positive, but there is still an error. So, I don’t not what else I can do. Maybe it’s very easy but I’m just learning to perform these analyses and I don’t know how to work with mxMatrix and mxAlgebra functions properly.

-------------------------------------------------------

resultsMZ <- mxEval(MZ.ExpCovMZ, ConCorModel, compute=TRUE, defvar.row = 426, cache = new.env(parent = emptyenv()), cacheBack=TRUE)[[1]]

resultsMZ

eigen(resultsMZ)$values

resultsDZ <- mxEval(DZ.ExpCovDZ, ConCorModel, compute=TRUE, defvar.row = 426, cache = new.env(parent = emptyenv()), cacheBack=TRUE)[[1]]

resultsDZ

> eigen(resultsMZ)$values

[1] 7.88408528 1.55588235 1.35711621 1.13959604 0.35107748 0.28868578 0.27918799 0.14503848 0.11393085 0.01127608 -0.05759287

[12] -0.51108037

> eigen(resultsDZ)$values

[1] 5.875105837 2.952744327 1.381547024 1.291366547 0.387578154 0.320333119 0.258804746 0.144368765 0.063997104 -0.006279386

[11] -0.009736677 -0.102626269

------------------------------------------------------------------

Some ideas?? Thanks a lot!