You are here

bivariate Cholesky descomposition for dummies

2 posts / 0 new
Last post
lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
bivariate Cholesky descomposition for dummies

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)
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
That looks pretty good to me,

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,

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 )

with this,

obj <- mxFitFunctionMultigroup(c("MZM","DZM","MZF","DZF"))
modelACE.Cholesky <- mxModel(model="ACE_Cholesky", modelMZM, modelDZM,modelMZF, modelDZF, obj, CI )

. 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 reference to mxCI(). Try one of the following: either create object

twinPaths <- list(pathAceTM1.6, pathAceTM1.9, pathAceTF1.6, pathAceTF1.9, 
                              pathAceTM2.6, pathAceTM2.9, pathAceTF2.6, pathAceTF2.9)

and include it in the MxModel object named "ACE_Cholesky", or change the mxCI() statement to (say)

CI <-mxCI(reference = c('MZM.aM.6','MZF.aF.6','MZM.cM.6','MZF.cF.6','MZM.eM.6','MZF.eF.6','MZM.aM.9','MZF.aF.9','MZM.cM.9','MZF.cF.9','MZM.eM.9','MZF.eF.9')) #default is to do two-sided 95% intervals
Log in or register to post comments