# Bivariate Cholesky model with paths specification

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

## Need unique labels for cross-paths

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

Log in or register to post comments

In reply to Need unique labels for cross-paths by AdminRobK

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

Log in or register to post comments

In reply to Thank you so much! by lior abramson

## two things

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