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
Printer-friendly version
The problem is here (and likewise for twin #2):
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 inpathAceT2_Profurther 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.
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
confidence intervals:
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 beaPos_temp^2, whereas the additive-genetic variance in "ObsPInitia" would beaPos_PB^2 + aPro^2. The genetic covariance between the two would beaPos_temp * aPos_PB.