You are here

Bivariate Cholesky model with paths specification

4 posts / 0 new
Last post
lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
Bivariate Cholesky model with paths specification

Dear Forum,
I wrote a syntax for bivariate cholesky model for two variables measuring two different traits at the same time-point. However, I think that I may have missed something. In the output I see the ACE coefficients for each observed variable, but I do not see the the ACE coefficients to the covariance between the two traits (i.e. the coefficients that explain what part of the variance is explained by shared genetic and environmental influences on the two traits).

Do you know what I should add to the syntax in order to receive this output?

I attach the syntax and the output I received.

Thank you so much,
Lior

Syntax:

#start values
meanFAC1 <-mean(na.omit(D6$POSFAC6_1))
meanFAC2 <-mean(na.omit(D6$POSFAC6_2))
StartMean <-mean(c(meanFAC1,meanFAC2))
 
FAC <-c(D6$POSFAC6_1,D6$POSFAC6_2) 
StartVar.ACE <-sqrt(var(na.omit(FAC))/3)
 
#Do the same for the AE and CE model- dividing the variance by 2 instead of 3
StartVar.AE <-sqrt(var(na.omit(FAC))/2)
StartVar.E <-sqrt(var(na.omit(FAC))/1)
 
require("OpenMx")
 
# Select Variables for Analysis
selVars <- c('POSFAC6_1','POSFAC6_2','ObsPInitia_a','ObsPInitia_b')
aceVars <- c("A1_Pos","C1_Pos","E1_Pos","A2_Pos","C2_Pos","E2_Pos","A1_Pro","C1_Pro","E1_Pro"
             ,"A2_Pro","C2_Pro","E2_Pro")
 
# Select Data for Analysis
mz <- subset(D6, zygosity==1,selVars)
dz <- subset(D6, 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.
obsMeans <- mxPath( from="one", to=selVars, arrows=1,
                     free=c(T,T,T,T), 
                     values=c(StartMean, StartMean,StartMean, StartMean), 
                     labels=c("meanFAC_Pos" ,"meanFAC_Pos","meanFAC_Pro","meanFAC_Pro"
                              ) )
 
 
# path coefficients for twin 1. 
pathAceT1_Pos <- mxPath( from=c("A1_Pos","C1_Pos","E1_Pos"), 
                         to="POSFAC6_1", arrows=1, free=c(T,T,T), 
                         values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE),
                         label=c("aPos","cPos" ,"ePos"))
 
 
 
pathAceT1_Pro <- mxPath( from=c("A1_Pos","C1_Pos","E1_Pos","A1_Pro","C1_Pro","E1_Pro"), 
                         to="ObsPInitia_a", 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("aPos","cPos" ,"ePos","aPro","cPro" ,"ePro") )
 
 
# path coefficients for twin 2
pathAceT2_Pos <- mxPath( from=c("A2_Pos","C2_Pos","E2_Pos"), 
                         to="POSFAC6_2", arrows=1, free=c(T,T,T), 
                         values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE),
                         label=c("aPos","cPos" ,"ePos"))
 
pathAceT2_Pro <- mxPath( from=c("A2_Pos","C2_Pos","E2_Pos","A2_Pro","C2_Pro","E2_Pro"), 
                         to="ObsPInitia_b", 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("aPos","cPos" ,"ePos","aPro","cPro" ,"ePro") )
 
# covariance between C1 & C2
covC1C2.Pos <- mxPath( from="C1_Pos", to="C2_Pos", arrows=2,free=FALSE, values=1 )
covC1C2.Pro <- mxPath( from="C1_Pro", to="C2_Pro", arrows=2,free=FALSE, values=1 )
 
 
# covariance between A1 & A2 in MZ's
covA1A2_MZ_Pos <-mxPath( from="A1_Pos", to="A2_Pos", arrows=2, free=FALSE,values=1 )
covA1A2_MZ_Pro <-mxPath( from="A1_Pro", to="A2_Pro", arrows=2, free=FALSE,values=1 )
# covariance between A1 & A2 in DZ's
covA1A2_DZ_Pos <-mxPath( from="A1_Pos", to="A2_Pos",arrows=2, free=FALSE,values=.5 )
covA1A2_DZ_Pro <-mxPath( from="A1_Pro", to="A2_Pro",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, obsMeans,pathAceT1_Pos,
              pathAceT2_Pos,pathAceT1_Pro,pathAceT2_Pro,covC1C2.Pos,covC1C2.Pro)
 
 
