You are here

Multivariate ACE error message in binary data

2 posts / 0 new
Last post
romidico's picture
Offline
Joined: 01/09/2013 - 03:10
Multivariate ACE error message in binary data

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 - 15:14
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