# ----------------------------------------------------------------------- # Program: MultivariateTwinAnalysis_MatrixRawCon.R # Author: Hermine Maes # Date: 01 13 2010 # # Multivariate Twin Saturated model to estimate means and (co)variances across multiple groups # Multivariate Cholesky ACE model to estimate genetic and environmental sources of variance # Multivariate Independent Pathway Model and Common Pathway Model # Matrix style model input - Raw data input # # Revision History # Hermine Maes -- 02 01 2010 updated & reformatted # Hermine Maes -- 02 18 2010 added algebras for equality constraints # ----------------------------------------------------------------------- require(OpenMx) require(psych) source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R") # Prepare Data # ----------------------------------------------------------------------- twin.wide$VarR_1<-twin.wide$RAxialLength.1 #Change here for variable twin.wide$VarL_1<-twin.wide$LAxialLength.1 #Change here for variable twin.wide$VarR_2<-twin.wide$RAxialLength.2 #Change here for variable twin.wide$VarL_2<-twin.wide$LAxialLength.2 #Change here for variable Vars <- c('VarR','VarL') nv <- 2 selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="_") ntv <- nv*2 mzData <- subset(twin.wide, Zygosity=='MZ', selVars) dzData <- subset(twin.wide, Zygosity=='DZ', selVars) # Print Descriptive Statistics # ----------------------------------------------------------------------- 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 Saturated Model # ----------------------------------------------------------------------- right <- c(1,3) #select R eye elements of mean vector left <- c(2,4) #select L eye elements of mean vector meanSV <- c(23,23) cholSV <- c( 1,0,0,0, 1,0,0, 1,0, 1) multiTwinSatModel <- mxModel("multiTwinSat", mxModel("MZ", mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=cholSV,name="CholMZ" ), mxAlgebra( expression=CholMZ %*% t(CholMZ), name="expCovMZ" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=meanSV,name="expMeanMZ" ), mxData( observed=mzData, type="raw" ), mxFIMLObjective( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars), # Algebra's needed for equality constraints mxAlgebra( expression=t(diag2vec(expCovMZ)), name="expVarMZ"), mxAlgebra( expression=expVarMZ[1,1:nv], name="expVarMZt1"), mxAlgebra( expression=expVarMZ[1,(nv+1):ntv], name="expVarMZt2"), mxAlgebra( expression=expMeanMZ[1,1:nv], name="expMeanMZt1"), mxAlgebra( expression=expMeanMZ[1,(nv+1):ntv], name="expMeanMZt2"), mxAlgebra( expression=expMeanMZ[1,right], name="expMeanMZR"), #R vs L mxAlgebra( expression=expMeanMZ[1,left], name="expMeanMZL") #R vs L ), mxModel("DZ", mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=cholSV,name="CholDZ" ), mxAlgebra( expression=CholDZ %*% t(CholDZ), name="expCovDZ" ), mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=meanSV,name="expMeanDZ" ), mxData( observed=dzData, type="raw" ), mxFIMLObjective( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars), # Algebra's needed for equality constraints mxAlgebra( expression=t(diag2vec(expCovDZ)), name="expVarDZ"), mxAlgebra( expression=expVarDZ[1,1:nv], name="expVarDZt1"), mxAlgebra( expression=expVarDZ[1,(nv+1):ntv], name="expVarDZt2"), mxAlgebra( expression=expMeanDZ[1,1:nv], name="expMeanDZt1"), mxAlgebra( expression=expMeanDZ[1,(nv+1):ntv], name="expMeanDZt2"), mxAlgebra( expression=expMeanDZ[1,right], name="expMeanDZR"), #R vs L mxAlgebra( expression=expMeanDZ[1,left], name="expMeanDZL") #R vs L ), mxAlgebra( MZ.objective + DZ.objective, name="2sumll" ), mxAlgebraObjective("2sumll") ) multiTwinSatFit <- mxRun(multiTwinSatModel) multiTwinSatSumm <- summary(multiTwinSatFit) multiTwinSatSumm # Generate Saturated Output # ----------------------------------------------------------------------- parameterSpecifications(multiTwinSatFit) expectedMeansCovariances(multiTwinSatFit) tableFitStatistics(multiTwinSatFit) # ?1 # Fit Multivariate Model with Equal Means & Variances across Twin Order,Zygosity and R vs L. Can comment out individual sections to constrain only twin order or R vs L, etc. # ----------------------------------------------------------------------- equateMeansVarsTwinZygModel <- mxModel(multiTwinSatFit, name="equateMeansVarsTwin", # mxConstraint( alg1="MZ.expVarMZt1", "=", alg2="MZ.expVarMZt2", name="VarMZt1_t2"), # mxConstraint( alg1="DZ.expVarDZt1", "=", alg2="DZ.expVarDZt2", name="VarDZt1_t2"), # mxConstraint( alg1="MZ.expMeanMZt1", "=", alg2="MZ.expMeanMZt2", name="MeanMZt1_t2"), # mxConstraint( alg1="DZ.expMeanDZt1", "=", alg2="DZ.expMeanDZt2", name="MeanDZt1_t2"), # mxConstraint( alg1="MZ.expVarMZ", "=", alg2="DZ.expVarDZ", name="VarMZ_DZ"), # mxConstraint( alg1="MZ.expMeanMZ", "=", alg2="DZ.expMeanDZ", name="MeanMZ_DZ"), mxConstraint( alg1="MZ.expMeanMZR", "=", alg2="MZ.expMeanMZL", name="MeanMZR_MZL"), #constrains R vs L mxConstraint( alg1="DZ.expMeanDZR", "=", alg2="DZ.expMeanDZL", name="MeanDZR_DZL") #constrains R vs L ) equateMeansVarsTwinZygFit <- mxRun(equateMeansVarsTwinZygModel) equateMeansVarsTwinZygSumm <- summary(equateMeansVarsTwinZygFit) equateMeansVarsTwinZygSumm tableFitStatistics(multiTwinSatFit, equateMeansVarsTwinZygFit) # Fit Multivariate ACE Model with RawData and Matrices Input # ----------------------------------------------------------------------- meanSVnv <- c(23,23) cholSVnv <- c(.5,0,.5) AFac <- c( "a11","a21", "a22" ) CFac <- c( "c11","c21", "c22" ) EFac <- c( "e11","e21", "e22" ) multiCholACEModel <- mxModel("multiCholACE", 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=cholSVnv, labels=AFac, name="a" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholSVnv, labels=CFac, name="c" ), mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=cholSVnv, labels=EFac, 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" ), # ?2 # Constrain elements of A and E matrices across R and L eyes to be the same mxAlgebra( expression=a[1,1], name="AR"), mxAlgebra( expression=a[2,1], name="AL"), mxAlgebra( expression=e[1,1], name="ER"), mxAlgebra( expression=e[2,1], name="EL"), # 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 bi/multivariate case # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=meanSVnv, name="Mean" ), mxAlgebra( expression= cbind(Mean,Mean), 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 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") ) multiCholACEFit <- mxRun(multiCholACEModel) multiCholACESumm <- summary(multiCholACEFit) multiCholACESumm # Generate Multivariate Cholesky ACE Output # ----------------------------------------------------------------------- parameterSpecifications(multiCholACEFit) expectedMeansCovariances(multiCholACEFit) tableFitStatistics(multiCholACEFit) # 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","stPathEst_a","stPathEst_c","stPathEst_e") formatOutputMatrices(multiCholACEFit,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(multiCholACEFit,ACEcovMatrices,ACEcovLabels,Vars,4) # Fit Multivariate AE Model with RawData and Matrices Input # ----------------------------------------------------------------------- multiCholAEModel <- mxModel(multiCholACEFit, name="multiCholAE", mxModel(multiCholACEModel$ACE, mxMatrix( type="Lower", nrow=nv, ncol=nv, free=FALSE, values=0, name="c" ) ) ) multiCholAEFit <- mxRun(multiCholAEModel) multiCholAESumm <- summary(multiCholAEFit) multiCholAESumm tableFitStatistics(multiCholACEFit,multiCholAEFit) formatOutputMatrices(multiCholAEFit,ACEpathMatrices,ACEpathLabels,Vars,4) AEcovMatrices <- c("ACE.A","ACE.E","ACE.V","ACE.A/ACE.V","ACE.E/ACE.V") AEcovLabels <- c("covComp_A","covComp_E","Var","stCovComp_A","stCovComp_E") formatOutputMatrices(multiCholAEFit,AEcovMatrices,AEcovLabels,Vars,4) # ?3 # Constrain R and L additive and environmental components # ----------------------------------------------------------------------- equateAEModel <- mxModel(multiCholAEModel, name="equateAE", mxConstraint( alg1="ACE.AR", "=", alg2="ACE.AL", name="AR_AL"), mxConstraint( alg1="ACE.ER", "=", alg2="ACE.EL", name="ER_EL") ) equateAEModelFit <- mxRun(equateAEModel) equateAEModelSumm <- summary(equateAEModelFit) equateAEModelSumm tableFitStatistics(multiTwinSatFit, equateAEModel) # ?4 # Stuck here - don't know how to convert constrained path coefficients to variance components # ----------------------------------------------------------------------- AEcovMatrices <- c("ACE.AR_AL","ACE.ER_EL","ACE.V","ACE.AR_AL/ACE.V","ACE.ER_EL/ACE.V") AEcovLabels <- c("covComp_A","covComp_E","Var","stCovComp_A","stCovComp_E") formatOutputMatrices(equateAEModelFit,AEcovMatrices,AEcovLabels,Vars,4)