Hi all,

The short form of my question is: is there a way to model a bivariate model using umxACEv, where the model for one variable is ACE (or its submodels) and the model for the other variable is ADE (or its submodels)?

It seems from this thread that this is theoretically possible, but I'm unsure if the umx syntax enables it.

To explain a bit more, I have tried to model a bivariate ACE model, but got a negative C variance estimate (see below). Following the recommendation here I have checked what happens when fixing the relevant C component to zero (as well as rC, using umxModify), but fixing it to zero results in a significantly worse model, meaning that the C component is "significantly negative". Indeed, when estimating the univariate models for each phenotype, the model that best fits the variable with negative C is an ADE model. How can I then test a bivariate model for these two different phenotypes?

Thanks,

Noam

selDVs <- c("WM", "Mnt") model <- umxACEv(selDVs = selDVs, sep = "_", zyg = "zyg", dzData = "DZ", mzData = "MZ",nSib = 2, data = TwinData.gen.wide, addCI = TRUE, type = "FIML") umxCI(model, run = "yes") Table: Standardized parameter estimates from a 2-factor Direct variance ACE model. A: additive genetic; C: common environment; E: unique environment. | | A1|A2 | C1|C2 | E1|E2 | |:----|-----:|:----|-----:|:----|----:|:----| |WM_ | 0.78|NA | -0.36|NA | 0.58|NA | |Mnt_ | -0.20|0.38 | 0.19|0.15 | 0.09|0.47 | Table: Means (from model$top$expMean) | | WM_1| Mnt_1| WM_2| Mnt_2| |:---------|----:|-----:|----:|-----:| |intercept | 0.02| -0.02| 0.02| -0.02| lbound estimate ubound lbound Code ubound Code top.A_std[1,1] 0.275 0.778 1.181 0 0 top.A_std[2,1] -0.516 -0.201 0.105 0 0 top.A_std[1,2] -0.516 -0.201 0.105 0 0 top.A_std[2,2] -0.092 0.379 0.731 0 0 top.C_std[1,1] -0.639 -0.356 -0.046 0 0 top.C_std[2,1] -0.010 0.189 0.387 0 0 top.C_std[1,2] -0.010 0.189 0.387 0 0 top.C_std[2,2] -0.108 0.152 0.430 0 0 top.E_std[1,1] 0.418 0.578 0.808 0 0 top.E_std[2,1] -0.034 0.094 0.232 0 0 top.E_std[1,2] -0.034 0.094 0.232 0 0 top.E_std[2,2] 0.331 0.470 0.691 0 0 expMean_WM_1 -0.195 0.018 0.231 0 0 expMean_Mnt_1 -0.194 -0.021 0.152 0 0 A_r1c1 2.313 6.661 10.482 0 0 A_r2c1 -3.305 -1.264 0.665 0 0 A_r2c2 -0.424 1.756 3.445 0 0 C_r1c1 -5.644 -3.050 -0.391 0 0 C_r2c1 -0.065 1.189 2.487 0 0 C_r2c2 -0.502 0.703 2.028 0 0 E_r1c1 3.604 4.947 6.940 0 0 E_r2c1 -0.221 0.591 1.488 0 0 E_r2c2 1.559 2.178 3.186 0 0 umxSummaryACEv(model, digits = 2, comparison = NULL, std = TRUE, showRg = TRUE, CIs = TRUE, report = c("markdown"), file = getOption("umx_auto_plot"), returnStd = FALSE, extended = FALSE, zero.print = ".", show = c("std", "raw") ) Table: Standardized parameter estimates from a 2-factor Direct variance ACE model. A: additive genetic; C: common environment; E: unique environment. | | A1|A2 | C1|C2 | E1|E2 | |:----|-----:|:----|-----:|:----|----:|:----| |WM_ | 0.78|NA | -0.36|NA | 0.58|NA | |Mnt_ | -0.20|0.38 | 0.19|0.15 | 0.09|0.47 | Table: Means (from model$top$expMean) | | WM_1| Mnt_1| WM_2| Mnt_2| |:---------|----:|-----:|----:|-----:| |intercept | 0.02| -0.02| 0.02| -0.02| Table: Genetic correlations | | rA1|rA2 |rC1 |rC2 | rE1|rE2 | |:----|-----:|:---|:---|:---|----:|:---| |WM_ | 1.00| | | | 1.00| | |Mnt_ | -0.37|1 | | | 0.18|1 |

Unfortunately, umxACEv has been programmed in a way that makes it difficult to convert for mixed variance component models (ACE for variable 1 and ADE for variable 2). The issue is in the calculation of the expected covariances for DZ twins, as can be seen when inspecting the bivariate model from the umxACEv examples.

Essentially, the Kronecker product forces all the C elements to be multiplied by 1. Changing dzCr to .25 would make both variables ADE modeled. Probably the dzCr matrix needs changing to a matrix that has elements c(1,0,0,.25) and the formula for expCovDZ changed to

`(dzAr %x% A) + (dzCr * C)`

(the off-diagonal zeroes are because C and D are assumed to be uncorrelated). It would also be necessary to change similarly the MZ group mzCr matrix and formula (to an identity matrix and again use * instead of %x%)Hi Mike!

Is this worth doing? Seems to me that as D and C can't be estimated simultaneously in these twin-only data, the consequences of not allowing dominance effects but not allowing them to enter into trait covariances (unlikely traits have covariance due to A but not to D?) and forcing the other traits to have zero C is a worse bargain than an ACE model where the trait supposedly with D and no C is forced to have its D estimated to some degree in the A partition? Has anyone argued the toss in print? Made me think of dear Lindon (RIP).

Hi Tim

The more variables in a multivariate ACE model, the more likely it is that C will be negative for some cases. If we really think that is due to genetic dominance or AxA non-additivity, then I don't think it's cool to have what is effectively C for some variables (rDZ > .5rMZ) and D or AxA interaction for others. In particular, rC between traits that have such different patterns should be fixed to zero on the assumption that there is no C for one trait or no D for the other. One might argue that in the event that such C-D covariances are found empirically (so a likelihood ratio test would significantly object to fixing the rC or rD to zero), they imply that there's both C and D operating for at least one of the traits. Possibly, an option to automatically fix rCs between variables that have positive-negative mismatch on the C component would be helpful here.

Thank you both; I appreciate the discussion on this topic.

I understand now that combining two phenotypes that each "responds" to another model (ACE / ADE ) cannot be fairly recommended (especially in my case, where my theoretical point of interest is the relations themselves; rA, rC, etc.)

Are there reading materials or other references to this issue? That would be a great help to me.