# Modification of oneACEvca.R # by Hermine Maes # # Twin Univariate ACE model to estimate causes of variation across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries & Options library(OpenMx) library(psych) # Create Output filename <- "oneACEvca" sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE) # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data data(twinData) dim(twinData) describe(twinData[,1:12], skew=F) # Select Variables for Analysis vars <- 'bmi' # 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 twinData[,'age'] <- twinData[,'age']/100 twinData <- twinData[-which(is.na(twinData$age)),] covVars <- 'age' # Select Data for Analysis mzData <- subset(twinData, zyg==1, c(selVars, covVars)) dzData <- subset(twinData, zyg==3, c(selVars, covVars)) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") cor(mzData,use="complete") cor(dzData,use="complete") # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 20 # start value for means svPa <- .2 # start value for path coefficient svPe <- .5 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defL" ) pathBl <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="b11", name="bl" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels="b0", name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(defL%*%bl,defL%*%bl), name="expMean" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="A11", name="A" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="C11", name="C" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="E11", 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" ) corMZ <- mxAlgebra( expression= A+C, name="rMZ" ) corDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="rDZ" ) 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" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( pathBl, meanG, covA, covC, covE, covP ) defs <- list( defL ) modelMZ <- mxModel( pars, defs, expMean, covMZ, corMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMean, covDZ, corDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) 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(c("VC","MZ.rMZ","DZ.rDZ")) # Build Model with Confidence Intervals modelACE <- mxModel( "oneACEvca", pars, modelMZ, modelDZ, multi, estVC, ciACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=F ) sumACE <- summary( fitACE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run AE model modelAE <- mxModel( fitACE, name="oneAEvca" ) modelAE <- omxSetParameters( modelAE, labels="C11", free=FALSE, values=0 ) fitAE <- mxRun( modelAE, intervals=F ) # Run CE model modelCE <- mxModel( fitACE, name="oneCEvca" ) modelCE <- omxSetParameters( modelCE, labels="A11", free=FALSE, values=0 ) modelCE <- omxSetParameters( modelCE, labels=c("E11","C11"), free=TRUE, values=.6 ) fitCE <- mxRun( modelCE, intervals=F ) # Run E model modelE <- mxModel( fitAE, name="oneEvca" ) modelE <- omxSetParameters( modelE, labels="A11", free=FALSE, values=0 ) fitE <- mxRun( modelE, intervals=T ) # Print Comparative Fit Statistics mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4) #Calculate confidence intervals whichever models' confidence intervals you want to report: fitACE <- omxRunCI(fitACE) #<--e.g., if you want to report CIs from 'fitACE' summary(fitACE) # ---------------------------------------------------------------------------------------------------------------------- sink()