ACE model with continuous covariate (different between twins)

Posted on
No user picture. Kelsey Joined: 01/31/2017
Helllo,

I’m very new at OpenMx but I’m trying my hand at creating an ACE model with one dichotomous variable and a continuous covariate, which differs between both twins. I am trying to adapt a script (http://ibg.colorado.edu/cdrom2016/maes/MultivariateAnalysis/oneACEba.R) which uses age as a covariate. But I can’t seem to adapt this script to include a continuous covariate for both twins without the following error:
“Error: The following error occurred while evaluating the subexpression 'MZ.b %*% MZ.BMI' during the evaluation of 'MZ.expMean' in model 'oneACEba' : non-conformable arguments.”

Could anyone help me to correctly adapt the script below?

Thanks!
Kelsey


# Select Variables for Analysis
vars <- 'hs' # list of variables names
nv <- 1 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Covariates for Analysis
varscov <- 'bmi' # list of covariate names
covVars <- paste(varscov,c(rep(1,nv),rep(2,nv)),sep="")

# Select Data for Analysis
mzData <- subset(twinData, zyg==1, c(selVars, covVars))
dzData <- subset(twinData, zyg==2, c(selVars, covVars))
mzDataF <- cbind(mxFactor( x=mzData[,selVars], levels=c(0:1)), mzData[,covVars])
dzDataF <- cbind(mxFactor( x=dzData[,selVars], levels=c(0:1)), dzData[,covVars])

# Set Starting Values
svTh <- .8 # start value for thresholds
svPa <- .4 # start value for path coefficient
svPe <- .8 # start value for path coefficient for e
lbPa <- .0001 # start value for lower bounds

# ------------------------------------------------------------------------------
# PREPARE MODEL

# ACE Model
# Create Matrices for Covariates and linear Regression Coefficients
defBMI <- mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c("data.bmi1","data.bmi2"), name="BMI" )
pathB <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.01, label=c("b11", "b12"), name="b" )

# Create Algebra for expected Mean & Threshold Matrices
meanG <- mxMatrix( type="Zero", nrow=1, ncol=ntv, name="meanG" )
expMean <- mxAlgebra( expression= meanG + cbind(b%*%BMI),b%*%BMI), name="expMean" )
threG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svTh, labels="tob", name="threG" )

# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", lbound=lbPa, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="c11", lbound=lbPa, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="e11", lbound=lbPa, name="e" )

# Create Algebra for 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" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+C+E, name="V" )
covMZ <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" )

# Constrain Variance of Binary Variables
var1 <- mxConstraint( expression=diag2vec(V)==1, name="Var1" )

# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzDataF, type="raw" )
dataDZ <- mxData( observed=dzDataF, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG" )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG" )
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars <- list(pathB, meanG, threG, pathA, pathC, pathE, covA, covC, covE, covP)
defs <- list(defBMI)
modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Create Algebra for Variance Components
rowVC <- rep('VC',nv)
colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))

# Create Confidence Interval Objects
ciACE <- mxCI( "VC[1,1:3]" )

# Build Model with Confidence Intervals
modelACE <- mxModel( "oneACEba", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE )

# ------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

# Compare with Saturated Model
mxCompare( fit, fitACE )
#lrtSAT(fitACE,-2ll,df)

# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEsts(fitACE)

# ------------------------------------------------------------------------------
# RUN SUBMODELS

# Test Significance of Covariate
modelCov <- mxModel( fitACE, name="testCov" )
modelCov <- omxSetParameters( modelCov, label=c("b11","b12" ), free=FALSE, values=0 )
fitCov <- mxRun( modelCov )
mxCompare( fitACE, fitCov )

# Run AE model
modelAE <- mxModel( fitACE, name="oneAEba" )
modelAE <- omxSetParameters( modelAE, labels="c11", free=FALSE, values=0 )
fitAE <- mxRun( modelAE, intervals=T )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)

# Run CE model
modelCE <- mxModel( fitACE, name="oneCEba" )
modelCE <- omxSetParameters( modelCE, labels="a11", free=FALSE, values=0 )
fitCE <- mxRun( modelCE, intervals=T )
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
fitEsts(fitCE)

# Run E model
modelE <- mxModel( fitAE, name="oneEba" )
modelE <- omxSetParameters( modelE, labels="a11", free=FALSE, values=0 )
fitE <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)

# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitCov, fitAE, fitCE, fitE) )
round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4)

Replied on Fri, 09/08/2017 - 11:29
Picture of user. AdminRobK Joined: 01/24/2014

There's a problem in this line,

expMean <- mxAlgebra( expression= meanG + cbind(b%*%BMI),b%*%BMI), name="expMean" )

which I'm surprised even ran. Try this instead:

expMean <- mxAlgebra( expression= meanG + b*BMI, name="expMean" )

Also, are you sure want to allow the effect of the covariate to be different for twins #1 and twins #2? If you don't, redefine pathB as:

pathB <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.01, labels="b1", name="b" )
Replied on Tue, 10/03/2017 - 16:20
No user picture. FrankJ Joined: 09/11/2017

In this script, only the main effect was examined regarding if BMI (the covariate) affects the phenotype(hs in this case). If we also want to assess GxE for BMI, how to improve the script? Do we have to constrain the variance of twin 1 equals variance of twin 2? For example, (a+am*BMI_1)^2+(c+cm*BMI_1)^2+(e+em*BMI_1)^2 == (a+am*BMI_2)^2+(c+cm*BMI_2)^2+(e+em*BMI_2)^2

Best,
Frank

Replied on Wed, 10/04/2017 - 10:59
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by FrankJ

The go-to reference for biometrical moderation modeling when the putative moderator may differ for each twin in a pair is:
van der Sluis, S., Posthuma, D., & Dolan, C. V. (2012). A note on false positives and power in G × E modelling of twin data. Behavior Genetics, 42(1), 170-186. Also see here for relevant slides from the 2016 Boulder Workshop.

The actual implementation in OpenMx requires that you copy the two columns (one for each twin) containing the putative moderator into your dataset again. In your case, you'd initially have a column for BMI for twin #1, and BMI for twin #2. Copy each column so that it appears twice in your dataset. One copy of each column will be used as an endogenous variable of the model (sort of like a second phenotype), and the other copy will be used as a definition variable.

Replied on Thu, 10/05/2017 - 11:15
No user picture. FrankJ Joined: 09/11/2017

Thanks very much. I will read and try this and may come back if have questions.