#Build a model seperately for MZ and DZ and for each sex
modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars,
                    latentVars=aceVars, paths, covA1A2_MZ_Pos
                    ,covA1A2_MZ_Pro, dataMZ)
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars,
                    latentVars=aceVars, paths,covA1A2_DZ_Pos
                    ,covA1A2_DZ_Pro, dataDZ )
 
#Ask for confidence intervals
CI <-mxCI(reference = c('aPos','cPos','ePos','aPro','cPro','ePro'), interval = 0.95, type=c("both", "lower", "upper"))
 
obj <- mxFitFunctionMultigroup(c("MZ","DZ"))
modelACE.Cholesky <- mxModel(model="ACE_Cholesky", modelMZ, modelDZ, obj, CI )
 
fitACE.Cholesky <- mxRun(modelACE.Cholesky, intervals=TRUE)
sumACE.Cholesky <- summary(fitACE.Cholesky)

Output:

> sumACE.Cholesky
Summary of ACE_Cholesky 
 
free parameters:
         name matrix          row          col      Estimate  Std.Error A
1        aPos   MZ.A    POSFAC6_1       A1_Pos  3.367661e-01 0.15893533  
2        cPos   MZ.A    POSFAC6_1       C1_Pos -5.412375e-01 0.07774325  
3        ePos   MZ.A    POSFAC6_1       E1_Pos  6.267074e-01 0.03372648  
4        aPro   MZ.A ObsPInitia_a       A1_Pro -6.837012e-09 0.25882084  
5        cPro   MZ.A ObsPInitia_a       C1_Pro -5.235314e-01 0.03486734  
6        ePro   MZ.A ObsPInitia_a       E1_Pro  6.592771e-01 0.02251672  
7 meanFAC_Pos   MZ.M            1    POSFAC6_1  3.414822e+00 0.03480667  
8 meanFAC_Pro   MZ.M            1 ObsPInitia_a  5.931390e-01 0.04724715  
 
confidence intervals:
         lbound      estimate     ubound note
aPos -0.5662877  3.367661e-01  0.5662876     
cPos -0.6643146 -5.412375e-01 -0.3565894     
ePos  0.5653805  6.267074e-01  0.6928926     
aPro -0.3805740 -6.837012e-09  0.3805739     
cPro -0.5917948 -5.235314e-01 -0.3940754     
ePro  0.6099563  6.592771e-01  0.7059469     
 
Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:              8                   1853              4606.011
   Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 490/1861
 
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       900.0111               4622.011                       NA
BIC:     -6872.2221               4655.566                 4630.175
CFI: NA 
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Need unique labels for cross-paths

