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