interpreting a bivariate ACE model with negative variance using umxACEv
noamMai
Joined: 03/20/2023
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](https://openmx.ssri.psu.edu/thread/4172) 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](https://openmx.ssri.psu.edu/node/4748) 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
> m1$top$hAC
mxAlgebra 'hAC'
$formula: (dzAr %x% A) + (dzCr %x% C)
$result:
[,1] [,2]
[1,] 26.207 2.9785
[2,] 2.978 0.5264
dimnames: NULL
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%)Log in or register to post comments
a good bargain?
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).
Log in or register to post comments
In reply to a good bargain? by tbates
May be worthwhile, especially multivariate
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.
Log in or register to post comments
Thanks!
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.
Log in or register to post comments