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