# ----------------------------------------------------------------------- # Program: MultivariateTwinAnalysis_MatrixRaw.R # Author: Meike Bartels (Based on script H.M.) # Date: 02-03-2010 # # Multivariate Cholesky ACE model to estimate genetic and environmental sources of variance # Matrix style model input - Raw data input # # Data: GFF, Happiness, and Satisfaction with Life, Adolescent Twins, Dutch Health and Behavior Questionnaire # ----------------------------------------------------------------------- require(OpenMx) source("GenEpiHelperFunctions.R") FamData <- read.table ("DHBQ_bs.dat", header=T, na=-999) require(psych) summary (FamData) describe(FamData) names(FamData) Famvars <- colnames (FamData) Vars <- c('gff', 'hap', 'sat') nv <- 3 selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('gff1', 'hap1', 'sat1','gff2' 'hap2', 'sat2') ntv <- nv*2 mzData <- subset(FamData, zyg2==1, selVars) dzData <- subset(FamData, zyg2==2, selVars) # Print Descriptives summary(mzData) summary(dzData) describe(mzData) describe(dzData) 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") # Fit Multivariate ACE Model using Cholesky Decomposition # ----------------------------------------------------------------------- multACEModel <- mxModel("multACE", mxModel("ACE", # Matrices a, c, and e to store a, c, and e path coefficients mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("a11", "a21", "a31", "a22", "a32", "a33"), name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("c11", "c21", "c31", "c22", "c32", "c33"), name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6,label=c("e11", "e21", "e31", "e22", "e32", "e33"), name="e" ), # Matrices A, C, and E compute variance components mxAlgebra( expression=a %*% t(a), name="A" ), mxAlgebra( expression=c %*% t(c), name="C" ), mxAlgebra( expression=e %*% t(e), name="E" ), # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=A+C+E, name="V" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*V)), name="isd"), ## Note that the rest of the mxModel statements do not change for bivariate/multivariate case # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, name="M" ), mxAlgebra( expression= cbind(M,M), name="expMean"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra( expression= rbind ( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ), # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix mxAlgebra( expression= rbind ( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) ), mxModel("MZ", mxData( observed=mzData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovMZ", means="ACE.expMean", dimnames=selVars ) ), mxModel("DZ", mxData( observed=dzData, type="raw" ), mxFIMLObjective( covariance="ACE.expCovDZ", means="ACE.expMean", dimnames=selVars ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="-2sumll" ), mxAlgebraObjective("-2sumll") ) multACEFit <- mxRun(multACEModel) multACESumm <- summary(multACEFit) multACESumm # Generate Multivariate ACE Output # ----------------------------------------------------------------------- parameterSpecifications(multACEFit) expectedMeansCovariances(multACEFit) #tableFitStatistics(multTwinSatFit, multACEFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("ACE.a","ACE.c","ACE.e","ACE.isd","ACE.isd %*% ACE.a","ACE.isd %*% ACE.c","ACE.isd %*% ACE.e") ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","sd","stParEst_a","stPathEst_c","stPathEst_e") formatOutputMatrices(multACEFit,ACEpathMatrices,ACEpathLabels, Vars, 4) ACEcovMatrices <- c("ACE.A","ACE.C","ACE.E","ACE.V","ACE.A/ACE.V","ACE.C/ACE.V","ACE.E/ACE.V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(multACEFit,ACEcovMatrices,ACEcovLabels, Vars, 4) # Genetic and Environmental Correlations ACEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))", "solve(sqrt(ACE.I*ACE.C)) %*% ACE.C %*% solve(sqrt(ACE.I*ACE.C))", "solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))") ACEcorLabels <- c("corComp_A","corComp_C","corComp_E") formatOutputMatrices(multACEFit,ACEcorMatrices,ACEcorLabels, Vars, 4) # Fit Multivariate AE model # ----------------------------------------------------------------------- multAEModel <- mxModel(multACEFit, name="multAE", mxModel(multACEFit$ACE, mxMatrix( type="Lower", nrow=nv, ncol=nv, free=FALSE, values=0, name="c" ) # drop c at 0 ) ) multAEFit <- mxRun(multAEModel) multAESumm <- summary(multAEFit) multAESumm # Generate Multivariate AE Output # ----------------------------------------------------------------------- parameterSpecifications(multAEFit) tableFitStatistics(multACEFit, multAEFit) formatOutputMatrices(multAEFit,ACEpathMatrices,ACEpathLabels,Vars, 4) formatOutputMatrices(multAEFit,ACEcovMatrices,ACEcovLabels,Vars, 4) # Genetic and Environmental Correlations AEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))", "solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))") AEcorLabels <- c("corComp_A","corComp_E") formatOutputMatrices(multAEFit,AEcorMatrices,AEcorLabels, Vars, 4) # Fit Multivariate ACE-AE-AE model # ----------------------------------------------------------------------- multACEAEAEModel <- mxModel(multACEFit, name="multACEAEAE", mxModel(multACEFit$ACE, mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,F,F,F,F,F), values=c(.5,0,0,0,0,0),, name="c" ) ) ) multACEAEAEFit <- mxRun(multACEAEAEModel) multACEAEAESumm <- summary(multACEAEAEFit) multACEAEAESumm # Generate Multivariate ACE-AE-AE Output # ----------------------------------------------------------------------- parameterSpecifications(multACEAEAEFit) tableFitStatistics(multACEFit, multACEAEAEFit) formatOutputMatrices(multACEAEAEFit,ACEpathMatrices,ACEpathLabels,Vars, 4) formatOutputMatrices(multACEAEAEFit,ACEcovMatrices,ACEcovLabels,Vars, 4) # Genetic and Environmental Correlations ACEAEAEcorMatrices <- c("solve(sqrt(ACE.I*ACE.A)) %*% ACE.A %*% solve(sqrt(ACE.I*ACE.A))", "solve(sqrt(ACE.I*ACE.E)) %*% ACE.E %*% solve(sqrt(ACE.I*ACE.E))") ACEAEAEcorLabels <- c("corComp_A","corComp_E") formatOutputMatrices(multACEAEAEFit,ACEAEAEcorMatrices,ACEAEAEcorLabels, Vars, 4) # Print Comparative Fit Statistics # ----------------------------------------------------------------------- multACENested <- list(multAEFit, multACEAEAEFit) tableFitStatistics(multACEFit,multACENested)