The problem is here (and likewise for twin #2):

# path coefficients for twin 1. 
pathAceT1_Pos <- mxPath( from=c("A1_Pos","C1_Pos","E1_Pos"), 
                         to="POSFAC6_1", arrows=1, free=c(T,T,T), 
                         values=c(StartVar.ACE,StartVar.ACE,StartVar.ACE),
                         label=c("aPos","cPos" ,"ePos"))
 
 
 
pathAceT1_Pro <- mxPath( from=c("A1_Pos","C1_Pos","E1_Pos","A1_Pro","C1_Pro","E1_Pro"), 
                         to="ObsPInitia_a", 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("aPos","cPos" ,"ePos","aPro","cPro" ,"ePro") )

The "cross-paths" from "A1_Pos", "C1_Pos", and "E1_pos" to "ObsPInitia_a" have the same labels as the paths from those latents to "POSFAC6_1". Path coefficients with the same labels are the same free parameters. Change the first three labels in pathAceT1_Pro, and do likewise in pathAceT2_Pro further in your script.

Also, note that some of the one-headed path coefficients from latents to manifests are indeterminate with regard to sign (positive or negative). That's why some of your confidence intervals have lower and upper limits of the same absolute value but with opposite sign. It looks like you're giving all those paths positive start values, so you could put lbounds of zero on those paths. But, on the other hand, unnecessary parameter bounds can interfere with CI calculation. Really, the object of inference here ought to be the variance component, not the one-headed path coefficient. For example, the additive-genetic variance in "POSFAC6" would be the square of the coefficient labeled "aPos". You could add MxAlgebras to calculate the variance components, and request CIs for those algebras.

lior abramson's picture
Offline
Joined: 07/21/2017 - 13:13
Thank you so much!

Thank you so much!
First, regarding your second comment, I think that in order not to mess with the calculations I will leave it as it is. As you said, what is important is the square of the coefficients and I can get that regardless of the path's sign.

Second, I changed the labels as you suggested and it worked. I now have information about the contribution of one trait's ACE components to the second trait. However, the squares of the ACE components for each trait (which represent the percentage of explained variance) do not sum up to 100%. Is it reasonable?

I attach the output.

Thank you again for the prompt and detailed response

Summary of ACE_Cholesky

free parameters:
          name matrix          row          col      Estimate  Std.Error A
1    aPos_temp   MZ.A    POSFAC6_1       A1_Pos  4.379781e-01 0.10897105  
2      aPos_PB   MZ.A ObsPInitia_a       A1_Pos  1.486165e-01 0.04885548  
3    cPos_temp   MZ.A    POSFAC6_1       C1_Pos  4.932675e-01 0.07982890  
4      cPos_PB   MZ.A ObsPInitia_a       C1_Pos  2.741930e-02 0.04430686  
5    ePos_temp   MZ.A    POSFAC6_1       E1_Pos  6.162058e-01 0.03065396  
6      ePos_PB   MZ.A ObsPInitia_a       E1_Pos  3.329547e-02 0.01956340  
7         aPro   MZ.A ObsPInitia_a       A1_Pro -2.072819e-07 0.08053711  
8         cPro   MZ.A ObsPInitia_a       C1_Pro  1.297346e-01 0.03718076  
9         ePro   MZ.A ObsPInitia_a       E1_Pro  2.825964e-01 0.01157627  
10 meanFAC_Pos   MZ.M            1    POSFAC6_1  3.416948e+00 0.03527785  
11 meanFAC_Pro   MZ.M            1 ObsPInitia_a  6.066708e-01 0.01270764  

confidence intervals:

                lbound      estimate     ubound note
aPos_temp  0.132270347  4.379781e-01 0.62266662     
cPos_temp  0.286730999  4.932675e-01 0.62584616     
ePos_temp  0.559349116  6.162058e-01 0.67898025     
aPos_PB    0.026593393  1.486165e-01 0.22889358     
cPos_PB   -0.086990898  2.741930e-02 0.10476178     
ePos_PB   -0.003577813  3.329547e-02 0.07306987     
aPro      -0.140971666 -2.072819e-07 0.14097156     
cPro      -0.178035725  1.297346e-01 0.17803564     
ePro       0.259936529  2.825964e-01 0.30535018     
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
two things
Second, I changed the labels as you suggested and it worked. I now have information about the contribution of one trait's ACE components to the second trait. However, the squares of the ACE components for each trait (which represent the percentage of explained variance) do not sum up to 100%. Is it reasonable?

You seem to be confused about two things. First, since you haven't standardized the phenotypes ahead of time, it should not be the case that the ACE variance components will sum to 1.0 for each trait. They're raw variance components, not variance proportions. You may find mxStandardizeRAMpaths() useful. Second, it's not quite as simple as squaring the coefficients. Using your parameter labels, the (say) additive-genetic variance in "POSFAC6" would be aPos_temp^2, whereas the additive-genetic variance in "ObsPInitia" would be aPos_PB^2 + aPro^2. The genetic covariance between the two would be aPos_temp * aPos_PB.

Log in or register to post comments