---
title: 'GWAS 4: Bayesian inference'
author: "Matti Pirinen, University of Helsinki"
date: "Updated: 2-Mar-2023."
output:
html_document: default
urlcolor: blue
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
set.seed(19)
```
This document is licensed under a
[Creative Commons Attribution-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by-sa/4.0/).
The slide set referred to in this document is "GWAS 4".
A crude statistical inference in the GWAS framework
is based on a fixed significance threshold and only determines whether
the P-value falls below the GWS threshold (typically 5e-8).
However, we almost always make more efficient use of the data by reporting
effect sizes, their standard errors and numerical values for P-values.
Suppose that there are two variants that both pass the GWS threshold, first
yielding a P-value of 4e-8 and the second of 2e-200.
How could we quantify the probability that each of these variants is having a non-zero effect?
We have previously learned that
the P-value itself is not a probability on the null hypothesis.
Let's next look what information and assumptions, in addition to the P-value,
we would need in order to talk about a probability of a non-zero effect.
Here is a review on Bayesian methods in GWAS by [Stephens & Balding](https://www.nature.com/articles/nrg2615).
### From significance to the probability of association
In order to talk about probability of a hypothesis, we need to define the
set of all possible hypotheses whose combined probability is 1, i.e.,
we work under the assumption that one of the hypotheses is true.
Note that for P-value calculation, we only ever define
the null hypothesis and therefore we are unable to talk about the probability
of the null hypothesis because we haven't defined what else is possible:
Even if the data seemed unlikely under the null hypothesis,
if the data were even more unlikely under other possible hypotheses, then
the null hypothesis might still be quite probable, and such
considerations cannot be done from the P-value alone.
In our case, we consider only two hypotheses: $H_0: \beta =0$ and $H_1: \beta \neq 0.$
Next we need to quantify the **prior probabilities** of these hypotheses.
These answer to the question: What would I consider as probability of each
hypothesis before I have seen the data.
The phrase "What would _I_" is there on purpose: prior probabilities are subjective. They are
based on whatever knowledge _I_ have available. Therefore different persons may have different
prior probabilities for the same hypothesis and my prior can (and will!)
change as I learn more about the question.
For example,
$P(H_1) = 10^{-5}$ could be a reasonable prior for a non-zero effect based on
what I know about GWAS.
(Last week we saw how a magnitude more stringent assumption $P(H_1) = 10^{-6}$
led us to the common GWS threshold of 5e-8.)
Then we observe the data $\mathcal{D}$ and our interest is in the probabilities
of each hypothesis after we have seen the data. This is the core question of **Bayesian
inference**: How does observing the data update our beliefs from our current state of knowledge,
described by prior probabilities $P(H_i)$,
into our **posterior probabilities** $P(H_i|\mathcal{D})$?
In short: How do we learn from data?
Not surprisingly, the answer is the Bayes rule.
#### 4.1 Bayes rule
To write down the Bayes rule (also called Bayes theorem, Bayes formula),
we just remember that we are considering joint distributions of two variables.
Here the two variables are the hypothesis $H_i$ and the observed data $\mathcal{D}$.
We expand their joint probability $P(H_i,\mathcal{D})$ using conditional probability rule
in both ways possible
$$P(H_i \,|\, \mathcal{D}) P(\mathcal{D})= P(H_i,\mathcal{D}) = P(\mathcal{D}\, |\, H_i) P(H_i),$$
from which we can solve the conditional probability that the hypothesis
holds given that we have observed the data $\mathcal{D}$:
$$
P(H_i|\mathcal{D}) = \frac{P(\mathcal{D}\,|\,H_i)P(H_i)}{P(\mathcal{D})},\qquad\textrm{ for } i=0,1.
$$
This is the Bayes rule. The conditional probability on the
left hand side is also called the **posterior probability**
of the hypothesis given the data.
Since the term $P(\mathcal{D})$ (the marginal probability of the observed data) does
not depend on hypothesis $H_i$, we can get rid of it by taking the ratio of the posteriors of
the two hypotheses:
$$
\frac{P(H_1\,|\,\mathcal{D})}{P(H_0\,|\,\mathcal{D})} = \frac{P(\mathcal{D}\,|\,H_1)P(H_1)}{P(\mathcal{D}\,|\,H_0)P(H_0)}.
$$
Hence, in order to compute the posterior probability ratio for the hypotheses,
we will still need the terms $P(\mathcal{D}\,|\,H_i)$ in addition to the prior probabilities.
$P(\mathcal{D}\,|\,H_i)$ describes what kind of data sets
we are likely to see under each hypothesis and with which probability.
After these probability densities are specified, the inference
is about letting the possible hypotheses compete both in how well they explain the
observed data (terms $P(\mathcal{D}\,|\,H_1)$ and $P(\mathcal{D}\,|\,H_0)$)
and in how probable they are *a priori* (prior probablity terms $P(H_1)$ and $P(H_0)$).
**Example 4.1.** The Bayesian inference shows that both the observed data AND the
prior knowledge is crucial for a complete inference.
Suppose, for example, that a sequencing effort of my genome returns data
$\mathcal{D}$ that seem to miss chromosome 6 completely.
We have two hypotheses: $H_0$: "There is a technical error",
or $H_1$: "I don't carry any copies of chromosome 6 (in the cells involved)".
The observed result could have
resulted from either of these options
and hence under both hypotheses $P(\mathcal{D}\,|\,H_i)$ is relatively high,
indicating that the data are consistent with both hypotheses.
However, the prior odds of $P(H_1)/P(H_0)$ are likely small
if we think that it is much more likely that there is a technical error somewhere
than that I would be missing chr 6 and still be alive and even pretty healthy.
Hence, the posterior conclusion that combines the prior probabilities
and the observed data is that it is more probable that
we have an error somewhere in the process than that I don't carry chr 6,
even though the observation alone couldn't tell the two hypotheses apart.
**Example 4.2.** Interpretation of a medical test result is [the standard
example](https://en.wikipedia.org/wiki/Bayes%27_theorem#Drug_testing) of the use of Bayes rule.
Let's apply it to a case where we try to determine whether an individual has a disease given his genotype.
Suppose that each copy of *HLA-DRB1\*1501* allele on chromosome 6 increases the risk of *multiple sclerosis*
by OR=3. Prevalence of MS-diseases is $K=0.001$ and population frequency of DRB1\*1501 is 0.20.
What is probability of ever getting the disease for each genotype (i.e. 0,1 or 2 copies of DRB\*1501)?
**Answer.** Denote by $D$ the event of getting the disease and by $X$ the genotype.
Here $D$ has the role of a hypothesis and $X$ the role of observed data in the above formulation of Bayes rule.
Bayes rule says that for each genotype $x \in \{0,1,2\}$:
$$P(D\,|\, X = x) = \frac{P(D)P(X=x\, |\, D)}{P(X=x)}.$$
We know that $P(D)=K=0.001$ and we can assume that the control frequencies
are approximately the population frequencies since the diseases has so low prevalence.
Let's use function `case.control.freqs()` from GWAS 1 to
determine the case frequencies assuming HWE in control population:
```{r, echo = FALSE}
case.control.freqs <- function(q, or){
#if dimension of 'q' is 1 then 'q' is taken as allele 1 freq. in controls
#if dimension of 'q' is 3 then 'q' is taken as the genotype (0,1,2) frequencies in controls
#if dimension of 'or' is 1, then 'or' is per each copy of allele 1
#if dimension of 'or' is 2, then 'or[1]' is for genotype 1 vs.0 and 'or[2]' is geno 2 vs. 0
if(length(q) == 1) q = c((1-q)^2, 2*q*(1-q), q^2)
stopifnot(length(q) == 3)
if(length(or) == 1) or = c(or, or^2)
stopifnot(length(or) == 2)
f = q / q[1] * c(1, or)
f = f / sum(f)
data.frame(cases = f, controls = q, row.names = c(0,1,2))
}
```
```{r}
K = 0.001
or = 3
f = 0.2
cc.f = case.control.freqs(f, or)
```
The risk to get the disease given the genotype is, according to Bayes rule
```{r}
rbind(genotype = c(0, 1, 2), risk = K * cc.f$cases / cc.f$controls)
```
So even though the disease status has a large effect on the genotype distributions,
still even the high risk group has risk < 0.5% and the genotype doesn't give
a practically useful predictive accuracy at the level of an individual from the population because
the disease is so rare.
#### 4.2 Probability model for observed GWAS data
To use Bayes rule in GWAS setting, we need to define a probability density function for the
observed data under both the null hypothesis and the alternative hypothesis.
In a linear regression GWAS model
$y = \mu + x \beta + \varepsilon$,
the observed data $\mathcal{D}_k=(\pmb{y},\pmb{x}_k)$ consist of the phenotype vector
$\pmb{y}$ and the vector of genotypes $\pmb{x}_k$ at the tested variant $k$.
When we assume Gaussian errors
(i.e. each $\varepsilon_i$ having a Normal distribution with same variance $\sigma^2$)
the probability density of data for a fixed value of effect size $\beta$ and error variance
$\sigma^2$ is $$p\left(\mathcal{D}_k\, |\, \beta, \sigma^2 \right) =
\mathcal{N}(\pmb{y};\,\pmb{x}_k\beta,\,\sigma^2 I) \propto
\exp\left(-(\pmb{y}-\pmb{x}_k\beta)^T(\pmb{y}-\pmb{x}_k\beta)/(2\sigma^2)\right).
$$
We have ignored $\mu$ here by assuming that $y$ and $x_k$ are mean-centered values; this
just simplifies the notation but doesn't affect the inference on parameter $\beta$.
Under the null model, we set $\beta=0$ and
in the alternative model,
we can set $\beta$ to some other value $b_1$. If we do not want to specify
our model of true effect by a single value $b_1$, we can use a **prior
distribution** for $\beta$, for example, by saying that under the alternative model
$\beta \sim \mathcal{N}(b_1,\tau_1^2)$.
This means that if the alternative model holds, then the true effect sizes
are distributed around value $b_1$ with a standard deviation of $\tau_1$.
With this prior distribution, the probability of data under $H_1$
is given by weighting the above likelihood function by the
prior probability of each possible
value of $\beta$: $$p(\mathcal{D}_k\,|\,H_1) = \int_\beta
p\left(\mathcal{D}_k\,|\,\beta,\sigma^2\right) p(\beta\,|\,H_1) d\beta = \int_\beta
\mathcal{N}\left(\pmb{y};\,\pmb{x}_k\beta,\,\sigma^2\right)
\mathcal{N}\left(\beta;\, b_1,\, \tau_1^2\right) d\beta.$$
(For simplicity, in both models, we may fix $\sigma^2$ to its empirical maximum likelihood estimate.)
If we assume that the mean parameter $b_1=0$ in the Gaussian prior of $\beta$,
then the integral can be done analytically to give
$$p(\mathcal{D}_k\,|\,H_1) =
c\cdot \mathcal{N}\left(\widehat{\beta};\, 0, \, \tau_1^2 + \textrm{SE}^2
\right),$$
where $c$ is a constant and $\widehat{\beta}$ is the MLE of $\beta$
and SE the corresponding standard error. Note that by replacing $\tau_1$ with 0,
we have
$$p(\mathcal{D}_k\,|\,H_0) = c\cdot \mathcal{N}\left(\widehat{\beta};\, 0, \,
\textrm{SE}^2 \right).$$
These results show that quantity $p(\mathcal{D}_k\,|\,H_i)$ that tells how well
model $H_i$ explains the data can be evaluated (up to a constant)
by asking how well the model can explain the observed estimate
$\widehat{\beta}$.
The null model assumes that the estimate is a result of
sampling variation alone (only the SE contributes to the deviation from 0)
whereas the alternative model allows some additional variance $\tau_1^2$ around 0
due to the truely non-zero effect.
Let's demonstrate this.
**Example 4.3.**
```{r}
n = 3000 # sample size for SE calculation
f = 0.3 # MAF for SE calculation
sigma = 1 # error SD
se = sigma / sqrt(2*f*(1-f)*n) # SE for QT GWAS
tau = 0.1 # prior standard deviation for effect size beta under H1
# Let's draw probability densities of "data" under the two models, H0 and H1
# as a function of MLE estimate
x = seq(-0.3, 0.3, by = 0.01)
y1 = dnorm(x, 0, sqrt(tau^2 + se^2) )
y0 = dnorm(x, 0, se)
plot(x, y0, t = "l", col = "cyan", lwd = 2, xlab = expression(hat(beta)),
ylab = "probability density of data")
lines(x, y1, col = "magenta", lwd = 2)
legend("topright", c("H0", "H1"), col = c("cyan","magenta"), lwd = 2)
# We make a shortcut and don't simulate data at all, but we directly simulate estimates
# Suppose we have two cases, first follows null (beta = 0), second follows alternative (beta = 0.2)
b =c(0, 0.2)
b.est = rnorm(2, b, se) #these are simulated estimates: true means and Gaussian noise determined by SE
points(b.est, c(0, 0), pch = 19, col = c("blue", "red") )
#Next: log of Bayes factor of H1 vs H0, will be explained below
# use log-scale to avoid inaccuracies.
log.bf.10 = dnorm(b.est, 0, sqrt(tau^2 + se^2), log = T ) - dnorm(b.est, 0, se, log = T)
bf.10 = exp(log.bf.10) #then turn from log-scale to Bayes factor scale
text(b.est[1], 4, signif(bf.10[1], 2), col = "blue")
text(b.est[2], 4, signif(bf.10[2], 2), col = "red")
```
The distribution under $H_0$ describes what kinds of effect estimates we can get with this
sample size and MAF when the true effect is exactly 0. Any deviation from 0 is then by
statistical sampling effect, as quantified by the SE.
The distribution under $H_1$ describes what kinds of effect estimates we expect when we have
BOTH a non-zero effect size, whose expected range is described by the standard deviation $\tau_1$,
and we have ALSO a statistical sampling effect as quantified by the SE.
If the observed estimate $\widehat{\beta}$ is close to 0, then $H_0$ explains the data better than
$H_1$, whereas the opposite is true when $\widehat{\beta}$ is farther away from
0. With these parameters, $H_1$ starts to dominate roughly when $|\widehat{\beta}|
\geq 0.05$. Values close to zero are relatively more likely under the null
than under the alternative.
This is because the the null model can only explain well data sets with effect estimates near zero
whereas the alternative model can also explain observations
that are farther away from the zero, and hence the probability mass of the alternative
is less concentrated around 0 than that of the null. Consequently,
the null model gets higher probability densities near 0 than the alternative model
that has spread its probability density more widely away from 0.
#### 4.3 Bayes factor
Two points shown in the plot above are examples of possible estimates that could result
either under $H_0$ (blue) or $H_1$ (red). The values are ratios
$P(\mathcal{D}|H_1)/P(\mathcal{D}|H_0)$ of probability densities computed at these two points.
When this ratio $<1$, the null model $H_0$ explains the data better than the alternative model
and when it is $>1$, the opposite is true. This ratio is called **Bayes factor (BF)** that is the
factor that multiplies the prior odds to get to the posterior odds: $$
\underbrace{\frac{P(H_1\,|\,\mathcal{D})}{P(H_0\,|\,\mathcal{D})}}_{\textrm{posterior odds}} =
\underbrace{\frac{P(\mathcal{D}\,|\,H_1)}{P(\mathcal{D}\,|\,H_0)}}_{\textrm{Bayes factor}} \times
\underbrace{\frac{P(H_1)}{P(H_0)}}_{\textrm{prior odds}}. $$
To interpret the Bayes factor we can think that in order that the
posterior probability of the alternative model would be at least 10-fold
compared to that of the null, the BF needs to be at least 10 times
the inverse of the prior odds.
If the prior odds are 1e-5, then we would need a BF of at least 1e+6 to result in posterior odds > 10.
We are almost there having calculated a proper probability for the null hypothesis.
We still need to agree on the prior probability of the null model.
Last week we saw that current genome-wide significance level can be thought to
correspond to a prior probability of about $\text{Pr}(H_1) = 10^{-6}$.
We know that this seems very stringent compared to
the true genetic architecture behind complex traits and diseases,
but it gives us a conservative reference point.
With this prior probability,
the posterior odds and posterior probabilities for the alternative model are:
```{r}
post.odds = bf.10 * 1e-6/(1-1e-6) # P(H_1|D) / P(H_0|D)
post.prob = post.odds / (1 + post.odds) # P(H_1|D)
paste(post.prob)
```
For an illustration, let's check the P-values corresponding to these two data sets
using a Wald statistic $\widehat{\beta}/\textrm{SE}$:
```{r}
pchisq( (b.est/se)^2, df = 1, lower = F )
```
So P-value of the first one is quite close to 0.2, whereas the second is
way beyond 5e-8. And the Bayesian analysis said that
the first one is almost certain to be null whereas the second one is
almost certain to have a non-zero effect.
Thus, there is no difference between the conclusion of the standard *P*-value-based GWAS inference
and the Bayesian inference, and this is typically the case when there is enough data available.
Conceptually, however, there is a large difference between the P-value,
which is a probability of seeing at least as extreme data under the null,
and the posterior probability of the hypothesis itself.
There were several assumptions made in the Bayesian analysis about the
effect sizes under $H_1$ and also on the prior probabilities of the models, and
the posterior probabilities will change when these assumptions are changed.
Therefore, *P*-values remain useful simple summaries of data that can
be computed easily. The important thing
is to know what *P*-values are and what they are not, and that what
kind of additional pieces of information would be needed in order to
properly quantify the probabilities of the hypotheses.
#### 4.4 Approximate Bayes factor in GWAS
The calculation of the Bayes factor above, that was
based on the maximum likelihood estimate $\widehat{\beta}$
and its SE, was proposed by
[Jon Wakefield in 2009](https://onlinelibrary.wiley.com/doi/abs/10.1002/gepi.20359). The formula is
$$
\textrm{ABF} \approx \frac{P(\mathcal{D}\,|\,H_1)}{P(\mathcal{D}\,|\,H_0)}
\approx \frac{\mathcal{N}\left(\widehat{\beta};\, 0, \, \tau_1^2 + \textrm{SE}^2
\right)}{\mathcal{N}\left(\widehat{\beta};\, 0, \, \textrm{SE}^2
\right)} = \frac{(2\pi)^{-0.5}(\tau_1^2 + \textrm{SE}^2)^{-0.5}\exp\left(-\frac{1}{2}\frac{\widehat{\beta}^2}{\tau_1^2 + \textrm{SE}^2}\right) }
{(2\pi)^{-0.5}(\textrm{SE}^2)^{-0.5}\exp\left(-\frac{1}{2}\frac{\widehat{\beta}^2}{ \textrm{SE}^2}\right) }
$$
$$
=\sqrt{\frac{\textrm{SE}^2}{\tau_1^2 + \textrm{SE}^2}}
\exp\left(\frac{1}{2} \frac{\widehat{\beta}^2}{\textrm{SE}^2} \frac{\tau_1^2}{\tau_1^2 + \textrm{SE}^2} \right)
,
$$
where the alternative model is specified by an effect size prior
$H_1: \beta \sim \mathcal{N}(0,\tau_1^2).$ We have presented the ABF in the form
where the alternative is in the numerator and null in the denominator.
Hence large ABF means strong evidence in favor of the alternative model.
(Wakefield's paper used the inverse of this quantity as Bayes factor, i.e.,
it computes Bayes factor comparing null to the alternative whereas
we compare alternative to null.)
In R, this is easy to compute using `dnorm()` function, as we did above.
It is always good to do the ratio of densities on log-scale to avoid possible numerical
underflows/overflows. In `dnorm()` this happens by adding `log = TRUE` to the command.
Then the ratio of densities becomes a difference between log-densities:
```
log.bf = dnorm(b.est, 0, sqrt(tau^2 + se^2), log = T ) - dnorm(b.est, 0, se, log = T)
bf = exp(log.bf) #turn from log-scale to Bayes factor scale
```
The same formula can be used when $\widehat{\beta}$ and its SE originate
from a disease study analyzed by logistic regression. In that case, the formula is
an approximation based on the assumption that the logistic regression likelihood has
a shape of a Gaussian density function. Therefore, this approach is generally called
**Approximate Bayes Factor (ABF)**.
ABF can be computed from the observed GWAS data ($\widehat{\beta},\textrm{SE}$)
once we have chosen the variance $\tau_1^2$
of the effect size distribution under the alternative. How should we do that?
**Example 4.4.** Let's assume that the non-zero effects have a distribution
$\mathcal{N}(0,\tau_1^2)$ and we want to determine $\tau_1$ in such a way that
with 95% probability the effect (of a SNP with MAF = 0.25)
explains less than proportion $p$ of
the phenotypic variance $v = \textrm{Var}(y)$.
We will first compute the effect size $\beta_p$ that explains exactly phenotypic
variance of $p\cdot v$, and then we will find the sd parameter $\tau_1$
for which 95% of the probability mass of $\mathcal{N}(0,\tau_1^2)$ is within the
region $(-\beta_p, \beta_p)$.
```{r}
v = 1 # Set this to the phenotypic variance
p = 0.01 # Effect explains less than 1% of the trait variance,
target.prob = 0.95 # with this probability
maf = 0.25 # minor allele frequency
# 2*maf*(1 - maf)*b^2 = p * v --> b = +/- sqrt(p*v/(2*maf*(1-maf)))
b = sqrt(p * v / ( 2 * maf * (1 - maf) ) )
tau.seq = seq(0, 1, 0.001) #grid to evaluate tau
x = pnorm(b, 0, tau.seq, lower = F) #what is the upper tail prob. at b for each value of tau?
tau.1 = tau.seq[which.min( abs(x - (1 - target.prob)/2) )] #which is closest to target prob?
#Check that the probability mass in (-b,b) is indeed close to target
print(paste0("tau.1=",tau.1," has mass ",signif(1 - 2*pnorm(b, 0, tau.1, lower = F),3),
" in (-",signif(b,4),", ",signif(b,4),")."))
```
**Example 4.5.** For case-control GWAS, we want to find such $\tau_1$ that with 95%
probability a variant can increase risk by at most OR of 1.30.
Now we get the critical point $\beta_p$ directly as $\log(1.30)$ and then proceed
as above.
```{r}
or = 1.30 # Effect is at most this large
target.prob = 0.95 # with this probability
b = log(or)
tau.seq = seq(0, 1, 0.001) # grid to evaluate tau
x = pnorm(b, 0, tau.seq, lower = F) # what is the upper tail prob. at b for each value of tau?
tau.1 = tau.seq[which.min( abs(x - (1 - target.prob)/2) )] # which is closest to target prob?
# Check that the probability mass in (-b,b) is indeed close to target
print(paste0("tau.1=",tau.1," has mass ",signif(1 - 2*pnorm(b, 0, tau.1, lower = F),3),
" in (-",signif(b,4),", ",signif(b,4),")."))
```
Note that the above choices of $\tau_1$ do not model very large effect sizes.
There are variants that explain, say, over 10% of the variation of a quantitative trait
or that have an OR of 3 for some disease. To properly model them in the Bayesian framework,
one would need to use several prior distributions and average the results
(Bayesian model averaging).