# ------------------------------------------------------------------------------ # Program: twins.R # Author: Hermine Maes # Date: 03 03 2014 # # Univariate Twin Saturated model to estimate means and (co)variances across multiple groups # Univariate Twin Analysis model to estimate causes of variation _ ADE model # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------| # Load Libraries require(OpenMx) require(psych) # -------|---------|---------|---------|---------|---------|---------|---------| # PREPARE DATA # Load Data data(twinData) dim(twinData) describe(twinData, skew=F) # Select Variables for Analysis Vars <- 'bmi' nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('bmi1','bmi2') # 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)) # Set Starting Values svMe <- 20 # start value for means svVa <- .8 # start value for variance lbVa <- .0001 # start value for lower bounds svVas <- diag(svVa,ntv,ntv) # start values on diagonal of covariance matrix lbVas <- diag(lbVa,ntv,ntv) # lower bound for variances laMeMZ <- c("m1MZ","m2MZ") # labels for means for MZ twins laMeDZ <- c("m1DZ","m2DZ") # labels for means for DZ twins laVaMZ <- c("v1MZ","c21MZ","v2MZ") # labels for (co)variances for MZ twins laVaDZ <- c("v1DZ","c21DZ","v2DZ") # labels for (co)variances for DZ twins # ------------------------------------------------------------------------------ # PREPARE SATURATED MODEL # mxMatrix( type, nrow, ncol, free, values, labels, lbound, ubound, name) TRiCkFiVeLaBsNow # Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="b11", name="b" ) # Algebra for expected Mean Matrices in MZ & DZ twins meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=laMeMZ, name="meanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=laMeDZ, name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(b%*%Age,b%*%Age), name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(b%*%Age,b%*%Age), name="expMeanDZ" ) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins expCovMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=laVaMZ, name="expCovMZ" ) expCovDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=svVas, lbound=lbVas, labels=laVaDZ, name="expCovDZ" ) # 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="expMeanMZ", dimnames=selVars ) objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars ) # Combine Groups modelMZ <- mxModel( "MZ", defAge, pathB, meanMZ, expMeanMZ, expCovMZ, dataMZ, objMZ ) modelDZ <- mxModel( "DZ", defAge, pathB, meanDZ, expMeanDZ, expCovDZ, dataDZ, objDZ ) minus2ll <- mxAlgebra( MZ.objective+ DZ.objective, name="minus2sumloglikelihood" ) obj <- mxAlgebraObjective( "minus2sumloglikelihood" ) ciCov <- mxCI( c('MZ.expCovMZ','DZ.expCovDZ' )) ciMean <- mxCI( c('MZ.expMeanMZ','DZ.expMeanDZ' )) twinSatModel <- mxModel( "twinSat", modelMZ, modelDZ, minus2ll, obj, ciCov, ciMean ) # ------------------------------------------------------------------------------ # RUN SATURATED MODELS # Run Saturated Model twinSatFit <- mxRun( twinSatModel, intervals=F ) twinSatSum <- summary( twinSatFit ) twinSatSum # Generate Saturated Model Output twinSatFit$MZ.expMeanMZ@result twinSatFit$DZ.expMeanDZ@result twinSatFit$MZ.expCovMZ@values twinSatFit$DZ.expCovDZ@values # Print Goodness-of-Fit Statistics twinSatSum$observedStatistics length(twinSatSum$parameters[[1]]) twinSatFit@output$Minus2LogLikelihood twinSatSum$degreesOfFreedom twinSatSum$AIC round(twinSatFit@output$estimate,4) # Test Significance of Covariate # Copy model, provide new name testCovModel <- mxModel(twinSatFit, name="testCov") # Change parameter by changing attributes for label testCovModel <- omxSetParameters( testCovModel, label="b11", free=FALSE, values=0 ) # Fit Nested Model testCovFit <- mxRun(testCovModel) # Compare Nested Model with 'Full' Model mxCompare(twinSatFit, testCovFit) # ------------------------------------------------------------------------------ # RUN SUBMODELS # Constrain expected Means to be equal across twin order eqMeansTwinModel <- mxModel(twinSatFit, name="eqMeansTwin" ) eqMeansTwinModel <- omxSetParameters( eqMeansTwinModel, label=c("m1MZ","m2MZ"), free=TRUE, values=svMe, newlabels='mMZ' ) eqMeansTwinModel <- omxSetParameters( eqMeansTwinModel, label=c("m1DZ","m2DZ"), free=TRUE, values=svMe, newlabels='mDZ' ) eqMeansTwinFit <- mxRun( eqMeansTwinModel, intervals=F ) eqMeansTwinSum <- summary( eqMeansTwinFit ) eqMeansTwinLLL <- eqMeansTwinFit@output$Minus2LogLikelihood twinSatLLL <- twinSatFit@output$Minus2LogLikelihood chi2Sat_eqM <- eqMeansTwinLLL-twinSatLLL pSat_eqM <- pchisq( chi2Sat_eqM, lower.tail=F, 2 ) chi2Sat_eqM; pSat_eqM mxCompare(twinSatFit, eqMeansTwinFit) # Constrain expected Means and Variances to be equal across twin order eqMVarsTwinModel <- mxModel(eqMeansTwinFit, name="eqMVarsTwin" ) eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1MZ","v2MZ"), free=TRUE, values=svMe, newlabels='vMZ' ) eqMVarsTwinModel <- omxSetParameters( eqMVarsTwinModel, label=c("v1DZ","v2DZ"), free=TRUE, values=svMe, newlabels='vDZ' ) eqMVarsTwinFit <- mxRun( eqMVarsTwinModel, intervals=F ) subs <- list(eqMeansTwinFit, eqMVarsTwinFit) mxCompare(twinSatFit, subs) # Constrain expected Means and Variances to be equal across twin order and zygosity eqMVarsZygModel <- mxModel(eqMVarsTwinModel, name="eqMVarsZyg" ) eqMVarsZygModel <- omxSetParameters( eqMVarsZygModel, label=c("mMZ","mDZ"), free=TRUE, values=svMe, newlabels='mZ' ) eqMVarsZygModel <- omxSetParameters( eqMVarsZygModel, label=c("vMZ","vDZ"), free=TRUE, values=svMe, newlabels='vZ' ) eqMVarsZygFit <- mxRun( eqMVarsZygModel, intervals=F ) mxCompare(eqMVarsTwinFit, eqMVarsZygFit) # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics SatNested <- list(eqMeansTwinFit, eqMVarsTwinFit, eqMVarsZygFit) mxCompare(twinSatFit, SatNested) # ------------------------------------------------------------------------------ # PREPARE ADE MODEL # Set Starting Values svMe <- 20 # start value for means svPa <- .6 # start value for path coefficients (sqrt(variance/#ofpaths)) # ADE Model # Matrices declared to store a, d, and e Path Coefficients pathA <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", name="a" ) pathD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="d11", name="d" ) pathE <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="e11", name="e" ) # Matrices generated to hold A, D, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covD <- mxAlgebra( expression=d %*% t(d), name="D" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Matrices for covariates and linear regression coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values= .01, label="b11", name="b" ) # Algebra for expected Mean Matrices in MZ & DZ twins meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels="xbmi", name="mean" ) expMean <- mxAlgebra( expression= mean + cbind(b%*%Age,b%*%Age), name="expMean" ) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+D+E, name="V" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, A+D), cbind(A+D, V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+ 0.25%x%D), cbind(0.5%x%A+ 0.25%x%D , V)), name="expCovDZ" ) # 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 ) objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars ) # Combine Groups pars <- list( pathA, pathD, pathE, covA, covD, covE, covP, pathB ) modelMZ <- mxModel( pars, defAge, meanG, expMean, expCovMZ, dataMZ, objMZ, name="MZ" ) modelDZ <- mxModel( pars, defAge, meanG, expMean, expCovDZ, dataDZ, objDZ, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxAlgebraObjective( "m2LL" ) AdeModel <- mxModel( "ADE", pars, modelMZ, modelDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN ADE MODEL # Run ADE model AdeFit <- mxRun(AdeModel) AdeSumm <- summary(AdeFit) AdeSumm mxCompare(twinSatFit,AdeFit) round(AdeFit@output$estimate,4) # Generate Table of Parameter Estimates using mxEval pathEstimatesADE <- print(round(mxEval(cbind(a,d,e), AdeFit),4)) varComponentsADE <- print(round(mxEval(cbind(A/V,D/V,E/V), AdeFit),4)) rownames(pathEstimatesADE) <- 'pathEstimates' colnames(pathEstimatesADE) <- c('a','d','e') rownames(varComponentsADE) <- 'varComponents' colnames(varComponentsADE) <- c('a^2','d^2','e^2') pathEstimatesADE varComponentsADE # Generate ADE Model Output estMean <- mxEval(expMean, AdeFit$MZ) # expected mean estCovMZ <- mxEval(expCovMZ, AdeFit$MZ) # expected covariance matrix for MZ's estCovDZ <- mxEval(expCovDZ, AdeFit$DZ) # expected covariance matrix for DZ's estVA <- mxEval(a*a, AdeFit) # additive genetic variance, a^2 estVD <- mxEval(d*d, AdeFit) # dominance variance, d^2 estVE <- mxEval(e*e, AdeFit) # unique environmental variance, e^2 estVP <- (estVA+estVD+estVE) # total variance estPropVA <- estVA/estVP # standardized additive genetic variance estPropVD <- estVD/estVP # standardized dominance variance estPropVE <- estVE/estVP # standardized unique environmental variance estADE <- rbind(cbind(estVA,estVD,estVE), # table of estimates cbind(estPropVA,estPropVD,estPropVE)) LL_ADE <- mxEval(objective, AdeFit) # likelihood of ADE model # ------------------------------------------------------------------------------ # RERUN MODEL with calculations # Algebras to hold unstandardized and standardized Variance Components rowVars <- rep('vars',nv) colVars <- rep(c('A','D','E','SA','SD','SE'),each=nv) estVars <- mxAlgebra( expression=cbind(A,D,E,A/V,D/V,E/V), name="Vars", dimnames=list(rowVars,colVars)) ciVars <- mxCI("Vars") # Algebra to compute standard deviations (diagonal only) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra to compute standardized Path Coefficients (squared to produce standardizd Variance Components) aS <- mxAlgebra( (iSD %*% a)^2, name="aS" ) dS <- mxAlgebra( (iSD %*% d)^2, name="dS" ) eS <- mxAlgebra( (iSD %*% e)^2, name="eS" ) AdeModel <- mxModel( "ADE", pars, modelMZ, modelDZ, minus2ll, obj, estVars, ciVars, matI, invSD, aS, dS, eS ) AdeFit <- mxRun(AdeModel, intervals=T) AdeSumm <- summary(AdeFit) AdeSumm mxCompare(twinSatFit,AdeFit) round(AdeFit@output$estimate,4) round(AdeFit$Vars@result,4) AdeFit$aS@result # ------------------------------------------------------------------------------ # FIT SUBMODELS # Run AE model AeModel <- mxModel( AdeFit, name="AE" ) AeModel <- omxSetParameters( AeModel, labels="d11", free=FALSE, values=0 ) AeFit <- mxRun(AeModel) mxCompare(AdeFit, AeFit) round(AeFit@output$estimate,4) round(AeFit$Vars@result,4) # Run E model eModel <- mxModel( AeFit, name="E" ) eModel <- omxSetParameters( eModel, labels="a11", free=FALSE, values=0 ) eFit <- mxRun(eModel) mxCompare(AeFit, eFit) round(eFit@output$estimate,4) round(eFit$Vars@result,4) # Print Comparative Fit Statistics AdeNested <- list(AeFit, eFit) mxCompare(AdeFit,AdeNested) round(rbind(AdeFit$Vars@result,AeFit$Vars@result,eFit$Vars@result),4)