Hi, I'm trying to run a simple ACE Model that controls for age and sex. Searched the boards and found some recent posts that seem similar to what I want to do: https://openmx.ssri.psu.edu/thread/4102
I tried to adapt it and this is what I've come up with:
# Select Variables for Analysis Vars <- c('caster') nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- c("caster1", "caster2") mydata$ageT1[ is.na(mydata$ageT1) ] <- -999 mydata$ageT2[ is.na(mydata$ageT2) ] <- -999 mydata$sexT1[ is.na(mydata$sexT1) ] <- -999 mydata$sexT2[ is.na(mydata$sexT2) ] <- -999 mydata[mydata$ageT1==-999,]$caster1 <- NA mydata[mydata$ageT2==-999,]$caster2 <- NA mzData <- subset(mydata, ZYGOSITY=="0") dzData <- subset(mydata, ZYGOSITY=="1") twinACEModel <- mxModel("twinACE", mxModel("ACE", # ace path coefficients for males and females mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.6, label="a11", name="a"), mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="c11", name="c"), mxMatrix("Lower", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="e11", name="e"), # ACE variance components for males and females mxAlgebra(a %*% t(a), name="A"), mxAlgebra(c %*% t(c), name="C"), mxAlgebra(e %*% t(e), name="E"), # Declare a matrix for the definition variable regression parameters mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("bage_caster"), name="betaage"), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= 0, label=c("bsex_caster"), name="betasex"), # Total variation mxAlgebra(A+C+E, name="V"), # Standard deviations for males and females mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="isd"), # Standardized path estimates mxAlgebra(a %*% isd, name="sta"), mxAlgebra(c %*% isd, name="stc"), mxAlgebra(e %*% isd, name="ste"), # Standardized variance components mxAlgebra(A/V, name="h2"), mxAlgebra(C/V, name="c2"), mxAlgebra(E/V, name="e2"), # Expected means vectors: male-male pairs; fem-fem pairs; male-fem opposite-sex pairs mxMatrix("Full", nrow=1, ncol=nv, free=TRUE, values= .6, label="mean", name="M"), mxAlgebra(cbind(M,M), name="expMean"), # Expected covariance matrix in MZ mxAlgebra( rbind( cbind(A+C+E, A+C), cbind(A+C, A+C+E) ), name="expCovMZ" ), # Expected covariance matrix in DZ mxAlgebra( rbind( cbind(A+C+E, (0.5 %x% A)+C), cbind((0.5 %x% A)+C, A+C+E) ), name="expCovDZ" ) ), mxModel("MZ", mxData(observed=mzData , type="raw"), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression= ACE.expMean + ACE.betasex %x%obs_sex + ACE.betaage %x%obs_age, name="expMean"), mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars ) ), mxModel("DZ", mxData(observed=dzData , type="raw"), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ), mxAlgebra( expression= ACE.expMean + ACE.betasex %x%obs_sex + ACE.betaage %x%obs_age, name="expMean"), mxFIMLObjective(covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars) ), mxAlgebra(MZ.objective + DZ.objective, name="minus2sumll"), mxAlgebraObjective("minus2sumll"), mxCI(c('ACE.h2', 'ACE.c2', 'ACE.e2')), mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')), mxCI(c('ACE.expCovMZ', 'ACE.expCovMZ')) ) twinACEFit <- mxRun(twinACEModel) twinACEFit <- mxRun(twinACEModel, intervals=T)
The issue is that when I look at the output, this is what it says:
Summary of twinACE free parameters: name matrix row col Estimate Std.Error A 1 a11 ACE.a 1 1 0.57664153 39.03038 2 c11 ACE.c 1 1 0.63788121 81.91912 3 e11 ACE.e 1 1 0.51015713 NA 4 bage_caster ACE.betaage 1 1 0.10000000 NA * 5 bsex_caster ACE.betasex 1 1 0.10000000 175.59616 * 6 mean ACE.M 1 1 0.00770848 78.21110
Eh? Why is the value exactly .1 for bage and bsex? And no standard error for one...?
I've looked over this multiple times and have to be missing something obvious.
The problem is with these two lines:
and
In the namespace of your MxModel object (which has R symbol
twinACEModel
),"ACE.expMean"
refers to the object created bywhich, you'll notice, contains only intercepts, and does not involve the regression coefficients for age and sex. Because you define the correct algebra for the phenotypic means,
mxAlgebra( expression= ACE.expMean + ACE.betasex %x%obs_sex + ACE.betaage %x%obs_age, name="expMean")
, in the same submodels as themxFIMLObjective()
statements, change those statements respectively toand
Thank you!!