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)