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