Multivariate Analysis: Saturated model with restrictions

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)
Difficult
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.
Log in or register to post comments
In reply to Difficult by neale
Multivariate analysis, saturated model with restrictions
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!
Log in or register to post comments