# Bivariate Cholesky model with paths specification

4 posts / 0 new
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

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 
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.

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     
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.