bivariate heritability

biv.rmd

---
title: "Bivariate heritabilities greater than 1"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: pdf_document
---

the model -

$$
\begin{aligned}
y_a &= g_a + u_a + e_a \\
y_b &= g_b + u_b + e_b \\
g_a &\sim g_b \sim MVN(0, \Sigma^2_g) \\
u_a &\sim u_b \sim MVN(0, \Sigma^2_u) \\
e_a &\sim N(0, 1 - \sigma^2_{g_a} - \sigma^2_{u_a}) \\
e_b &\sim N(0, 1 - \sigma^2_{g_b} - \sigma^2_{u_b}) \\
\Sigma^2_g &= \begin{bmatrix}
\sigma^2_{g_a} & \rho_g \\ 
\rho_g & \sigma^2_{g_b}
\end{bmatrix} \\
\Sigma^2_u &= \begin{bmatrix}
\sigma^2_{u_a} & \rho_u \\ 
\rho_u & \sigma^2_{u_b}
\end{bmatrix}
\end{aligned}
$$

$y_a$ is phenotype A, $y_b$ is phenotype B, etc.

there are two ways in which the confounder is being specified

1. how much of the trait variance it explains in $y_a$ and $y_b$ (call it $\sigma^2_{u_a}$ and $\sigma^2_{u_b}$ - variance of a and b, respectively, explained by the confounder)
2. the correlation between $u_a$ and $u_b$ (call it $\rho_u$, you have called it `U` in the results table)

imagine if $\sigma^2_{u_a}$ and $\sigma^2_{u_b}$ are large, i.e. $u$ has a large impact on the correlation between $y_a$ and $y_b$, as a consequence the $\rho_u$ can be smaller because it has a larger weight. vice versa, if $\sigma^2_u$ is small then $\rho_u$ has to be bigger to influence the phenotypic correlation more.

you have a situation where there are actually an infinite number of parameters for $\sigma^2_{u_a}$, $\sigma^2_{u_b}$ and $\rho_u$ that can satisfy your results. the phenotypic correlation can be described as

$$
\begin{aligned}
r_p &= \sqrt{h^2_a}\sqrt{h^2_b}r_g + r_u \\
\frac{r_p - \sqrt{h^2_a}\sqrt{h^2_b}r_g}{\sqrt{u^2_a}\sqrt{u^2_b}} &= r_u
\end{aligned}
$$

you have values for $r_p$, $h^2_a$, $h^2_b$ and $r_g$, so you can find the range of values that $\sigma^2_{u_a}$, $\sigma^2_{u_b}$ and $r_u$ can take to satisfy this equation. for example:


```{r }
rp <- 0.45
ha <- 0.53
hb <- 0.6
rg <- 0.94

dat <- expand.grid(
	ua = seq(0, 1, by=0.01),
	ub = seq(0, 1, by=0.01)
)

dat$lhs <- rp - sqrt(ha) * sqrt(hb) * rg
dat$ru <- dat$lhs / dat$ua / dat$ub

# Any value of ru that is outside -1,1 is not possible so set to missing
dat$ru[dat$ru < -1] <- NA
```

Ok so this can be plotted to get an idea of the permissible parameter space -

```{r fig.cap="The surface of this plot represents possible paramters that u2a, u2b and ru can take that would satisfy the estimates of rp, h2a, h2b and rg"}

library(lattice)
wireframe(ru ~ ua + ub, 
	dat=dat, 
	drape=TRUE, 
	scales = list(arrows=FALSE), 
	ylab=expression(u[b]^2), xlab=expression(u[a]^2), zlab=expression(rho[u]))
```