bivariate heritability and coheritability

biv.R

library(MASS)


#' Simulate two traits with shared genetic and environmental effects
#'
#' @param nid sample size
#' @param  nsnp number of snps
#' @param  hsq1 h2 of trait 1
#' @param  hsq2 h2 of trait 2
#' @param  rg genetic correlation
#' @param  re non-genetic correlation
#' @param  vp1 variance of trait 1
#' @param  vp2 variance of trait 2
#' @export 
#' @return list
simulate_bivar <- function(nid, nsnp, hsq1, hsq2, rg, re, vp1, vp2)
{
	# Create SNPs
	snps <- matrix(rbinom(nid*nsnp, 2, 0.5), nid, nsnp)

	# Create effects
	eff <- mvrnorm(nsnp, c(0,0), matrix(c(1,rg,rg,1), 2))
	cor(eff)

	# create genetic values
	g1 <- scale(snps %*% eff[,1]) * sqrt(hsq1) * sqrt(vp1)
	g2 <- scale(snps %*% eff[,2]) * sqrt(hsq2) * sqrt(vp2)

	# create non-genetic values
	e <- mvrnorm(nid, c(0,0), matrix(c(1,re,re,1), 2))
	e1 <- scale(e[,1]) * sqrt(1-hsq1) * sqrt(vp1)
	e2 <- scale(e[,2]) * sqrt(1-hsq2) * sqrt(vp2)

	# create phenotypes
	y1 <- g1 + e1
	y2 <- g2 + e2

	return(list(dat = data.frame(y1=y1, y2=y2, g1=g1, g2=g2, e1=e1, e2=e2), eff = eff))
}

detach(dat$dat)
nid <- 1000
nsnp <- 100
hsq1 <- 0.5
hsq2 <- 0.7
rg <- 0.8
re <- 0.2
vp1 <- 10
vp2 <- 5
dat <- simulate_bivar(nid, nsnp, hsq1, hsq2, rg, re, vp1, vp2)
attach(dat$dat)


# Coheritability
coher <- cor(g1, g2) * sqrt(hsq1*hsq2)

# Bivariate heritability
biher <- coher / cor(y1, y2)

# Coheritability as defined by Janssens 1978 is the same as bivariate heritability
cov(g1,g2) / cov(y1,y2)
cov(g1,g2) / cov(y1,y2)
coher
biher

# How to interpret either of these values?
# Variance of trait x explainable by gx = heritability of x
cor(g1, y1)^2
cor(g2, y2)^2

# How much of trait y can be explained by genetics of trait x
cor(g1, y2)^2
coher^2 / hsq1

cor(g2, y1)^2
coher^2 / hsq2
biher^2 / hsq2 * cor(y1,y2)^2

cor(g1,g2) * sd(g1) * sd(g2) / (sd(y1)*sd(y2))
cor(g1,g2) * sd(g1) * sd(g2) / sqrt(vp1*vp2)
coher

# Bivariate heritability is proportion of phenotypic similarity that is due to genetic similarity

# Coheritability is proportion of phenotypic variance of both traits that is due to shared genetic effects

# Interpretation:
# With perfect knowledge of g1 the maximum prediction accuracy of trait 2 is coher^2 / hsq2


nid <- 4000
nsnp <- 100
hsq1 <- 0.5
hsq2 <- 0.7
rg <- 0.8
re <- 0.2
vp1 <- 1
vp2 <- 1
nsim <- 1000
run_sim <- function(nsim, nid, nsnp, hsq1, hsq2, rg, ve, vp1, vp2)
{
	l <- list()
	for(i in 1:nsim)
	{
		cat(i, "\n")
		l[[i]] <- simulate_bivar(nid, nsnp, hsq1, hsq2, rg, re, vp1, vp2)$dat
	}
	return(l)
}

sim1 <- run_sim(nsim, nid, nsnp, hsq1, hsq2, rg, ve, vp1, vp2)

var(sapply(sim1, function(x)
{
	cov(x$g1, x$g2)
}))

coher <- function(x)
{
	return(with(x, cor(g1, g2) * sqrt(var(g1)/var(y1)*var(g2)/var(y2))))
}

biher <- function(x)
{
	return(with(x, cor(g1, g2) * sqrt(var(g1)/var(y1)*var(g2)/var(y2)) / cor(y1,y2) ))
}


cov(sapply(sim1, coher), sapply(sim1, function(x) cor(x$y1, x$y2)))

var(sapply(sim1, biher))


cov(sapply(sim1, function(x)cov(x$g1,x$g2)), sapply(sim1, function(x)cov(x$e1,x$e2)))



vcg <- var(sapply(sim1, function(x) cov(x$g1, x$g2)))
rp <- with(sim1[[1]], cor(y1,y2))
cg <- with(sim1[[1]], cov(g1,g2))

vcg / rp^2 * (1 - 2*cg/rp)

Source  Variance        SE
V(G)_tr1        0.001934        0.000246
V(G)_tr2        58.393773       12.997064
C(G)_tr12       0.134715        0.041772
V(e)_tr1        0.001893        0.000203
V(e)_tr2        152.310709      12.407324
C(e)_tr12       0.127962        0.036919
Vp_tr1  0.003827        0.000131
Vp_tr2  210.704482      6.983511
V(G)/Vp_tr1     0.505424        0.056253
V(G)/Vp_tr2     0.277136        0.059087
rG      0.400858        0.105958
logL    -1568.523
n       3840