# ------------------------------------------------------------------------------ # Program: twinAceOrd.R # Author: Hermine Maes # Date: 02 20 2012 # # Univariate Twin Analysis model to estimate causes of variation # Matrix style model - Raw data - Ordinal data # -------|---------|---------|---------|---------|---------|---------|---------| # Load Library require(OpenMx) require(psych) # ------------------------------------------------------------------------------ # PREPARE DATA # Load Data data(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') nth <- 4 # number of thresholds # Select Data for Analysis mzData <- subset(twinData, zyg==1, selVars) dzData <- subset(twinData, zyg==3, selVars) mzDataOrd <- data.frame(mzData) for (i in 1:2) { mzDataOrd[,i] <- cut(mzData[,i], breaks=c(18,20.5,21,21.5,22,28), ordered_result=T, labels=c(0,1,2,3,4)) } dzDataOrd <- data.frame(dzData) for (i in 1:2) { dzDataOrd[,i] <- cut(dzData[,i], breaks=c(18,20.5,21,21.5,22,28), ordered_result=T, labels=c(0,1,2,3,4)) } mzDataOrdF <- mxFactor( x=mzDataOrd, levels=c(0:4) ) dzDataOrdF <- mxFactor( x=dzDataOrd, levels=c(0:4) ) # Set Starting Values lth <- -1.5 # start value for first threshold ith <- 1 # start value for increments thVal <- matrix(rep(c(lth,(rep(ith,nth-1)))),nrow=nth,ncol=nv) thLB <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) thUB <- 3 corVals <- .5 # start value for correlations lbrVal <- -0.99 # start value for lower bounds ubrVal <- 0.99 # start value for upper bounds thLab <- paste("t",1:nth,"Z",sep="") # Generate Descriptive Statistics #colMeans(mzDataOrd,na.rm=TRUE) #colMeans(dzDataOrd,na.rm=TRUE) cov(mzDataOrd,use="complete") cov(dzDataOrd,use="complete") table(mzDataOrd$bmi1,mzDataOrd$bmi2) table(dzDataOrd$bmi1,dzDataOrd$bmi2) # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="a11", name="a" ) pathC <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="c11", name="c" ) pathE <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.6, label="e11", name="e" ) # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Matrix & Algebra for expected means vector and expected thresholds meanG <- mxMatrix( type="Zero", nrow=1, ncol=nv, name="Mean" ) meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ) threG <- mxMatrix( type="Full", nrow=nth, ncol=nv, free=TRUE, values=thVal, lbound=thLB, ubound=thUB, labels=thLab, name="Thre" ) Inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="Inc" ) threT <- mxAlgebra( expression= cbind(Inc %*% Thre,Inc %*% Thre), name="expThre" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) # Algebras generated to hold Parameter Estimates and Derived Variance Components rowVars <- rep('vars',nv) colVars <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Vars", dimnames=list(rowVars,colVars)) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) # Constraint on variance of Ordinal variables matUnv <- mxMatrix( type="Unit", nrow=nv, ncol=1, name="Unv1" ) var1 <- mxConstraint( expression=diag2vec(V)==Unv1, name="Var1" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzDataOrdF, type="raw" ) dataDZ <- mxData( observed=dzDataOrdF, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="expThre" ) objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="expThre" ) # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, Inc, covP, matUnv, estVars ) modelMZ <- mxModel( pars, meanG, meanT, threG, threT, covMZ, dataMZ, objMZ, name="MZ" ) modelDZ <- mxModel( pars, meanG, meanT, threG, threT, covDZ, dataDZ, objDZ, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxAlgebraObjective( "m2LL" ) AceModel <- mxModel( "ACE", pars, var1, modelMZ, modelDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN MODELS # Run ACE Model AceFit <- mxRun(AceModel, unsafe=T) AceSumm <- summary(AceFit) AceSumm round(AceFit@output$estimate,4) round(AceFit$Vars@result,4) # Generate ACE Model Output source("GenEpiHelperFunctions.R") parameterSpecifications(AceFit) expectedMeansCovariances(AceFit) tableFitStatistics(AceFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatrices <- c("a","c",path"e","iSD","iSD %*% a","iSD %*% c","iSD %*% e") ACEpathLabels <- c("path_a","path_c","path_e","isd","spath_a","spath_c","spath_e") formatOutputMatrices(AceFit,ACEpathMatrices,ACEpathLabels,Vars,4) ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covs_A","covs_C","covs_E","Var","prop_A","prop_C","prop_E") formatOutputMatrices(AceFit,ACEcovMatrices,ACEcovLabels,Vars,4) # ------------------------------------------------------------------------------ # FIT SUBMODELS # Run AE model AeModel <- mxModel( AceFit, name="AE" ) AeModel <- omxSetParameters( AeModel, labels="c11", free=FALSE, values=0 ) AeFit <- mxRun(AeModel) round(AeFit@output$estimate,4) round(AeFit$Vars@result,4) # ------------------------------------------------------------------------------ # Run CE model CeModel <- mxModel( AceFit, name="CE" ) CeModel <- omxSetParameters( CeModel, labels="a11", free=FALSE, values=0 ) CeFit <- mxRun(CeModel) round(CeFit@output$estimate,4) round(CeFit$Vars@result,4) # ------------------------------------------------------------------------------ # Run E model eModel <- mxModel( AceFit, name="E" ) eModel <- omxSetParameters( eModel, labels="a11", free=FALSE, values=0 ) eFit <- mxRun(eModel) round(eFit@output$estimate,4) round(eFit$Vars@result,4) # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics AceNested <- list(AeFit, CeFit, eFit) tableFitStatistics(AceFit,AceNested) round(rbind(AceFit@output$estimate,AeFit@output$estimate,CeFit@output$estimate,eFit@output$estimate),4) round(rbind(AceFit$Vars@result,AeFit$Vars@result,CeFit$Vars@result,eFit$Vars@result),4)