You are here

A Simple ACE Model that Seems to Be Failing...?

3 posts / 0 new
Last post
Jeremy's picture
Offline
Joined: 03/25/2016 - 13:07
A Simple ACE Model that Seems to Be Failing...?

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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
betas not identified

The problem is with these two lines:

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 by

mxAlgebra(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 the mxFIMLObjective() statements, change those statements respectively to

mxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars )  ),

and

mxFIMLObjective(covariance="ACE.expCovDZ", means="expMean", dimnames=selVars)
Jeremy's picture
Offline
Joined: 03/25/2016 - 13:07
Thank you!!

Thank you!!