# # Copyright 2007-2016 The OpenMx Project # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. # ----------------------------------------------------------------------------- # Program: UnivariateTwinAnalysis_PathRaw.R # Author: Hermine Maes # Date: 2009.08.01 # # ModelType: ACE # DataType: Twin # Field: Human Behavior Genetics # # Purpose: # Univariate Twin Analysis model to estimate causes of variation # Path style model input - Raw data input # # RevisionHistory: # Hermine Maes -- 2009.10.08 updated & reformatted # Ross Gore -- 2011.06.06 added Model, Data & Field metadata # Hermine Maes -- 2014.11.04 piecewise specification # ----------------------------------------------------------------------------- require(OpenMx) require(foreign) require(psych) # Load Library # ----------------------------------------------------------------------------- setwd("D://R2/scripts") getwd() # Load Data data <- read.spss("jiang.sav", to.data.frame = TRUE) # Select Variables for Analysis selVars <- c('y1','y2') aceVars <- c("A1","C1","E1","A2","C2","E2") # Select Data for Analysis mzData <- subset(data, type=="1", selVars) dzData <- subset(data, type=="2", selVars) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Prepare Data # ----------------------------------------------------------------------------- # Path objects for Multiple Groups manifestVars=selVars latentVars=aceVars # variances of latent variables latVariances <- mxPath( from=aceVars, arrows=2, free=FALSE, values=1 ) # means of latent variables latMeans <- mxPath( from="one", to=aceVars, arrows=1, free=FALSE, values=0 ) # means of observed variables obsMeans <- mxPath( from="one", to=selVars, arrows=1, free=TRUE, values=1, labels="mean" ) # path coefficients for twin 1 pathAceT1 <- mxPath( from=c("A1","C1","E1"), to=selVars[1], arrows=1, free=TRUE, values=.5, label=c("a","c","e") ) # path coefficients for twin 2 pathAceT2 <- mxPath( from=c("A2","C2","E2"), to=selVars[2], arrows=1, free=TRUE, values=.5, label=c("a","c","e") ) # covariance between C1 & C2 covC1C2 <- mxPath( from="C1", to="C2", arrows=2, free=FALSE, values=1 ) # covariance between A1 & A2 in MZ twins covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2, free=FALSE, values=1 ) # covariance between A1 & A2 in DZ twins covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2, free=FALSE, values=.5 ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Combine Groups paths <- list( latVariances, latMeans, obsMeans, pathAceT1, pathAceT2, covC1C2 ) modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_MZ, dataMZ ) modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_DZ, dataDZ ) minus2ll <- mxAlgebra( expression=MZ.fitfunction + DZ.fitfunction, name="minus2loglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2loglikelihood" ) twinACEModel <- mxModel(model="ACE", modelMZ, modelDZ, minus2ll, obj ) # Run Model twinACEFit <- mxRun(twinACEModel) twinACESum <- summary(twinACEFit) twinACESum # Fit ACE Model with RawData and Path-Style Input # ----------------------------------------------------------------------------- # Generate & Print Output # mean M <- mxEval(mean, twinACEFit) # additive genetic variance, a^2 A <- mxEval(a*a, twinACEFit) # shared environmental variance, c^2 C <- mxEval(c*c, twinACEFit) # unique environmental variance, e^2 E <- mxEval(e*e, twinACEFit) # total variance V <- (A+C+E) # standardized A a2 <- A/V # standardized C c2 <- C/V # standardized E e2 <- E/V # table of estimates estACE <- rbind(cbind(A,C,E),cbind(a2,c2,e2)) # likelihood of ACE model LL_ACE <- mxEval(fitfunction, twinACEFit) # Get Model Output # ----------------------------------------------------------------------------- # Change Path Objects # path coefficients for twin 1 pathAceT1 <- mxPath( from=c("A1","C1","E1"), to=selVars[1], arrows=1, free=c(T,F,T), values=c(.6,0,.6), label=c("a","c","e") ) # path coefficients for twin 2 pathAceT2 <- mxPath( from=c("A2","C2","E2"), to=selVars[2], arrows=1, free=c(T,F,T), values=c(.6,0,.6), label=c("a","c","e") ) # Combine Groups paths <- list( latVariances, latMeans, obsMeans, pathAceT1, pathAceT2, covC1C2 ) modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_MZ, dataMZ ) modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_DZ, dataDZ ) twinAEModel <- mxModel(model="AE", modelMZ, modelDZ, minus2ll, obj ) # Run Model twinAEFit <- mxRun(twinAEModel) twinAESum <- summary(twinAEFit) # Fit AE model # ----------------------------------------------------------------------------- # Generate & Print Output M <- mxEval(mean, twinAEFit) A <- mxEval(a*a, twinAEFit) C <- mxEval(c*c, twinAEFit) E <- mxEval(e*e, twinAEFit) V <- (A+C+E) a2 <- A/V c2 <- C/V e2 <- E/V estAE <- rbind(cbind(A, C, E),cbind(a2, c2, e2)) LL_AE <- mxEval(fitfunction, twinAEFit) LRT_ACE_AE <- LL_AE - LL_ACE estACE estAE LRT_ACE_AE # Get Model Output # ----------------------------------------------------------------------------- # AE model details twinAESum # CE model pathAceT1 <- mxPath( from=c("A1","C1","E1"), to=selVars[1], arrows=1, free=c(F,T,T), values=c(0,.6,.6), label=c("a","c","e") ) # path coefficients for twin 2 pathAceT2 <- mxPath( from=c("A2","C2","E2"), to=selVars[2], arrows=1, free=c(F,T,T), values=c(0,.6,.6), label=c("a","c","e") ) # Combine Groups paths <- list( latVariances, latMeans, obsMeans, pathAceT1, pathAceT2, covC1C2 ) modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_MZ, dataMZ ) modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars, latentVars=aceVars, paths, covA1A2_DZ, dataDZ ) twinCEModel <- mxModel(model="CE", modelMZ, modelDZ, minus2ll, obj ) twinCEFit <- mxRun(twinCEModel) twinCESum <- summary(twinCEFit) M <- mxEval(mean, twinCEFit) A <- mxEval(a*a, twinCEFit) C <- mxEval(c*c, twinCEFit) E <- mxEval(e*e, twinCEFit) V <- (A+C+E) a2 <- A/V c2 <- C/V e2 <- E/V estCE <- rbind(cbind(A, C, E),cbind(a2, c2, e2)) LL_CE <- mxEval(fitfunction, twinCEFit) LRT_ACE_CE <- LL_CE - LL_ACE estACE estCE twinCESum