A Simple ACE Model that Seems to Be Failing...?
Jeremy
Joined: 03/25/2016
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.
betas not identified
mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars ) ),
and
mxFIMLObjective(covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars)
In the namespace of your MxModel object (which has R symbol
twinACEModel),"ACE.expMean"refers to the object created bymxAlgebra(cbind(M,M), name="expMean")
which, 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 tomxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars ) ),
and
mxFIMLObjective(covariance="ACE.expCovDZ", means="expMean", dimnames=selVars)
Log in or register to post comments
Thank you!!
Log in or register to post comments