checking total h2 at time 2 in a longitudinal choelsky model with paths script
lior abramson
Joined: 07/21/2017
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