320 likes | 459 Views
Moderating the Covariance between Family Member’s Behavior. Brad Verhulst Sarah Medland. Motivation. Substance Use of Family Members is highly correlated This may be a function of: Environmental Factors Special Twin Environment Sibling Interaction Vertical Cultural Transmission
E N D
Moderating the Covariance between Family Member’s Behavior Brad Verhulst Sarah Medland
Motivation Substance Use of Family Members is highly correlated This may be a function of: • Environmental Factors • Special Twin Environment • Sibling Interaction • Vertical Cultural Transmission • Genetic Factors • Additive & Non-Additive But twins may be a larger part of their co-twins environment not because of some magical special twin environment, but because they are the same age and are simply more likely to have the same interests, share the same friends, and have the same life experiences But if genes are expressed at different points in the lifespan, the differential expression will look like genetic non-additivity.
So What’s the Alternative? • The correlation between behavior decays may decay as a function of differences along a continuum. • In this case we will focus on age difference but realistically several possible moderating variable could be used. • The larger the age difference between relatives the lower the correlation • Two reasons: • As siblings are more different in age, they are less likely to share environments • Age specific gene expression
Background Methodology Both Genetic and Environmental factors have cumulative effects Some Effects are specific to the time of measurement Related: Hopper and Matthews (1982, 1983) - The longer siblings live together, the more similar they are.
Relationship between the negative exponential function and a linear time series with discrete time points
Existing Gene-Environment Interaction Models Purcell (2002) MZ=1 DZ = ½ Classical Twin Design 1 1 1 A A 1 1 C C a + βaM a + βaM 1 1 E E c c + βcM + βcM e + βeM e + βeM Basic Means and Variances Pt1 Pt2 Means Moderation Model μ + Mβm μ + Mβm 1
Our Model MZ= e-|Δage| *αa DZ = e-|Δage| * ½*αa 1 1 A A e-|Δage| * αc 1 1 C C a a 1 1 T T c c 1 1 E t E t e e Pt1 Pt2 μ + Mβm μ + Mβm 1 Cov = Acov* e-|Δage|*αa + Ccov* e-|Δage|*αc + Tcov
Estimated Variance Components for the A, C, T and E Parameters Note: Red Line Indicates the Simulated Value
The Relationship between the Simulated and Estimated Moderation Parameters
What does our model predict about the correlations between relatives? Exercise: • Using ModGraph.R • Run the function • Adjust the parameter values to fit your expectations about the decay in the correlations • How does age related gene expression affect the model? • What would happen if siblings shared less of a common environment as they were more discrepant in age? ModGraph(a=1,c=1,t=1,e=1, AgeRange = 10, aM = 1, cM = 1)
Data Requirements • 5 Types of Families (with twins) • Sex limitation model as the baseline model • Each Family has up to 6 people: • 2 twins • 2 brothers • 2 sisters • Missing data on the phenotypes is allowed • Code missing phenotypic observations NA • Missing data on the definition variables is NOT allowed • Code missing definition variables some arbitrarily large negative number: eg -9999
Reading the Data into R # ------------------------------------------------------------------------------ # Program: TaSMod.R # Author: Brad Verhulst # Date: 10 22, 2012 # # Twin and Sibling Model with Genetic and Environmental Moderation Parameters # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------| # Call required packages require(OpenMx) # ------------------------------------------------------------------------------ # PREPARE DATA setwd("C:/Users/Brad Verhulst/Desktop/NIDA mod Scripts") # Read in the Data TaSdata<-read.csv("TaSdata.csv")
Basic Model Parameters #Basic Parameters Required for the model nv <- 1 Vars <- c('tw','tw','br','br','si','si') selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") selAges <- paste('age',1:6, sep="") ntv <- nv*2 npeople <-6 # Subset the data MZMdata <- subset(TaSdata, zyg==1,c(selVars,selAges)) MZFdata <- subset(TaSdata, zyg==2,c(selVars,selAges)) DZMdata <- subset(TaSdata, zyg==3,c(selVars,selAges)) DZFdata <- subset(TaSdata, zyg==4,c(selVars,selAges)) DZOdata <- subset(TaSdata, zyg==5,c(selVars,selAges))
General Model Parameters ### Specifying the Parameters for the model # Matrices to store the Path Coefficients am <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="am11", lbound = .001, name="am") cm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cm11", lbound = .001, name="cm" ) tm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="tm11", lbound = .001, name="tm" ) em <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="em11", lbound = .001, name="em" ) af <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="af11", lbound = .001, name="af" ) cf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="cf11", lbound = .001, name="cf" ) tf <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="tf11", lbound = .001, name="tf" ) ef <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, label="ef11", lbound = .001, name="ef" ) rg <- mxMatrix( type="Diag", nrow=1, ncol=1, free=FALSE, values=1, label="rg11", lbound = -1, ubound = 1, name="rg" ) paths <- c(am, cm, tm ,em, af, cf, tf, ef, rg) # Matrices A, C, and E compute variance components Am <- mxAlgebra( expression=am %*% t(am), name="Am" ) Cm <- mxAlgebra( expression=cm %*% t(cm), name="Cm" ) Tm <- mxAlgebra( expression=tm %*% t(tm), name="Tm" ) Em <- mxAlgebra( expression=em %*% t(em), name="Em" ) Af <- mxAlgebra( expression=af %*% t(af), name="Af" ) Cf <- mxAlgebra( expression=cf %*% t(cf), name="Cf" ) Tf <- mxAlgebra( expression=tf %*% t(tf), name="Tf" ) Ef <- mxAlgebra( expression=ef %*% t(ef), name="Ef" ) VCs <- c(Am, Cm, Tm, Em, Af, Cf, Tf, Ef, Amf, Cmf, Tmf, Emf, Afm, Cfm, Tfm, Efm)
Specifying the Genetic Correlation Matrix mzmA <- mxAlgebra(expression= rbind(cbind(Am , Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(Am , Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, .5*Af , Af )), name = "mzmA") mzfA <- mxAlgebra(expression= rbind(cbind(Af , Af , .5*Amf, .5*Amf, .5*Af , .5*Af ), cbind(Af , Af , .5*Amf, .5*Amf, .5*Af , .5*Af ), cbind(.5*Afm, .5*Afm, Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Afm, .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Af , .5*Af , .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Af, .5*Af , .5*Afm, .5*Afm, .5*Af , Af )), name = "mzfA") dzmA <- mxAlgebra(expression= rbind(cbind(Am , .5*Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, .5*Af , Af )), name = "dzmA") dzfA <- mxAlgebra(expression= rbind(cbind(Af , .5*Af , .5*Amf, .5*Amf, .5*Af , .5*Af ), cbind(.5*Af , Af , .5*Amf, .5*Amf, .5*Af , .5*Af ), cbind(.5*Afm, .5*Afm, Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Afm, .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Af , .5*Af , .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Af , .5*Af , .5*Afm, .5*Afm, .5*Af , Af )), name = "dzfA") dzoA <- mxAlgebra(expression= rbind(cbind(Am , .5*Amf, .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Afm, Af , .5*Amf, .5*Amf, .5*Af , .5*Af ), cbind(.5*Am , .5*Afm, Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Afm, .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Af , .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Afm, .5*Af , .5*Afm, .5*Afm, .5*Af , Af )), name = "dzoA") mzmA <- mxAlgebra(expression= rbind(cbind(Am , Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(Am , Am , .5*Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , Am , .5*Am , .5*Amf, .5*Amf), cbind(.5*Am , .5*Am , .5*Am , Am , .5*Amf, .5*Amf), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, Af , .5*Af ), cbind(.5*Afm, .5*Afm, .5*Afm, .5*Afm, .5*Af , Af )), name = "mzmA")
Specifying the Common Environmental Correlation Matrix malC <- mxAlgebra(expression= rbind(cbind(Cm , Cm , Cm , Cm , Cmf, Cmf), cbind(Cm , Cm , Cm , Cm , Cmf, Cmf), cbind(Cm , Cm , Cm , Cm , Cmf, Cmf), cbind(Cm , Cm , Cm , Cm , Cmf, Cmf), cbind(Cfm , Cfm, Cfm, Cfm, Cf , Cf ), cbind(Cfm , Cfm, Cfm, Cfm, Cf , Cf )),name = "malC") femC <- mxAlgebra(expression= rbind(cbind(Cf , Cf , Cmf, Cmf, Cf , Cf ), cbind(Cf , Cf , Cmf, Cmf, Cf , Cf ), cbind(Cfm , Cfm, Cm , Cm , Cmf, Cmf), cbind(Cfm , Cfm, Cm , Cm , Cmf, Cmf), cbind(Cf , Cf , Cfm, Cfm, Cf , Cf ), cbind(Cf , Cf , Cfm, Cfm, Cf , Cf )),name = "femC") oppC <- mxAlgebra(expression= rbind(cbind(Cm , Cmf, Cm , Cm , Cmf, Cmf), cbind(Cfm , Cf , Cmf, Cmf, Cf , Cf ), cbind(Cm , Cfm, Cm , Cm , Cmf, Cmf), cbind(Cm , Cfm, Cm , Cm , Cmf, Cmf), cbind(Cfm , Cf , Cfm, Cfm, Cf , Cf ), cbind(Cfm , Cf , Cfm, Cfm, Cf , Cf )),name = "oppC")
Specifying the Special Twin Environment Correlation Matrix malT <- mxAlgebra(expression= rbind(cbind(Tm, Tm, 0 , 0 , 0 , 0 ), cbind(Tm, Tm, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Tm, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Tm, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Tf, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Tf)),name = "malT") femT <- mxAlgebra(expression= rbind(cbind(Tf, Tf, 0 , 0 , 0 , 0 ), cbind(Tf, Tf, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Tm, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Tm, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Tf, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Tf)),name = "femT") oppT <- mxAlgebra(expression= rbind(cbind(Tm, Tmf, 0 , 0 , 0 , 0 ), cbind(Tfm, Tf, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Tm, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Tm, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Tf, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Tf)),name = "oppT")
Specifying the Unique Environment Correlation Matrix malE <- mxAlgebra(expression= rbind(cbind(Em, 0 , 0 , 0 , 0 , 0 ), cbind(0 , Em, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Em, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Em, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Ef, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Ef)),name = "malE") femE <- mxAlgebra(expression= rbind(cbind(Ef, 0, 0 , 0 , 0 , 0 ), cbind(0 , Ef, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Em, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Em, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Ef, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Ef)),name = "femE") oppE <- mxAlgebra(expression= rbind(cbind(Em, 0 , 0 , 0 , 0 , 0 ), cbind(0 , Ef, 0 , 0 , 0 , 0 ), cbind(0 , 0 , Em, 0 , 0 , 0 ), cbind(0 , 0 , 0 , Em, 0 , 0 ), cbind(0 , 0 , 0 , 0 , Ef, 0 ), cbind(0 , 0 , 0 , 0 , 0 , Ef)),name = "oppE")
These are the interesting matrices ## Variance - Covariance Matrices (Unmoderated) MaleVar <- mxAlgebra( expression=Am+Cm+Tm+Em, name="MaleVar" ) ## Male Total Variance FemaVar <- mxAlgebra( expression=Af+Cf+Tf+Ef, name="FemaVar" ) ### Female Total Variance ## Means (Unmoderated) Meanm <- mxMatrix( type="Full", nrow=1, ncol=nv, free=T, values= 1, label="meanm", name="Meanm" ) Meanf <- mxMatrix( type="Full", nrow=1, ncol=nv, free=T, values= 1, label="meanf", name="Meanf" ) ## Matrices to store the effect of Age on the Mean bm <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=0, label=c("intM"), name="bm" ) bf <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=0, label=c("intF"), name="bf" ) # Matrices to store the moderation parameters aMOD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.2, label="aM11", name="aMOD" ) cMOD <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0 , label="cM11", name="cMOD" ) unit <- mxMatrix(type="Full",nrow=nv, ncol=npeople, free=FALSE, values=1, name="unit")
Matrices for the definition terms ### Matrices for the definition terms # MZM mzmDtw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="Dmt1") #twin1 mzmDtw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="Dmt2") #twin2 mzmDms1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age3"), name="Dms1") #Male Sib 1 mzmDms2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age4"), name="Dms2") #Male Sib 2 mzmDfs1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age5"), name="Dfs1") #Female Sib 1 mzmDfs2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age6"), name="Dfs2") #Female Sib 2 # MZF mzfDtw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="Dft1") #twin1 mzfDtw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="Dft2") #twin2 mzfDms1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age3"), name="Dms1") #Male Sib 1 mzfDms2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age4"), name="Dms2") #Male Sib 2 mzfDfs1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age5"), name="Dfs1") #Female Sib 1 mzfDfs2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age6"), name="Dfs2") #Female Sib 2 # DZM dzmDtw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="Dtw1") #twin1 dzmDtw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="Dtw2") #twin2 dzmDms1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age3"), name="Dms1") #Male Sib 1 dzmDms2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age4"), name="Dms2") #Male Sib 2 dzmDfs1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age5"), name="Dfs1") #Female Sib 1 dzmDfs2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age6"), name="Dfs2") #Female Sib 2 # DZF dzfDtw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="Dtw1") #twin1 dzfDtw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="Dtw2") #twin2 dzfDms1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age3"), name="Dms1") #Male Sib 1 dzfDms2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age4"), name="Dms2") #Male Sib 2 dzfDfs1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age5"), name="Dfs1") #Female Sib 1 dzfDfs2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age6"), name="Dfs2") #Female Sib 2 # DZO dzoDtw1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age1"), name="Dtw1") #twin1 dzoDtw2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age2"), name="Dtw2") #twin2 dzoDms1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age3"), name="Dms1") #Male Sib 1 dzoDms2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age4"), name="Dms2") #Male Sib 2 dzoDfs1 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age5"), name="Dfs1") #Female Sib 1 dzoDfs2 <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age6"), name="Dfs2") #Female Sib 2
Algebra for the definition terms ### Algebra for effect of the definition term on the mean # MZM mzmdefmt1 <- mxAlgebra( expression= bm %*% Dmt1, name="defmt1") mzmdefmt2 <- mxAlgebra( expression= bm %*% Dmt2, name="defmt2") mzmdefms1 <- mxAlgebra( expression= bm %*% Dms1, name="defms1") mzmdefms2 <- mxAlgebra( expression= bm %*% Dms2, name="defms2") mzmdeffs1 <- mxAlgebra( expression= bf %*% Dfs1, name="deffs1") mzmdeffs2 <- mxAlgebra( expression= bf %*% Dfs2, name="deffs2") # MZF mzfdefmt1 <- mxAlgebra( expression= bf %*% Dft1, name="defft1") mzfdefmt2 <- mxAlgebra( expression= bf %*% Dft2, name="defft2") mzfdefms1 <- mxAlgebra( expression= bm %*% Dms1, name="defms1") mzfdefms2 <- mxAlgebra( expression= bm %*% Dms2, name="defms2") mzfdeffs1 <- mxAlgebra( expression= bf %*% Dfs1, name="deffs1") mzfdeffs2 <- mxAlgebra( expression= bf %*% Dfs2, name="deffs2") # DZM dzmdefmt1 <- mxAlgebra( expression= bm %*% Dtw1, name="defmt1") dzmdefmt2 <- mxAlgebra( expression= bm %*% Dtw2, name="defmt2") dzmdefms1 <- mxAlgebra( expression= bm %*% Dms1, name="defms1") dzmdefms2 <- mxAlgebra( expression= bm %*% Dms2, name="defms2") dzmdeffs1 <- mxAlgebra( expression= bf %*% Dfs1, name="deffs1") dzmdeffs2 <- mxAlgebra( expression= bf %*% Dfs2, name="deffs2") # DZF dzfdefmt1 <- mxAlgebra( expression= bf %*% Dtw1, name="defft1") dzfdefmt2 <- mxAlgebra( expression= bf %*% Dtw2, name="defft2") dzfdefms1 <- mxAlgebra( expression= bm %*% Dms1, name="defms1") dzfdefms2 <- mxAlgebra( expression= bm %*% Dms2, name="defms2") dzfdeffs1 <- mxAlgebra( expression= bf %*% Dfs1, name="deffs1") dzfdeffs2 <- mxAlgebra( expression= bf %*% Dfs2, name="deffs2") # DZO dzodefmt1 <- mxAlgebra( expression= bm %*% Dtw1, name="defmt1") dzodefmt2 <- mxAlgebra( expression= bf %*% Dtw2, name="defft2") dzodefms1 <- mxAlgebra( expression= bm %*% Dms1, name="defms1") dzodefms2 <- mxAlgebra( expression= bm %*% Dms2, name="defms2") dzodeffs1 <- mxAlgebra( expression= bf %*% Dfs1, name="deffs1") dzodeffs2 <- mxAlgebra( expression= bf %*% Dfs2, name="deffs2")
Calculating the Mean ## Algebra to calculate the Expected Mean expMeanMZm <- mxAlgebra( expression= cbind((Meanm + defmt1),(Meanm + defmt2),(Meanm + defms1),(Meanm + defms2),(Meanf + deffs1),(Meanf + deffs2)), name="expMeanMZm") expMeanMZf <- mxAlgebra( expression= cbind((Meanf + defft1),(Meanf + defft2),(Meanm + defms1),(Meanm + defms2),(Meanf + deffs1),(Meanf + deffs2)), name="expMeanMZf") expMeanDZm <- mxAlgebra( expression= cbind((Meanm + defmt1),(Meanm + defmt2),(Meanm + defms1),(Meanm + defms2),(Meanf + deffs1),(Meanf + deffs2)), name="expMeanDZm") expMeanDZf <- mxAlgebra( expression= cbind((Meanf + defft1),(Meanf + defft2),(Meanm + defms1),(Meanm + defms2),(Meanf + deffs1),(Meanf + deffs2)), name="expMeanDZf") expMeanDZmf <- mxAlgebra( expression= cbind((Meanm + defmt1),(Meanf + defft2),(Meanm + defms1),(Meanm + defms2),(Meanf + deffs1),(Meanf + deffs2)), name="expMeanDZmf")
Setting up the Moderation Matrix ### Setting up the moderation matrix # Data for Ages AgeMZm <- mxMatrix(type="Full",nrow=npeople, ncol=nv, free=FALSE, labels=c("data.age1","data.age2","data.age3","data.age4","data.age5","data.age6"), name="AgeMZm") AgeMZf <- mxMatrix(type="Full",nrow=npeople, ncol=nv, free=FALSE, labels=c("data.age1","data.age2","data.age3","data.age4","data.age5","data.age6"), name="AgeMZf") AgeDZm <- mxMatrix(type="Full",nrow=npeople, ncol=nv, free=FALSE, labels=c("data.age1","data.age2","data.age3","data.age4","data.age5","data.age6"), name="AgeDZm") AgeDZf <- mxMatrix(type="Full",nrow=npeople, ncol=nv, free=FALSE, labels=c("data.age1","data.age2","data.age3","data.age4","data.age5","data.age6"), name="AgeDZf") AgeDZmf <- mxMatrix(type="Full",nrow=npeople, ncol=nv, free=FALSE, labels=c("data.age1","data.age2","data.age3","data.age4","data.age5","data.age6"), name="AgeDZmf") # Age matrix (pre) ageMatMZm <- mxAlgebra(AgeMZm %x% unit, name = "ageMatMZm") ageMatMZf <- mxAlgebra(AgeMZf %x% unit, name = "ageMatMZf") ageMatDZm <- mxAlgebra(AgeDZm %x% unit, name = "ageMatDZm") ageMatDZf <- mxAlgebra(AgeDZf %x% unit, name = "ageMatDZf") ageMatDZmf <- mxAlgebra(AgeDZmf %x% unit, name = "ageMatDZmf") # Age Difference Matrix deltaMZm <- mxAlgebra(abs(ageMatMZm - t(ageMatMZm)), name = "deltaMZm" ) deltaMZf <- mxAlgebra(abs(ageMatMZf - t(ageMatMZf)), name = "deltaMZf" ) deltaDZm <- mxAlgebra(abs(ageMatDZm - t(ageMatDZm)), name = "deltaDZm" ) deltaDZf <- mxAlgebra(abs(ageMatDZf - t(ageMatDZf)), name = "deltaDZf" ) deltaDZmf <- mxAlgebra(abs(ageMatDZmf - t(ageMatDZmf)), name = "deltaDZmf" ) %x%
Algebras for the Moderation Matrices # Algebra for the Genetic Moderation Parameter aMODmzm <- mxAlgebra(exp(-deltaMZm %x% aMOD), name = "aMODmzm") aMODmzf <- mxAlgebra(exp(-deltaMZf %x% aMOD), name = "aMODmzf") aMODdzm <- mxAlgebra(exp(-deltaDZm %x% aMOD), name = "aMODdzm") aMODdzf <- mxAlgebra(exp(-deltaDZf %x% aMOD), name = "aMODdzf") aMODdzmf <- mxAlgebra(exp(-deltaDZmf %x% aMOD), name = "aMODdzmf") # Algebra for the Common Environmental Moderation Parameter cMODmzm <- mxAlgebra(exp(-deltaMZm %x% cMOD), name = "cMODmzm") cMODmzf <- mxAlgebra(exp(-deltaMZf %x% cMOD), name = "cMODmzf") cMODdzm <- mxAlgebra(exp(-deltaDZm %x% cMOD), name = "cMODdzm") cMODdzf <- mxAlgebra(exp(-deltaDZf %x% cMOD), name = "cMODdzf") cMODdzmf <- mxAlgebra(exp(-deltaDZmf %x% cMOD), name = "cMODdzmf")
Specifying the Expected Covariance Matrices ## Expected Covariance Matrices expCovMZm <- mxAlgebra(mzmA * aMODmzm + malC* cMODmzm + malT + malE, name="expCovMZm") expCovMZf <- mxAlgebra(mzfA * aMODmzf + femC* cMODmzf + femT + femE , name="expCovMZf") expCovDZm <- mxAlgebra(dzmA * aMODdzm + malC* cMODdzm + malT + malE, name="expCovDZm") expCovDZf <- mxAlgebra(dzfA * aMODdzf + femC* cMODdzf + femT + femE, name="expCovDZf") expCovDZmf <- mxAlgebra(dzoA * aMODdzmf+ oppC* cMODdzmf+ oppT + oppE, name="expCovDZmf")
Specifying the Data and Objective Functions # Specifying the Data Frames dataMZM <- mxData( observed=MZMdata, type="raw" ) dataMZF <- mxData( observed=MZFdata, type="raw") dataDZM <- mxData( observed=DZMdata, type="raw" ) dataDZF <- mxData( observed=DZFdata, type="raw" ) dataDZO <- mxData( observed=DZOdata, type="raw" ) # Specifying the Objective Functions MZMobj <- mxFIMLObjective( covariance="expCovMZm", means="expMeanMZm", dimnames=selVars ) MZFobj <- mxFIMLObjective( covariance="expCovMZf", means="expMeanMZf", dimnames=selVars ) DZMobj <- mxFIMLObjective( covariance="expCovDZm", means="expMeanDZm", dimnames=selVars ) DZFobj <- mxFIMLObjective( covariance="expCovDZf", means="expMeanDZf", dimnames=selVars ) DZOobj <- mxFIMLObjective( covariance="expCovDZmf", means="expMeanDZmf", dimnames=selVars )
Specifying the Submodels MZm <- mxModel("MZm", paths, VCs, aMOD,cMOD, bm,bf, MaleVar,FemaVar, mzmA,aMODmzm,malC,cMODmzm,malT,malE, mzmDtw1,mzmDtw2,mzmDms1,mzmDms2,mzmDfs1,mzmDfs2, mzmdefmt1,mzmdefmt2,mzmdefms1,mzmdefms2,mzmdeffs1,mzmdeffs2, Meanm,Meanf,expMeanMZm, unit,AgeMZm,ageMatMZm,deltaMZm, aMODmzm,cMODmzm,expCovMZm, dataMZM,MZMobj ) MZf <- mxModel("MZf", paths, VCs, aMOD,cMOD, bm,bf, MaleVar,FemaVar, mzfA,aMODmzf,femC,cMODmzf,femT,femE, mzfDtw1,mzfDtw2,mzfDms1,mzfDms2,mzfDfs1,mzfDfs2, mzfdefmt1,mzfdefmt2,mzfdefms1,mzfdefms2,mzfdeffs1,mzfdeffs2, Meanm,Meanf,expMeanMZf, unit,AgeMZf,ageMatMZf,deltaMZf, aMODmzf,cMODmzf,expCovMZf, dataMZF,MZFobj ) DZm <- mxModel("DZm", paths, VCs, aMOD,cMOD, bm,bf, MaleVar,FemaVar, dzmA,aMODdzm,malC,cMODdzm,malT,malE, dzmDtw1,dzmDtw2,dzmDms1,dzmDms2,dzmDfs1,dzmDfs2, dzmdefmt1,dzmdefmt2,dzmdefms1,dzmdefms2,dzmdeffs1,dzmdeffs2, Meanm,Meanf,expMeanDZm, unit,AgeDZm,ageMatDZm,deltaDZm, aMODdzm,cMODdzm,expCovDZm, dataDZM,DZMobj ) DZf <- mxModel("DZf", paths, VCs, aMOD,cMOD, bm,bf, MaleVar,FemaVar, dzfA,aMODdzf,femC,cMODdzf,femT,femE, dzfDtw1,dzfDtw2,dzfDms1,dzfDms2,dzfDfs1,dzfDfs2, dzfdefmt1,dzfdefmt2,dzfdefms1,dzfdefms2,dzfdeffs1,dzfdeffs2, Meanm,Meanf,expMeanDZf, unit,AgeDZf,ageMatDZf,deltaDZf, aMODdzf,cMODdzf,expCovDZf, dataDZF,DZFobj ) DZmf <- mxModel("DZmf", paths, VCs, aMOD,cMOD, bm,bf, MaleVar,FemaVar, dzoA,aMODdzmf,oppC,cMODdzmf,oppT,oppE, dzoDtw1,dzoDtw2,dzoDms1,dzoDms2,dzoDfs1,dzoDfs2, dzodefmt1,dzodefmt2,dzodefms1,dzodefms2,dzodeffs1,dzodeffs2, Meanm,Meanf,expMeanDZmf, unit,AgeDZmf,ageMatDZmf,deltaDZmf, aMODdzmf,cMODdzmf,expCovDZmf, dataDZO,DZOobj )
Running the Model TaSModModel <- mxModel("TaSHetACE", MZm, MZf, DZm, DZf, DZmf, mxAlgebra( expression=MZm.objective + MZf.objective +DZm.objective + DZf.objective + DZmf.objective, name="minus2sumloglikelihood" ), mxAlgebraObjective("minus2sumloglikelihood") ) TaSModFit <- mxRun(TaSModModel ) TaSModSumm <- summary(TaSModFit ) TaSModSumm
Looking at the Output free parameters: name matrix row col Estimate Std.Errorlboundubound 1 am11 MZm.am 1 1 0.69356237 0.020621014 0.001 2 cm11 MZm.cm 1 1 0.45026574 0.017106106 0.001 3 tm11 MZm.tm 1 1 0.27076106 0.035208626 0.001 4 em11 MZm.em 1 1 0.49112887 0.010067007 0.001 5 af11 MZm.af 1 1 0.44093396 0.025880749 0.001 6 cf11 MZm.cf 1 1 0.70389297 0.014464905 0.001 7 tf11 MZm.tf 1 1 0.03834467 0.061631613 0.001 8 ef11 MZm.ef 1 1 0.50301198 0.009813318 0.001 9 aM11 MZm.aMOD 1 1 0.35600616 0.684446181 10 cM11 MZm.cMOD 1 1 0.31349403 0.253406803 11 intM MZm.bm 1 1 -0.03724208 0.156987583 12 intF MZm.bf 1 1 0.14033247 0.106558294 13 meanmMZm.Meanm 1 1 0.02990390 0.084356156 14 meanfMZm.Meanf 1 1 -0.07092350 0.059893376 observed statistics: 30000 estimated parameters: 14 degrees of freedom: 29986 -2 log likelihood: 71950.61 saturated -2 log likelihood: NA number of observations: 5000 chi-square: NA p: NA Information Criteria: df Penalty Parameters Penalty Sample-Size Adjusted AIC: 11978.61 71978.61 NA BIC: -183445.95 72069.85 72025.36 CFI: NA TLI: NA RMSEA: NA timestamp: 2012-10-18 11:36:44 frontend time: 3.249828 secs backend time: 1.830538 mins independent submodels time: 0 secs wall clock time: 113.0821 secs cpu time: 113.0821 secs openmx version number: 1.3.0-2169
Looking at the Output -2 log likelihood: 71950.61 saturated -2 log likelihood: NA number of observations: 5000 chi-square: NA p: NA Information Criteria: df Penalty Parameters Penalty Sample-Size Adjusted AIC: 11978.61 71978.61 NA BIC: -183445.95 72069.85 72025.36 CFI: NA TLI: NA RMSEA: NA timestamp: 2012-10-18 11:36:44 frontend time: 3.249828 secs backend time: 1.830538 mins independent submodels time: 0 secs wall clock time: 113.0821 secs cpu time: 113.0821 secs openmx version number: 1.3.0-2169
Communicating the Results Males Females
SpecialAcknowledgements Thank You R25DA026119 (PI: Neale) Mike Neale Lindon Eaves HermineMaes