You are here

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

3 posts / 0 new
Last post
lior abramson's picture
Joined: 07/21/2017 - 13:13
checking total h2 at time 2 in a longitudinal choelsky model with paths script

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))
  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"
# Select Data for Analysis
  mz <- subset(DfsMPLUS, zygosity==1,selVars)
  dz <- subset(DfsMPLUS, zygosity==2,selVars)
# Path objects for Multiple Groups
# 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,
                     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,
                       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,
                       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")
  dataDZ <- mxData(observed=dz, 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
                   ,covA1A2_DZ.9, dataDZ )
#ask for confidence intervals-not standardized
  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,
  fitACE.Cholesky <- mxRun(modelACE.Cholesky, intervals=TRUE)
  sumACE.Cholesky <- summary(fitACE.Cholesky)
AdminNeale's picture
Joined: 03/01/2013 - 14:09


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?

AdminHunter's picture
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.

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)),
runci <- mxRun(modci, intervals=TRUE)
# 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]