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") 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)
Printer-friendly version
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?
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.