# ----------------------------------------------------------------------------------------------------------- # Script: 2ACE_Chol.R # Fits a MV Cholesky ACE Decomposition # Data File: MVasb.dat (E-Risk Study) # Phenotypes: mothers-, teachers-, examiners-, and children's self-report on antisocial behavior # Type of data: Raw continuous data # ZYG: 1=MZ 2=DZ # Assumptions: means and variances can be equated across birth order & zygosity groups # ------------------------------------------------------------------------------------------------------------ require(OpenMx) require(psych) # Read in data file W5Data <- read.csv("W5OpenMXF.csv") #header=TRUE, na.strings='-99') W5S <- subset(W5Data, randomTwin == 1) #tells R to only select cases in which the variable = 1 source("GenEpiHelperFunctions.R") summary(W5S) describe(W5S) nv <- 3 # number of variables ntv <- 2*nv # number of variables*max family size Vars <- c('dfswanhim','dfswanin','LogSCT') # This yields a list of the 4 variables for twin 1 and twin 2: # Mother1 Teacher1 Examiner1 Self1 Mother2 Teacher2 Examiner2 Self2 selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") mzData <- subset(W5S, zyg==1, selVars) dzData <- subset(W5S, zyg==2, selVars) # To create start values Stmean <-colMeans(mzData[,1:3],na.rm=TRUE) # To create Labels for Lower Triangular Matrices aLabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") cLabs <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") eLabs <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") # To create Labels for Column and Diagonal Matrices mLabs <- paste("m",1:nv,sep="") # ------------------------------------------------------------------------------------------------- # (model 2) Specify Cholesky ACE Decomposition # ------------------------------------------------------------------------------------------------- # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.4,.4,.4,.4,.4,.4), labels=c("a11", "a21", "a31", "a22", "a32", "a33"), name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.1,.1,.1,.1,.1,.1), labels=c("c11", "c21", "c31", "c22", "c32", "c33"), name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c(T,T,T,T,T,T), values=c(.2,.2,.2,.2,.2,.2), labels=c("e11", "e21", "e31", "e22", "e32", "e33"), 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" ) covP <- mxAlgebra( expression=A+C+E, name="V" ) StA <- mxAlgebra( expression=A/V, name="h2" ) StC <- mxAlgebra( expression=C/V, name="c2" ) StE <- mxAlgebra( expression=E/V, name="e2" ) # Algebra to compute Correlations matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") Rph <- mxAlgebra( expression=solve(sqrt(I*V)) %&% V , name="Phcor") Rg <- mxAlgebra( expression=solve(sqrt(I*A)) %&% A , name="Acor") Rc <- mxAlgebra( expression=solve(sqrt(I*C)) %&% C , name="Ccor") Re <- mxAlgebra( expression=solve(sqrt(I*E)) %&% E , name="Ecor") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins Mean <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=T, values=Stmean, labels=c(mLabs,mLabs), name="ExpMean" ) 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( pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, Rph, Rg, Rc, Re) modelMZ <- mxModel( pars, covMZ, Mean, dataMZ, objMZ, fitFunction, name="MZ" ) modelDZ <- mxModel( pars, covDZ, Mean, dataDZ, objDZ, fitFunction, name="DZ" ) minus2ll<- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) conf1 <- mxCI (c ('MZ.h2[1,1]', 'MZ.h2[2,2]', 'MZ.h2[3,3]', 'MZ.h2[4,4]') ) # h2 conf2 <- mxCI (c ('MZ.c2[1,1]', 'MZ.c2[2,2]', 'MZ.c2[3,3]', 'MZ.c2[4,4]') ) # c2 conf3 <- mxCI (c ('MZ.e2[1,1]', 'MZ.e2[2,2]', 'MZ.e2[3,3]', 'MZ.e2[4,4]') ) # e2 conf4 <- mxCI (c ('MZ.Acor[2,1]','MZ.Acor[3,1]','MZ.Acor[4,1]','MZ.Acor[3,2]','MZ.Acor[4,2]','MZ.Acor[4,3]'))#Rg conf5 <- mxCI (c ('MZ.Ccor[2,1]','MZ.Ccor[3,1]','MZ.Ccor[4,1]','MZ.Ccor[3,2]','MZ.Ccor[4,2]','MZ.Ccor[4,3]'))#Rc conf6 <- mxCI (c ('MZ.Ecor[2,1]','MZ.Ecor[3,1]','MZ.Ecor[4,1]','MZ.Ecor[3,2]','MZ.Ecor[4,2]','MZ.Ecor[4,3]'))#Re CholAceModel<- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj, conf1, conf2, conf3, conf4, conf5, conf6) #--------------------------------------------------------------------------------------------------------------------- # (model 2) RUN Cholesky Decomposition ACE Model CholAceFit <- mxRun(CholAceModel, intervals=F) (CholAceSum <- summary(CholAceFit)) #round(CholAceFit@output\$estimate,4) # Print the h2, c2, e2, rG, rC, rE matrices (note: in the summary output you will find them with 95% CI) mxEval(CholACE.h2, CholAceFit) mxEval(CholACE.c2, CholAceFit) mxEval(CholACE.e2, CholAceFit) mxEval(CholACE.Acor, CholAceFit) mxEval(CholACE.Ccor, CholAceFit) mxEval(CholACE.Ecor, CholAceFit) mxEval(CholACE.Phcor, CholAceFit)