Hi,
I have a longitudinal dataset on twins at the age of 6 and 9 years. For the purpose of examining change and stability in genetic and environmental effects, I would like to conduct a Cholesky decomposition analysis, allowing girls and boys to differ in their ACE components and including one categorical covariate in the analysis.
Since I didn't find any example for this kind of analysis with path specification, I decided to adjust a univariate analysis syntax according to my theoretical understanding of the Cholesky model. I would like to consult and ask whether this adjustment is sufficient.
Specifically:
*I estimated the means of the observed variables in both ages, allowing them to differ between age 6 and 9 (different labels)
*For each twin, I specified paths from the age 6 ACE components and covariate to the observed variable of age 6
*For each twin, I specified another set of paths from the age 6 and age 9 ACE and covariates to the observed variable at age 9
*I constrained the C component covariances at age 6 and at age 9 to be equal to 1, but I did not put any constraints on the covariances between the C component across ages (e.g., C of twin 1 at age 6 with C of twin 2 at age 9). I also did not estimate these covariances as free parameters.
*Similarly, I constrained the covariances of the A components according to the twins model assumptions within each age but not across ages.
I attach the script in case it would be useful (and in the case it may help others who would want an example in the future).
Is this the entire Cholesky decomposition model? Or are there any other paths I should specify to make it right?
Thank you very much
# Select Variables for Analysis selVars <- c('POSFAC6_1_Res','POSFAC6_2_Res','LabHome6','POSFAC9_1_Res','POSFAC9_2_Res','LabHome9') 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 mzM <- subset(D.impWide, zygosity==1 & sex_a=='sex1',selVars) mzF <- subset(D.impWide, zygosity==1 & sex_a=='sex2',selVars) dzM <- subset(D.impWide, zygosity==2 & sex_a=='sex1',selVars) dzF <- subset(D.impWide, zygosity==2 & sex_a=='sex2',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 ) varCovar <- mxPath( from= c('LabHome6','LabHome9' ), arrows=2, free=(c(FALSE)), values=1, label='placeVar' ) # 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,F,T,T,F), values=c(StartMean, StartMean, 0,StartMean, StartMean,0), labels=c("meanFAC6","meanFAC6","meanPlace6","meanFAC9","meanFAC9","meanPlace9") ) obsMeansF <- mxPath( from="one", to=selVars, arrows=1, free=c(T,T,F,T,T,F), values=c(StartMean, StartMean, 0,StartMean, StartMean,0), labels=c("meanFAC6","meanFAC6","meanPlace6","meanFAC9","meanFAC9" ,"meanPlace9") ) # path coefficients for twin 1. pathAceTM1.6 <- mxPath( from=c("A1.6","C1.6","E1.6", "LabHome6"), to="POSFAC6_1_Res", arrows=1, free=c(T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c ("aM.6","cM.6" ,"eM.6","place.6") ) pathAceTM1.9 <- mxPath( from=c("A1.6","C1.6","E1.6", "LabHome6","A1.9","C1.9","E1.9", "LabHome9"), to="POSFAC9_1_Res", arrows=1, free=c(T,T,T,T,T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c("aM.6","cM.6" ,"eM.6","place.6","aM.9","cM.9" ,"eM.9","place.9") ) pathAceTF1.6 <- mxPath( from=c("A1.6","C1.6","E1.6", "LabHome6"), to="POSFAC6_1_Res", arrows=1, free=c(T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c("aF.6","cF.6" ,"eF.6","place.6") ) pathAceTF1.9 <- mxPath( from=c("A1.6","C1.6","E1.6", "LabHome6","A1.9","C1.9","E1.9", "LabHome9"), to="POSFAC9_1_Res", arrows=1, free=c(T,T,T,T,T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1,StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c("aF.6","cF.6" ,"eF.6","place.6","aF.9","cF.9" ,"eF.9","place.9") ) # path coefficients for twin 2 pathAceTM2.6 <- mxPath( from=c("A2.6","C2.6","E2.6", "LabHome6"), to="POSFAC6_2_Res", arrows=1, free=c(T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c ("aM.6","cM.6" ,"eM.6","place.6") ) pathAceTM2.9 <- mxPath( from=c("A2.6","C2.6","E2.6", "LabHome6","A2.9","C2.9","E2.9", "LabHome9"), to="POSFAC9_2_Res", arrows=1, free=c(T,T,T,T,T,T,T,T), values= c(StartVar.ACE,StartVar.ACE, StartVar.ACE,1),label=c("aM.6","cM.6" ,"eM.6","place.6","aM.9","cM.9" ,"eM.9","place.9") ) pathAceTF2.6 <- mxPath( from=c("A2.6","C2.6","E2.6", "LabHome6"), to="POSFAC6_2_Res", arrows=1, free=c(T,T,T,T), values= c(StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c("aF.6","cF.6" ,"eF.6","place.6") ) pathAceTF2.9 <- mxPath( from=c("A2.6","C2.6","E2.6", "LabHome6","A2.9","C2.9","E2.9", "LabHome9"), to="POSFAC9_2_Res", arrows=1, free=c(T,T,T,T,T,T,T,T), values= c(StartVar.ACE,StartVar.ACE, StartVar.ACE,1,StartVar.ACE,StartVar.ACE,StartVar.ACE,1),label=c("aF.6","cF.6" ,"eF.6","place.6","aF.9","cF.9","eF.9","place.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 dataMZM <- mxData( observed=mzM, type="raw" ) dataMZF <- mxData( observed=mzF, type="raw" ) dataDZM <- mxData( observed=dzM, type="raw" ) dataDZF <- mxData( observed=dzF, type="raw" ) # Combine Groups #arange together all the paths that are the same for MZ and DZ pathsM <- list(latVariances, latMeans, obsMeansM,pathAceTM1.6,pathAceTM1.9, pathAceTM2.6,pathAceTM2.9, covC1C2.6,covC1C2.9, varCovar ) pathsF <- list(latVariances, latMeans, obsMeansF,pathAceTF1.6,pathAceTF1.9, pathAceTF2.6,pathAceTF2.9, covC1C2.6,covC1C2.9, varCovar ) #Build a model seperately for MZ and DZ and for each sex modelMZM <- mxModel(model="MZM", type="RAM", manifestVars=selVars,latentVars=aceVars, pathsM, threshold, covA1A2_MZ.6, covA1A2_MZ.9, dataMZM ) modelDZM <- mxModel(model="DZM", type="RAM", manifestVars=selVars, latentVars=aceVars, pathsM, threshold, covA1A2_DZ.6,covA1A2_DZ.9, dataDZM ) modelMZF <- mxModel(model="MZF", type="RAM", manifestVars=selVars,latentVars=aceVars, pathsF, threshold, covA1A2_MZ.6 ,covA1A2_MZ.9, dataMZF ) modelDZF <- mxModel(model="DZF", type="RAM", manifestVars=selVars, latentVars=aceVars, pathsF, threshold, covA1A2_DZ.6,covA1A2_DZ.9, dataDZF ) #Ask for confidence intervals CI <-mxCI(reference = c('aM.6','aF.6','cM.6','cF.6','eM.6','eF.6','aM.9','aF.9','cM.9','cF.9','eM.9','eF.9'), interval = 0.95, type=c("both", "lower", "upper")) minus2ll <- mxAlgebra(expression=MZM.fitfunction + DZM.fitfunction+MZF.fitfunction + DZF.fitfunction, name="minus2loglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2loglikelihood" ) modelACE.Cholesky <- mxModel(model="ACE_Cholesky", modelMZM, modelDZM,modelMZF, modelDZF, minus2ll, obj, CI ) fitACE.Cholesky <- mxRun(modelACE.Cholesky, intervals=TRUE) sumACE.Cholesky <- summary(fitACE.Cholesky)
Printer-friendly version
That looks pretty good to me, except for three things.
First, in OpenMx, a period means "member of" (a la C/C++). I don't offhand know if using a period in latent-variable names is going to cause namespace issues for you. I'll assume in the rest of this post that it does not.
Second, replace this,
with this,
. It won't make any difference to your results, but there are some OpenMx features you might want to use later on that are compatible with MxFitFunctionMultigroups but aren't compatble with MxFitFunctionAlgebras.
Third, I don't think your syntax for confidence intervals is going to work, because the MxModel object named "ACE_Cholesky", which contains object 'CI', doesn't contain any of the named quantities you provided in argument
referencetomxCI(). Try one of the following: either create objectand include it in the MxModel object named "ACE_Cholesky", or change the
mxCI()statement to (say)