Multivariate ACE error message in binary data

1 reply [Last post]
romidico's picture
Offline
Joined: 01/09/2013

Hi all,

I am trying to run multivariate ACE model for binary data, but I keep on getting error message...

This is my script:

require(OpenMx)

# ------------------------------------------------------------------------------
# PREPARE DATA

OMdataset <- read.csv("twins-only.csv", header=T,)
# sep = ",", dec = ".", quote="\"", na.strings="-99"

OMdataset$TONSa -> OMdataset$T_01
OMdataset$TONSb -> OMdataset$T_02

OMdataset$ADENa -> OMdataset$A_01
OMdataset$ADENb -> OMdataset$A_02

OMdataset$GROMa -> OMdataset$G_01
OMdataset$GROMb -> OMdataset$G_02

# Set the Binary variable up for OpenMx
OMdataset$A_01 <- mxFactor(OMdataset$A_01, levels=c(0:1) )
OMdataset$A_02 <- mxFactor(OMdataset$A_02, levels=c(0:1) )
OMdataset$G_01 <- mxFactor(OMdataset$G_01, levels=c(0:1) )
OMdataset$G_02 <- mxFactor(OMdataset$G_02, levels=c(0:1) )
OMdataset$T_01 <- mxFactor(OMdataset$T_01, levels=c(0:1) )
OMdataset$T_02 <- mxFactor(OMdataset$T_02, levels=c(0:1) )
OMdataset$agea[OMdataset$agea==0] <- 12
OMdataset$ageb[OMdataset$ageb==0] <- 12
OMdataset$SEXa[ OMdataset$Zyg==6] <- 1
OMdataset$agea <- OMdataset$agea/100
OMdataset$ageb <- OMdataset$ageb/100
# Select Variables for Analysis
selVars <- c('A_01' ,'G_01' ,'T_01' , 'A_02' ,'G_02' ,'T_02' )
defVars <- c('SEXa','SEXb','agea','ageb')
useVars <- c('A_01' ,'G_01' ,'T_01' , 'A_02' ,'G_02' ,'T_02' ,
'SEXa','SEXb','agea','ageb')
# ------------------------------------------------------------------

# PREPARE MODEL
stThresh <-1.1
nv <-3
ntv <-6
# Sex and Age corrections
# Regression effects
beta1 <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE,
values=.1, labels=c('betaSexA','betaSexG','betaSexT'), name="beta1" )
# Independent variables
obsSex <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE,
labels=c("data.SEXa","data.SEXa","data.SEXa","data.SEXb","data.SEXb","data.SEXb" ), name="Sex")

#Means - fixed
mean <- mxMatrix( type="Zero", nrow=1, ncol=6, name="mean" )
thresholds <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=1,
labels=c("threshA","threshG","threshT"), name="threshold" )

expthresh <- mxAlgebra( threshold + beta1 * Sex, name="expThresh")

inclusions <- list (beta1, obsSex, mean, expthresh, thresholds)

# Select Data for Analysis
mzData <- subset(OMdataset, Zyg<=2, useVars)

dzData <- subset(OMdataset, Zyg>=3, useVars)
# ------------------------------------------------------------------------------

# ACE Model
# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.8,.5,.5,.4,.3,.4),
labels=c("a11","a21","a31","a22","a32","a33"), name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.2,
labels=c("c11","c21","c31","c22","c32","c33"), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(.4,.1,.1,.4,.1,.4),
labels=c("e11","e21","e31","e22","e32","e33"), name="e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
totalvariance <- mxAlgebra( A+C+E , name="V" )

# Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covMZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )

covDZ <- mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )
one <- mxMatrix( type="Unit", nrow=3, ncol=1, name="one" )
var1 <- mxConstraint( diag2vec(V)==one, name="Var1" )

# 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="mean", dimnames=selVars , thresholds="expThresh")
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="mean", dimnames=selVars , thresholds="expThresh")

# Combine Groups
pars <- list(one, pathA, pathC, pathE, covA, covC, covE, totalvariance, var1 )
modelMZ <- mxModel( pars, inclusions, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, inclusions, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
CI <- mxCI(c('ACE.a','ACE.c', 'ACE.e','ACE.A','ACE.C', 'ACE.E'))
AceModel <- mxModel( "ACE", pars, modelMZ, modelDZ, minus2ll, obj, CI )

# ------------------------------------------------------------------------------
# RUN MODEL

# Run ACE model
AceFit <- mxRun(AceModel, intervals=FALSE)
try2 <- AceFit
try2fit <-mxRun(try2, intervals=FALSE)
AceSumm <- summary(AceFit)
AceSumm
round(AceFit@output$estimate,4)

When it gets to
# Run ACE model
AceFit <- mxRun(AceModel, intervals=FALSE)

It takes ages to run, and at the end this is the error message I get:

Running ACE
Warning message:
In model 'ACE' NPSOL returned a non-zero status code 6. The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)
> AceSumm <- summary(AceFit)AceSumm
Error: unexpected symbol in "AceSumm <- summary(AceFit)AceSumm"

I think it has something to do with threshold values, but I don't even know what threshold values to put...

Someone, plz help!!!!!!!

Thanks

Rob

neale's picture
Offline
Joined: 07/31/2009
Check Model Identification, and rerun

This example takes ages because there are 12 variables (6 per twin) and speedups such as grouping like observations for single evaluation don't help due to the use of definition variables. Using a cluster or multiple cores on a laptop/desktop system might help.

There has been much discussion of Code 6 - which may be a false alarm - on this site. Perhaps one of the clearest descriptions was by Ryne in this thread http://openmx.psyc.virginia.edu/thread/1659 A good strategy here is to rerun the script from different starting values, or simply to re-run from the solution:

AceFit2 <- mxRun(AceFit)

and see if the error code goes away.

The Error, unexpected symbol... is just an R error because the AceSumm should not be on the end of the line
AceSumm <- summary(AceFit)AceSumm