# checking total h2 at time 2 in a longitudinal choelsky model with paths script

3 posts / 0 new
Offline
Joined: 07/21/2017 - 13:13
checking total h2 at time 2 in a longitudinal choelsky model with paths script

Hello,
I am conducting a longitudinal Cholesky model with one variable in two time points.
I use the script written below with paths.
How can I modify this script to find for the total h2 at time 2? that is, I want to find the heritability explained by both the genetic effect from time 1 and the genetic effect from time 2, and the confidence interval of this estimate.
Thank you very much for your help

#start values based on all data from both ages (NA=variable at time 1. NA9=variable at time 2)
meanFAC1 <-mean(na.omit(DfsMPLUS$NA_1)) meanFAC2 <-mean(na.omit(DfsMPLUS$NA_2))
meanFAC3 <-mean(na.omit(DfsMPLUS$NA9_1)) meanFAC4 <-mean(na.omit(DfsMPLUS$NA9_2))
StartMean <-mean(c(meanFAC1,meanFAC2,meanFAC3,meanFAC4))

FAC <-c(DfsMPLUS$NA_1,DfsMPLUS$NA_2,DfsMPLUS$NA9_1,DfsMPLUS$NA9_2)
StartVar.ACE <-sqrt(var(na.omit(FAC))/3) #3 because there are 3 ACE components

selVars <- c('NA_1','NA_2','NA9_1','NA9_2')
aceVars <- c("A1_6","C1_6","E1_6","A2_6","C2_6","E2_6","A1_9","C1_9","E1_9"
,"A2_9","C2_9","E2_9")

# Select Data for Analysis
mz <- subset(DfsMPLUS, zygosity==1,selVars)
dz <- subset(DfsMPLUS, zygosity==2,selVars)

# Path objects for Multiple Groups
manifestVars=selVars
latentVars=aceVars

# specify paths for the variances and means of the 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 )

#specify paths for the means of the observed variables.
obsMeansM <- mxPath( from="one", to=selVars, arrows=1,
free=c(T,T,T,T),
values=c(StartMean, StartMean,StartMean, StartMean),
labels=c("meanFAC6" ,"meanFAC6","meanFAC9","meanFAC9" ) )

# path coefficients for twin 1.
pathAceT1.6 <- mxPath( from=c("A1_6","C1_6","E1_6"),to="NA_1", arrows=1, free=c(T,T,T),
values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE),
label=c("a_6","c_6" ,"e_6") )

pathAceT1.9 <- mxPath( from=c("A1_6","C1_6","E1_6","A1_9","C1_9","E1_9"),
to="NA9_1", arrows=1, free=c(T,T,T,T,T,T),
values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,
StartVar.ACE,StartVar.ACE,StartVar.ACE),
label=c("a_69","c_69" ,"e_69","a_9","c_9" ,"e_9") )

# path coefficients for twin 2
pathAceT2.6 <- mxPath( from=c("A2_6","C2_6","E2_6"),to="NA_2", arrows=1, free=c(T,T,T),
values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE),
label=c("a_6","c_6" ,"e_6") )

pathAceT2.9 <- mxPath( from=c("A2_6","C2_6","E2_6","A2_9","C2_9","E2_9"),
to="NA9_2", arrows=1, free=c(T,T,T,T,T,T),
values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,
StartVar.ACE,StartVar.ACE,StartVar.ACE),
label=c("a_69","c_69" ,"e_69","a_9","c_9" ,"e_9") )

# covariance between C1 & C2
covC1C2.6 <- mxPath( from="C1_6", to="C2_6", arrows=2,free=FALSE, values=1 )
covC1C2.9 <- mxPath( from="C1_9", to="C2_9", arrows=2,free=FALSE, values=1 )

