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):(2nv)],use="complete")) +jiggle
StBetweenDZ <-vech(cov(dzData[,1:nv],dzData[,(nv+1):(2nv)],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)
Printer-friendly version
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!