# ----------------------------------------------------------------------- # Program: HeterogeneityTwinAnalysis_MatrixRaw.R # Author: Hermine Maes # Date: 12 01 2009 # # Univariate Heterogeneity Twin Analysis model to estimate causes of variation (ADE) # 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) dzfmData <- subset(twinData, zyg==5, selVars) colMeans(mzfData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(mzmData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) colMeans(dzfmData,na.rm=TRUE) cor(mzfData,use="complete") cor(dzfData,use="complete") cor(mzmData,use="complete") cor(dzmData,use="complete") cor(dzfmData,use="complete") # Fit Heterogeneity ADE Model with RawData and Matrices Input # ----------------------------------------------------------------------- univHetADE5Model <- 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=FALSE, values=0, 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=FALSE, values=0, label="df11", name="df" ), mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ef11", name="ef" ), mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=1, label="rg11", lbound=0, ubound=1, name="rg"), # 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"), mxAlgebra( expression= cbind(Mf,Mm), name="expMeanZfm"), # 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" ), mxAlgebra( expression= rbind ( cbind(Af+Df+Ef , 0.5%*%rg%x%(af%*%t(am)) + 0.25%x%(df%*%t(dm))), cbind(0.5%*%rg%x%(am%*%t(af)) + 0.25%x%(dm%*%t(df)) , Am+Dm+Em)), name="expCovDZfm" ) ), 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 ) ), mxModel("DZfm", mxData( observed=dzfmData, type="raw" ), mxFIMLObjective( covariance="ADE.expCovDZfm", means="ADE.expMeanZfm", dimnames=selVars ) ), mxAlgebra( expression=MZm.objective + DZm.objective + MZf.objective + DZf.objective + DZfm.objective, name="-2sumll" ), mxAlgebraObjective("-2sumll") ) univHetADE5Fit <- mxRun(univHetADE5Model) univHetADE5Summ <- summary(univHetADE5Fit) # Generate Heterogeneity ADE Output # ----------------------------------------------------------------------- parameterSpecifications(univHetADE5Fit) expectedMeansCovariances(univHetADE5Fit) tableFitStatistics(univHetADE5Fit) # Generate Table of Parameter Estimates using mxEval pathEstimatesADE <- round(mxEval(cbind(ADE.am,ADE.dm,ADE.em,ADE.af,ADE.df,ADE.ef,ADE.rg), univHetADE5Fit),4) rownames(pathEstimatesADE) <- 'pathEstimates' colnames(pathEstimatesADE) <- c('am','dm','em','af','df','ef','rg') 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), univHetADE5Fit),4) rownames(varComponentsADE) <- 'varComponents' colnames(varComponentsADE) <- c('%am^2','%dm^2','%em^2','%af^2','%df^2','%ef^2') varComponentsADE # Fit Qualitative Sex Differences ACE Model with RawData and Matrices Input # First test for no shared genes between m, f # ----------------------------------------------------------------------- univQualADEModel <- mxModel(univHetADE5Fit, name="univQualADE", mxModel(univHetADE5Fit$ADE, mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="rg11", name="rg") ) ) univQualADEFit <- mxRun(univQualADEModel) univQualADESumm <- summary(univQualADEFit) # Print Comparative Fit Statistics # ----------------------------------------------------------------------- tableFitStatistics(univHetADE5Fit, univQualADEFit) # Fit Qualitative Sex Differences ACE Model with RawData and Matrices Input # Now test for completely shared genes between m, f # ----------------------------------------------------------------------- univQualADEModel <- mxModel(univHetADE5Fit, name="univQualADE", mxModel(univHetADE5Fit$ADE, mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, values=1, label="rg11", name="rg") ) ) univQualADEFit <- mxRun(univQualADEModel) univQualADESumm <- summary(univQualADEFit) # Print Comparative Fit Statistics # ----------------------------------------------------------------------- tableFitStatistics(univHetADE5Fit, univQualADEFit)