##input for each of the metabolites data <- read.csv("met1.csv", header=T) require(OpenMx) require(psych) names (data) summary(data) describe(data) str(data) nv <- 1 # number of variables for a twin = 1 in Univariate ntv <- 2*nv # number of variables for a pair = 2* 1 for Univariate Vars <- 'metabolite' selVars <- c('met1','met2') mzData <- subset(data, zyg==1, selVars) dzData <- subset(data, zyg==2, selVars) # Descriptive Statistics describe(data) colMeans(mzData,na.rm=TRUE) cov(mzData,use="complete") colMeans(dzData,na.rm=TRUE) cov(dzData,use="complete") cor(mzData) cor(dzData) sd(mzData$met1, na.rm=T) sd(mzData$met2, na.rm=T) # ------------------------------------------------------------------------------- # 1 Specify and run a fully Saturated Model (Cholesky Decomposition) # ------------------------------------------------------------------------------- # Matrix & Algebra for expected means and covariances CholMZ <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=.04, name="lowMZ" ) CholDZ <-mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=T, values=.04, name="lowDZ" ) MZCov <-mxAlgebra( expression=lowMZ %*% t(lowMZ), name="expCovMZ" ) DZCov <-mxAlgebra( expression=lowDZ %*% t(lowDZ), name="expCovDZ" ) MZMeans <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=.03, labels=c("Mmz1","Mmz2"), name="expMeanMZ" ) DZMeans <-mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=.03, labels=c("Mdz1","Mdz2"), name="expMeanDZ" ) # Algebra's needed for standardizing the covariances matI <-mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I") MZcor <-mxAlgebra( expression= solve(sqrt(I*expCovMZ)) %*% expCovMZ %*% solve(sqrt(I*expCovMZ)), name="rMZ" ) DZcor <-mxAlgebra( expression= solve(sqrt(I*expCovDZ)) %*% expCovDZ %*% solve(sqrt(I*expCovDZ)), name="rDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups # mxExpectationNormal: Objective functions which uses Full–Information maximum # likelihood, the preferred method for raw data. # Objective functions are functions for which free parameter values are chosen # such that the value of the objective function is minimized objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars) fitFunction <- mxFitFunctionML() # Combine Groups modelMZ <- mxModel( CholMZ, MZCov, MZMeans, matI, MZcor, dataMZ, objMZ, fitFunction, name="MZ") modelDZ <- mxModel( CholDZ, DZCov, DZMeans, matI, DZcor, dataDZ, objDZ, fitFunction, name="DZ") minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) Conf <- mxCI (c ('MZ.rMZ[2,1]', 'DZ.rDZ[2,1]') ) SatModel <- mxModel( "Sat", modelMZ, modelDZ, minus2ll, obj, Conf) # ----------------------------------------------------------------------------------- # RUN Saturated Model SatFit <- mxRun(SatModel, intervals=T) (SatSumm <- summary(SatFit)) # Generate some output SatFit$MZ$expCovMZ SatFit$DZ$expCovDZ SatFit$MZ$rMZ SatFit$DZ$rDZ mxEval(MZ.expMeanMZ, SatFit) mxEval(MZ.expCovMZ, SatFit) mxEval(DZ.expMeanDZ, SatFit) mxEval(DZ.expCovDZ, SatFit) # ----------------------------------------------------------------------- # 2 Specify & Run full ACE Model # ----------------------------------------------------------------------- # Matrix & Algebra for expected means vector meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=.03, label="M", name="expMean" ) # Matrices to store a, c, and e path coefficients pathA <-mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.08, label="a11", name="a" ) pathC <-mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.08, label="c11", name="c" ) pathE <-mxMatrix( type="Lower", nrow=nv, ncol=nv, free=T, values=.08, label="e11", name="e" ) # Algebra to generate Matrices 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" ) # Algebra to compute Total Variance covP <-mxAlgebra( expression=A+C+E, name="V" ) # Algebra to get '1/sd' Id <-mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") SD <-mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra to compute standardized path estimates stA <-mxAlgebra( expression=a%*%iSD, name="sta") stC <-mxAlgebra( expression=c%*%iSD, name="stc") stE <-mxAlgebra( expression=e%*%iSD, name="ste") # Algebra to compute standardized variance components h2 <-mxAlgebra( expression=A/V, name="h2") c2 <-mxAlgebra( expression=C/V, name="c2") e2 <-mxAlgebra( expression=E/V, name="e2") # Algebra to generate a Matrix to hold (standardized)estimates rowVars <-Vars colVars <-rep(c('A','C','E','h2','c2','e2'),each=nv) estVars <-mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="Est", dimnames=list(rowVars,colVars)) # Algebra for expected 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" ) # Data objects for Multiple Groups dataMZ <-mxData( observed=mzData, type="raw" ) dataDZ <-mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars) fitFunction <- mxFitFunctionML() # Combine Groups pars <-list(meanG,pathA,pathC,pathE,covA,covC,covE,covP,Id,SD,stA,stC,stE,h2,c2,e2,estVars,fitFunction ) modelMZ <-mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" ) modelDZ <-mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" ) minus2ll <-mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <-mxFitFunctionAlgebra( "m2LL" ) Conf <-mxCI (c ('MZ.h2[1,1]', 'MZ.c2[1,1]', 'MZ.e2[1,1]') ) ACEModel <-mxModel( "ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf) # ------------------------------------------------- # RUN full ACE MODEL ACEFit <-mxRun(ACEModel, intervals=T) (ACESum <-summary(ACEFit)) # Generate ACE Output ACEFit$ACE.h2 ACEFit$ACE.c2 ACEFit$ACE.e2 round(ACEFit@output$estimate,4) round(ACEFit$Est@result,4) # 3 Specify and Run sub-model: AE # ----------------------------------------------------------------------- AEModel <-mxModel( ACEFit, name="AE") AEModel <-omxSetParameters( AEModel, labels="c11", free=F, values=0 ) AEFit <-mxRun(AEModel, intervals=F) (AESum <-summary (AEFit)) round(AEFit@output$estimate,4) round(AEFit$Est@result,4) # 4 Specify and Run sub-model: CE # ----------------------------------------------------------------------- CEModel <-mxModel( ACEFit, name="CE") CEModel <-omxSetParameters( CEModel, labels="a11", free=F, values=0 ) CEFit <-mxRun(CEModel, intervals=F) (CESum <-summary (CEFit)) round(CEFit@output$estimate,4) round(CEFit$Est@result,4) # 5 Specify and Run sub-model: E # ----------------------------------------------------------------------- EModel <- mxModel( ACEFit, name="E") EModel <- omxSetParameters( EModel, labels=c("a11","c11"), free=F, values=0 ) EFit <- mxRun(EModel) (ESum <- summary(EFit)) round(EFit@output$estimate,4) round(EFit$Est@result,4) # COMPARE ALL MODELS: Print Comparative Fit Statistics # --------------------------------------------------------------------------------------------- (Nested.fit <- rbind( mxCompare(SatFit, ACEFit), mxCompare(ACEFit, AEFit)[2,], mxCompare(ACEFit, CEFit)[2,], mxCompare(ACEFit, EFit)[2,] ) ) # Note, the [2,1] within the compare command will suppress the fit of the first model (ACE) to be # printed for that comparison (to avoid repeats in the output table)