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)
There's a problem in this line,
which I'm surprised even ran. Try this instead:
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: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+amBMI_1)^2+(c+cmBMI_1)^2+(e+emBMI_1)^2 == (a+amBMI_2)^2+(c+cmBMI_2)^2+(e+emBMI_2)^2
Best,
Frank
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.
You mind find the attached demo script useful.
Thanks very much. I will read and try this and may come back if have questions.