# ----------------------------------------------------------------------- # Program: 07_Analysis_Independent1ACE_DATE.R # ----------------------------------------------------------------------- # Datasets # Source: SXtwinDat.txt (converted to wide from SXlong from MOAFTSData from IG) # DataType: Raw, Twin Data # ----------------------------------------------------------------------- # Tasks: # 1) Load Packages, Data # 2) Prepare Data # 3) Build ACE Model # 4) Run ACE Model, Output Fit, Parameter Estimates # 5) Build/Run 2ACE Model (from ACE model), Compare Fit # 6) Run 2ACE Model, Output Fit, Parameter Estimates # ----------------------------------------------------------------------- # Revision History # Hermine Maes 2014 03 05 (from Boulder Workshop, for OpenMx 1.4) # Jarrod Ellingson 2014 07 17 (adapated for current study) # ----------------------------------------------------------------------- # Comments # # ----------------------------------------------------------------------- # References # citation("OpenMx") # http://ibg.colorado.edu/cdrom2014/maes/Multivariate/Multivariate.R # ----------------------------------------------------------------------- # ----------------------------------------------------------------------- # 1) Load Packages & Functions # ----------------------------------------------------------------------- # A) Packages (set working directory uncomment to install OpenMx, psych) # ----------------------------------------------------------------------------- #source('http://openmx.psyc.virginia.edu/getOpenMx.R') require(OpenMx) require(psych) # B) Functions (from Bartels [link above] -- set working directory to folder with .R files # ----------------------------------------------------------------------------- source("myFunctions.R") source("GenEpiHelperFunctions.R") # ----------------------------------------------------------------------------- # 2) Prepare Data # ----------------------------------------------------------------------------- # A) Load Data # ----------------------------------------------------------------------------- twinData <- read.table(file="SXtwinR.txt", header=TRUE, na=-999) Vars <- c("SX1_s", "SX2_s", "SX3_s", "SX4_s", "SX8_s", "SX9_s", "SX10_s", "SX14_s", "SX152_s") nv <- length(Vars) # number of variables ntv <- nv*2 # number of total variables selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="_") # B) Select Data for Analysis, Convert data to numeric matrix # ----------------------------------------------------------------------------- mzData <- as.matrix(subset(twinData, twinzyg==1, selVars)) dzData <- as.matrix(subset(twinData, twinzyg==3, selVars)) MZrows <- nrow(mzData) DZrows <- nrow(dzData) NumMZdata <- matrix(as.numeric(unlist(mzData)),nrow=nrow(mzData)) NumDZdata <- matrix(as.numeric(unlist(dzData)),nrow=nrow(dzData)) colnames(NumMZdata) <- selVars colnames(NumDZdata) <- selVars class(NumMZdata[, 1]) class(NumDZdata[, 1]) # C) Descriptive Statistics # ----------------------------------------------------------------------------- # describe(mzData, skew=F) # describe(dzData, skew=F) # cov(NumMZdata,use="complete") # cov(NumDZdata,use="complete") # cor(NumMZdata,use="complete") # cor(NumDZdata,use="complete") # ----------------------------------------------------------------------------- # 3) Build ACE Model # ----------------------------------------------------------------------------- # A) Set Paramater Start Values, Labels # ----------------------------------------------------------------------------- svMe <- c(rep(0,nv)) # start value for means svPa <- valDiag(nv,.6) # start values for parameters on diagonal lbPa <- valLUDiag(nv,.0001,-10,NA) # lower bounds for parameters on diagonal nf <- 1 # number of factors # B) Matrices ac, cc, and ec to store a, c, and e path coefficients for common factors # ----------------------------------------------------------------------------- pathAc1 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("ac1",nv,nf), name="ac1" ) pathAc2 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=FALSE, values=0, labels=labFull("ac2",nv,nf), name="ac2" ) # set to 0, only modeling 1 common A factor in this model pathCc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("cc",nv,nf), name="cc" ) pathEc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=labFull("ec",nv,nf), name="ec" ) # C) Matrices as, cs, and es to store a, c, and e path coefficients for specific factors # ----------------------------------------------------------------------------- pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=labDiag("as",nv), lbound=.00001, name="as" ) pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=labDiag("cs",nv), lbound=.00001, name="cs" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=labDiag("es",nv), lbound=.00001, name="es" ) pathRac <- mxMatrix (type="Diag", nrow=1, ncol=1, free=FALSE, values=0, label="rac1", lbound=0, ubound=1, name="rac") # D) Matrices A, C, and E compute variance components # ----------------------------------------------------------------------------- covAc1 <- mxAlgebra( expression=ac1 %*% t(ac1), name="Ac1" ) covAc2 <- mxAlgebra( expression=ac2 %*% t(ac2), name="Ac2" ) covAc <- mxAlgebra( expression=Ac1 + Ac2, name="Ac" ) covCc <- mxAlgebra( expression=cc %*% t(cc), name="Cc" ) covEc <- mxAlgebra( expression=ec %*% t(ec), name="Ec" ) covAs <- mxAlgebra( expression=as %*% t(as), name="As" ) covCs <- mxAlgebra( expression=cs %*% t(cs), name="Cs" ) covEs <- mxAlgebra( expression=es %*% t(es), name="Es" ) covA <- mxAlgebra( expression=Ac + As, name="A" ) covC <- mxAlgebra( expression=Cc + Cs, name="C" ) covE <- mxAlgebra( expression=Ec + Es, name="E" ) # E) Algebra to compute total variances and standard deviations (diagonal only) # ----------------------------------------------------------------------------- covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # F) Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins # ----------------------------------------------------------------------------- meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labFull("me",1,nv), name="expMean" ) covMZ <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" ) # G) Data objects for Multiple Groups # ----------------------------------------------------------------------------- dataMZ <- mxData( observed=NumMZdata, type="raw" ) dataDZ <- mxData( observed=NumDZdata, type="raw" ) # H) FIML objects for Multiple Groups # ----------------------------------------------------------------------------- objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars ) objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars ) # I) Combine components, for each group (MZ, DZ), into full model # ----------------------------------------------------------------------------- pars <- list( pathAc1, pathAc2, pathCc, pathEc, pathAs, pathCs, pathEs, pathRac, covAc1, covAc2, covAc, covCc, covEc, covAs, covCs, covEs, covA, covC, covE, covP, matI, invSD, meanG ) 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 <- mxAlgebraObjective( "m2LL" ) IndAceModel <- mxModel( "IndACE", pars, modelMZ, modelDZ, minus2ll, obj ) # ----------------------------------------------------------------------------- # 4) Run ACE Model, Output Fit, Parameter Estimates # ----------------------------------------------------------------------------- # A) Fit Cholesky Decomposition ACE model # ----------------------------------------------------------------------------- IndAceFit <- mxRun(IndAceModel) summary(IndAceFit) # ----------------------------------------------------------------------------- # 5) Build/Run 2ACE Model (from ACE model), Compare Fit # ----------------------------------------------------------------------------- # A) Set Paramater Start Values, Labels # ----------------------------------------------------------------------------- # Create Free and Values for 2 Additive Genetic Factors # free values # A1 A2 A1 A2 # SX1_s F T 0 .6 # SX2_s F T 0 .6 # SX3_s F T 0 .6 # SX4_s F T 0 .6 # SX8_s F T 0 .6 # SX9_s F T 0 .6 # SX10_s T F .6 0 # SX14_s T F .6 0 # SX152_s T F .6 0 Ac1Free <- c(rep(T,6),rep(F,3)) Ac2Free <- c(rep(F,6),rep(T,3)) Ac1Values <- c(rep(1,6),rep(0,3)) Ac2Values <- c(rep(0,6),rep(1,3)) # B) Create 2ACE model with same structure as ACE # ----------------------------------------------------------------------------- Ind2A1C1EeModel <- mxModel(IndAceFit, name="Ind2Ace") # C) Replace Old matrix named 'ac' by Newly Created Object 'pathAc' # Needs to be done in every MxModel in which it appears (MZ, DZ and overall model) # ----------------------------------------------------------------------------- pathAc1 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac1Free, values=Ac1Values, labels=labFull("ac1",nv,nf), name="ac1" ) pathAc2 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=Ac2Free, values=Ac2Values, labels=labFull("ac2",nv,nf), name="ac2" ) pathRac <- mxMatrix (type="Diag", nrow=1, ncol=1, free=TRUE, values=0.5, label="rac1", lbound=0, ubound=1, name="rac") # didn't work #covAc <- mxAlgebra( expression=ac1 %x% rac %x% t(ac2) + ac2 %x% rac %x% t(ac1), name="Ac" ) # didn't work #covAc <- mxAlgebra( expression=ac1 %x% t(ac2) + ac2 %x% t(ac1), name="Ac" ) # didn't work #covAc <- mxAlgebra( expression=Ac1 %x% rac %x% t(Ac2), name="Ac" ) # didn't work covAc <- mxAlgebra( expression=ac1 %*% rac %*% t(ac2) + ac2 %*% rac %*% t(ac1), name="Ac" ) # covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E, A+C), # cbind(A+C, A+C+E)), name="expCovMZ" ) # don't need to change? # covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E, 0.5%x%A+C), # cbind(0.5%x%A+C, A+C+E)), name="expCovDZ") # don't need to change? Ind2A1C1EeModel$MZ$ac1 <- pathAc1 Ind2A1C1EeModel$DZ$ac1 <- pathAc1 Ind2A1C1EeModel$ac1 <- pathAc1 Ind2A1C1EeModel$MZ$ac2 <- pathAc2 Ind2A1C1EeModel$DZ$ac2 <- pathAc2 Ind2A1C1EeModel$ac2 <- pathAc2 #Ind2A1C1EeModel$MZ$rac <- pathRac #Ind2A1C1EeModel$DZ$rac <- pathRac #Ind2A1C1EeModel$rac <- pathRac Ind2A1C1EeModel$MZ$Ac <- covAc Ind2A1C1EeModel$DZ$Ac <- covAc Ind2A1C1EeModel$Ac <- covAc # D) Run 2ACE model # ----------------------------------------------------------------------------- Ind2A1C1EFit <- mxRun(Ind2A1C1EeModel) summary(Ind2A1C1EFit) mxCompare(Ind2A1C1EFit, IndAceFit)