# covariance between A1 & A2 in MZ's
covA1A2_MZ.6 <-mxPath( from="A1_6", to="A2_6", arrows=2, free=FALSE,values=1 )
covA1A2_MZ.9 <-mxPath( from="A1_9", to="A2_9", arrows=2, free=FALSE,values=1 )
# covariance between A1 & A2 in DZ's
covA1A2_DZ.6 <-mxPath( from="A1_6", to="A2_6",arrows=2, free=FALSE,values=.5 )
covA1A2_DZ.9 <-mxPath( from="A1_9", to="A2_9",arrows=2, free=FALSE,values=.5 )

#call the data.frame with the MZ raw data, mzData, and the DZ raw data
dataMZ <- mxData(observed=mz, type="raw")

# Combine Groups
#arange together all the paths that are the same for MZ and DZ
paths <- list(latVariances, latMeans, obsMeansM,pathAceT1.6,pathAceT1.9,
pathAceT2.6,pathAceT2.9, covC1C2.6,covC1C2.9 )

#Build a model seperately for MZ and DZ
modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars,
latentVars=aceVars, paths, covA1A2_MZ.6
,covA1A2_MZ.9, dataMZ )
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars,
latentVars=aceVars, paths, covA1A2_DZ.6

CI <-mxCI(reference = c('a_6','c_6','e_6','a_69','c_69','e_69','a_9','c_9','e_9'), interval = 0.95, type=c("both", "lower", "upper"))

#ask for confidence interval for the variance components-age 6
aceMat6 <-mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=StartVar.ACE,labels=c("a_6","c_6","e_6"),name="ace6")
mxal6 <-mxAlgebra( (ace6%^%2) %x% solve(t(ace6)%*%ace6), name="StdVarComp6",dimnames=list(c("a62","c62","e62"),NULL) )

#confidence interval for the variance components
aceMat9 <-mxMatrix(type="Full",nrow=6,ncol=1,free=T,values=StartVar.ACE,labels=c("a_69","c_69","e_69","a_9","c_9","e_9"),name="ace9")
mxal9 <-mxAlgebra( (ace9%^%2) %x% solve(t(ace9)%*%ace9), name="StdVarComp9",dimnames=list(c("a692","c692","e692","a92","c92","e92"),NULL) )

obj <- mxFitFunctionMultigroup(c("MZ","DZ"))
modelACE.Cholesky <- mxModel(model="ACE_Cholesky", modelMZ, modelDZ,obj,CI,
aceMat6,mxal6,mxCI("StdVarComp6"),
aceMat9,mxal9,mxCI("StdVarComp9"))

fitACE.Cholesky <- mxRun(modelACE.Cholesky, intervals=TRUE)
sumACE.Cholesky <- summary(fitACE.Cholesky)
Offline
Joined: 03/01/2013 - 14:09
Data?

Hi

It's difficult to trace the model without some data to run it on. Could you please provide (possibly by generating with mxGenerateData() data of similar structure to the real data that you may not want to share) some data?

Offline
Joined: 03/01/2013 - 11:03
Try the EasyMx way?

I'm not familiar with the way you're setting up this model. In general, you want to compute the variance at time 2 due to additive genetics, then divide this by the total variance from all ACE components at time 2.

If I were doing this, I'd use EasyMx.

require(EasyMx)
require(OpenMx)
data(twinData)
twinVar = names(twinData)
selVars <- c('ht1', 'bmi1','ht2','bmi2')
mzdzData <- subset(twinData, zyg %in% c(1, 3), c(selVars, 'zyg'))
mzdzData$RCoef <- c(1, NA, .5)[mzdzData$zyg]

mod <- emxTwinModel(model='Cholesky', relatedness='RCoef',
data=mzdzData, use=selVars, run=FALSE, name='TwCh')

# Say ht is time 1, and bmi is time 2
# ht1 is twin 1 at time 1
# ht2 is twin 2 and time 1
# etc

modci <- mxModel(mod,
mxAlgebra(name='Time12Herit', diag2vec(A)/diag2vec(A + C + E)),
mxCI('Time12Herit'))
runci <- mxRun(modci, intervals=TRUE)
summary(runci)
# Note in this model the CIs are running into box constraints/bounds

# Total variance due to A at time 2
mxEval(Time12Herit, runci)[2,1]