# 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))

```
``` 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")

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,

aceMat6,mxal6,mxCI("StdVarComp6"),

aceMat9,mxal9,mxCI("StdVarComp9"))

` fitACE.Cholesky <- mxRun(modelACE.Cholesky, intervals=TRUE)`

sumACE.Cholesky <- summary(fitACE.Cholesky)

## Data?

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?

Log in or register to post comments

## Try the EasyMx way?

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]

Log in or register to post comments