# ----------------------------------------------------------------------- # Program: HeterogeneityTwinAnalysis_MatrixRaw.R # Author: Hermine Maes # Date: 12 01 2009 # # Univariate Heterogeneity Twin Analysis model to estimate causes of variation (ACE) # Matrix style model input - Raw data input # # Revision History # Hermine Maes -- 02 24 2010 updated & reformatted # Danielle Dick & Tim York -- 02 24 2010 modified for Mx workshop # ----------------------------------------------------------------------- #make sure GenEpiHelperFunctions in in the proper working directory setwd("C:/Documents and Settings/Danielle Dick/My Documents/R/OpenMxClass/heterogeneity") require(OpenMx) source("GenEpiHelperFunctions.R") # Prepare Data and Print Summary Statistics # ----------------------------------------------------------------------- data(twinData) describe(twinData) Vars <- 'bmi' nv <- 1 selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") ntv <- nv*2 #CREATE MULTIPLE GROUPS mzfData <- subset(twinData, zyg==1, selVars) dzfData <- subset(twinData, zyg==3, selVars) mzmData <- subset(twinData, zyg==2, selVars) dzmData <- subset(twinData, zyg==4, selVars) colMeans(mzfData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(mzmData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) cor(mzfData,use="complete") cor(dzfData,use="complete") cor(mzmData,use="complete") cor(dzmData,use="complete") # Fit Heterogeneity ADE Model with RawData and Matrices Input # ----------------------------------------------------------------------- univHetADEModel <- mxModel("univHetADE", mxModel("ADE", # Matrices a, d, and e to store a, d, and e path coefficients mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="am11", name="am" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="dm11", name="dm" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="em11", name="em" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="af11", name="af" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="df11", name="df" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ef11", name="ef" ), # Matrices A, D, and E compute variance components mxAlgebra( expression=am %*% t(am), name="Am" ), mxAlgebra( expression=dm %*% t(dm), name="Dm" ), mxAlgebra( expression=em %*% t(em), name="Em" ), mxAlgebra( expression=af %*% t(af), name="Af" ), mxAlgebra( expression=df %*% t(df), name="Df" ), mxAlgebra( expression=ef %*% t(ef), name="Ef" ), # Algebra to compute total variances and standard deviations (diagonal only) mxAlgebra( expression=Am+Dm+Em, name="Vm" ), mxAlgebra( expression=Af+Df+Ef, name="Vf" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"), mxAlgebra( expression=solve(sqrt(I*Vm)), name="isdm"), mxAlgebra( expression=solve(sqrt(I*Vf)), name="isdf"), # Matrix & Algebra for expected means vector mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, label="meanm", name="Mm" ), mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 20, label="meanf", name="Mf" ), mxAlgebra( expression= cbind(Mm,Mm), name="expMeanZm"), mxAlgebra( expression= cbind(Mf,Mf), name="expMeanZf"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra( expression= rbind ( cbind(Am+Dm+Em , Am+Dm), cbind(Am+Dm , Am+Dm+Em)), name="expCovMZm" ), mxAlgebra( expression= rbind ( cbind(Af+Df+Ef , Af+Df), cbind(Af+Df , Af+Df+Ef)), name="expCovMZf" ), # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix mxAlgebra( expression= rbind ( cbind(Am+Dm+Em , 0.5%x%Am+0.25%x%Dm), cbind(0.5%x%Am+0.25%x%Dm , Am+Dm+Em)), name="expCovDZm" ), mxAlgebra( expression= rbind ( cbind(Af+Df+Ef , 0.5%x%Af+0.25%x%Df), cbind(0.5%x%Af+0.25%x%Df , Af+Df+Ef)), name="expCovDZf" ) ), mxModel("MZm", mxData( observed=mzmData, type="raw" ), mxFIMLObjective( covariance="ADE.expCovMZm", means="ADE.expMeanZm", dimnames=selVars ) ), mxModel("DZm", mxData( observed=dzmData, type="raw" ), mxFIMLObjective( covariance="ADE.expCovDZm", means="ADE.expMeanZm", dimnames=selVars ) ), mxModel("MZf", mxData( observed=mzfData, type="raw" ), mxFIMLObjective( covariance="ADE.expCovMZf", means="ADE.expMeanZf", dimnames=selVars ) ), mxModel("DZf", mxData( observed=dzfData, type="raw" ), mxFIMLObjective( covariance="ADE.expCovDZf", means="ADE.expMeanZf", dimnames=selVars ) ), mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective, name="-2sumll" ), mxAlgebraObjective("-2sumll") ) univHetADEFit <- mxRun(univHetADEModel) univHetADESumm <- summary(univHetADEFit) # Generate Heterogeneity ADE Output # ----------------------------------------------------------------------- parameterSpecifications(univHetADEFit) expectedMeansCovariances(univHetADEFit) tableFitStatistics(univHetADEFit) # Generate Table of Parameter Estimates using mxEval pathEstimatesADE <- round(mxEval(cbind(ADE.am,ADE.dm,ADE.em,ADE.af,ADE.df,ADE.ef), univHetADEFit),4) rownames(pathEstimatesADE) <- 'pathEstimates' colnames(pathEstimatesADE) <- c('am','dm','em','af','df','ef') pathEstimatesADE varComponentsADE <- round(mxEval(cbind(ADE.Am/ADE.Vm,ADE.Dm/ADE.Vm,ADE.Em/ADE.Vm,ADE.Af/ADE.Vf,ADE.Df/ADE.Vf,ADE.Ef/ADE.Vf), univHetADEFit),4) rownames(varComponentsADE) <- 'varComponents' colnames(varComponentsADE) <- c('%am^2','%dm^2','%em^2','%af^2','%df^2','%ef^2') varComponentsADE # Fit Homogeneity ADE Model with RawData and Matrices Input # ----------------------------------------------------------------------- univHomADEModel <- mxModel(univHetADEFit, name="univHomADE", mxModel(univHetADEFit$ADE, mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="am" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="d11", name="dm" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="em" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="a11", name="af" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="d11", name="df" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="e11", name="ef" ) ) ) univHomADEFit <- mxRun(univHomADEModel) univHomADESumm <- summary(univHomADEFit) # Generate Homogeneity ADE Output # ----------------------------------------------------------------------- parameterSpecifications(univHomADEFit) expectedMeansCovariances(univHomADEFit) tableFitStatistics(univHomADEFit) # Generate Table of Parameter Estimates using mxEval pathEstimatesADE <- round(mxEval(cbind(ADE.am,ADE.dm,ADE.em), univHomADEFit),4) varComponentsADE <- round(mxEval(cbind(ADE.Am/ADE.Vm,ADE.Dm/ADE.Vm,ADE.Em/ADE.Vm), univHomADEFit),4) rownames(pathEstimatesADE) <- 'pathEstimates' colnames(pathEstimatesADE) <- c('a','d','e') rownames(varComponentsADE) <- 'varComponents' colnames(varComponentsADE) <- c('%a^2','%d^2','%e^2') pathEstimatesADE varComponentsADE # Print Comparative Fit Statistics # ----------------------------------------------------------------------- tableFitStatistics(univHetADEFit, univHomADEFit)