# ------------------------------------------------------------------------------ # Program: twinAceSibs.R # Author: Sarah Medland # Date: 13.04.2012 # # Univariate Twin Analysis model to estimate causes of variation # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------| # Load Library require(OpenMx) # ------------------------------------------------------------------------------ # PREPARE DATA # Load Data mzData <- read.csv("mzData.csv") dzData <- read.csv("dzData.csv") # Select Variables for Analysis nv <- 2 # number of variables nind <- 3 ntv <- nv*nind # number of total variables selVars <- c("Score_t1", "Disease_t1", "Score_t2", "Disease_t2", "Score_s1", "Disease_s1") defVars <- c("ageTwins", "ageSib") useVars <- c(selVars,defVars) # Set the Binary variable up for OpenMx mzData$Disease_t1 <- mxFactor(mzData$Disease_t1, levels=c(0:1) ) mzData$Disease_t2 <- mxFactor(mzData$Disease_t2, levels=c(0:1) ) mzData$Disease_s1 <- mxFactor(mzData$Disease_s1, levels=c(0:1) ) dzData$Disease_t1 <- mxFactor(dzData$Disease_t1, levels=c(0:1) ) dzData$Disease_t2 <- mxFactor(dzData$Disease_t2, levels=c(0:1) ) dzData$Disease_s1 <- mxFactor(dzData$Disease_s1, levels=c(0:1) ) # ------------------------------------------------------------------ # PREPARE MODEL # Age correction # Regression effects Beta1 <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=.1, labels=c('betaAgeD', 'betaAgeS'), name="beta1" ) # Independent variable obsAge <- mxMatrix( type="Full", nrow=1, ncol=3, free=FALSE, labels=c("ageTwins", "ageTwins", "ageSib"), name="Age") #setting up the regression intercept <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=c(T,F), values=c(3,0), labels=c("meanP","binary"), name="intercept" ) threshold <-mxMatrix( type="Full", nrow=1, ncol=nind, free=T, values=1, labels="thresh", name="Threshold" ) expMean <- mxAlgebra( intercept + (Age %x% beta1), name="expMean") inclusions <- list (Beta1, obsAge, intercept, expMean, threshold) # ------------------------------------------------------------------------------ # ACE Model # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label=c("a11","a21","a22"), name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label=c("c11","c21","c22"), name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,F), values=1, label=c("e11","e21","e22"), name="e" ) # Matrices generated to hold A, C, and E computed 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" ) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C, .5%x%A+C), cbind(A+C, A+C+E, .5%x%A+C), cbind(.5%x%A+C, .5%x%A+C, A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , .5%x%A+C, .5%x%A+C), cbind(.5%x%A+C, A+C+E, .5%x%A+C), cbind(.5%x%A+C, .5%x%A+C, A+C+E)), name="expCovDZ" ) Imat <- mxMatrix(name= "I", type="Iden", nrow = nv, ncol = nv) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="Threshold", threshnames=c('Disease_t1','Disease_t2','Disease_s1') ) objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars , thresholds="Threshold", threshnames=c('Disease_t1','Disease_t2','Disease_s1') ) # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE , Imat) modelMZ <- mxModel( pars, inclusions, covMZ, dataMZ, objMZ, name="MZ" ) modelDZ <- mxModel( pars, inclusions, covDZ, dataDZ, objDZ, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxAlgebraObjective( "m2LL" ) AceModel <- mxModel( "ACE", pars, modelMZ, modelDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE model AceFit <- mxRun(AceModel) AceSumm <- summary(AceFit) AceSumm round(AceFit@output$estimate,4) # To get standardized results iSDA <- mxEval( solve (sqrt(I*A) ), AceFit) (corA <- mxEval( iSDA %*% A %*% iSDA, AceFit)) iSDC <- mxEval( solve (sqrt(I*C) ), AceFit) (corC <- mxEval( iSDC %*% C %*% iSDC , AceFit)) iSDE <- mxEval( solve (sqrt(I*E) ), AceFit) (corE <- mxEval( iSDE %*% E %*% iSDE, AceFit)) covP <- mxEval( A+C+E , AceFit ) iSDP <- mxEval( solve (sqrt(I*covP) ), AceFit) (corP <- mxEval( iSDP %*% covP %*% iSDP, AceFit )) (standA <- mxEval( iSDP %*% a , AceFit)) (standC <- mxEval( iSDP %*% c , AceFit)) (standE <- mxEval( iSDP %*% e , AceFit)) # Generate ACE Model Output estVA <- mxEval(A, AceFit) # additive genetic variance, a^2 estVC <- mxEval(C, AceFit) # shared environmental variance, c^2 estVE <- mxEval(E, AceFit) # unique environmental variance, e^2 estVP <- (estVA+estVC+estVE) # total variance estPropVA <- estVA/estVP # standardized additive genetic variance estPropVC <- estVC/estVP # standardized shared environmental variance estPropVE <- estVE/estVP # standardized unique environmental variance LL_ACE <- mxEval(objective, AceFit) # likelihood of ACE model # ------------------------------------------------------------------------------ # FIT SUBMODELS # Run AE model AeModel <- mxModel( AceFit, name="AE" ) AeModel <- omxSetParameters( AeModel, labels="c11", free=FALSE, values=0 ) AeModel <- omxSetParameters( AeModel, labels="c21", free=FALSE, values=0 ) AeModel <- omxSetParameters( AeModel, labels="c22", free=FALSE, values=0 ) AeFit <- mxRun(AeModel) round(AeFit@output$estimate,4) # To get standardized results iSDA <- mxEval( solve (sqrt(I*A) ), AeFit) corA <- mxEval( iSDA %*% A %*% iSDA, AeFit) iSDE <- mxEval( solve (sqrt(I*E) ), AeFit) corE <- mxEval( iSDE %*% E %*% iSDE, AeFit) covP <- mxEval( A+C+E , AeFit ) iSDP <- mxEval( solve (sqrt(I*covP) ), AeFit) corP <- mxEval( iSDP %*% covP %*% iSDP , AeFit) standA <- mxEval( iSDP %*% a , AeFit) standE <- mxEval( iSDP %*% e , AeFit) # Generate AE Model Output estVA <- mxEval(a*a, AeFit) # additive genetic variance, a^2 estVE <- mxEval(e*e, AeFit) # unique environmental variance, e^2 estVP <- (estVA+estVE) # total variance estPropVA <- estVA/estVP # standardized additive genetic variance estPropVE <- estVE/estVP # standardized unique environmental variance LL_AE <- mxEval(objective, AeFit) # likelihood of AE model LRT_ACE_AE <- LL_AE - LL_ACE # likelihood ratio test ACE-AE model # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics source("GenEpiHelperFunctions.R") tableFitStatistics(AceFit,AeFit)