You are here

ACE model with continuous covariate (different between twins)

2 posts / 0 new
Last post
Kelsey's picture
Offline
Joined: 01/31/2017 - 11:34
ACE model with continuous covariate (different between twins)

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)
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
There's a problem in this

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" )
Log in or register to post comments