User:Timothee Flutre/Notebook/Postdoc/2011/11/10: Difference between revisions

Revision as of 15:41, 3 February 2013

Project name

<html><img src="/images/9/94/Report.png" border="0" /></html> Main project page
<html><img src="/images/c/c3/Resultset_previous.png" border="0" /></html>Previous entry<html>      </html>Next entry<html><img src="/images/5/5c/Resultset_next.png" border="0" /></html>

Bayesian model of univariate linear regression for QTL detection

This page aims at helping people like me, interested in quantitative genetics, to get a better understanding of some Bayesian models, most importantly the impact of the modeling assumptions as well as the underlying maths. It starts with a simple model, and gradually increases the scope to relax assumptions. See references to scientific articles at the end.

Data: let's assume that we obtained data from N individuals. We note [math]\displaystyle{ y_1,\ldots,y_N }[/math] the (quantitative) phenotypes (e.g. expression levels at a given gene), and [math]\displaystyle{ g_1,\ldots,g_N }[/math] the genotypes at a given SNP (encoded as allele dose: 0, 1 or 2).

Goal: we want to assess the evidence in the data for an effect of the genotype on the phenotype.

Assumptions: the relationship between genotype and phenotype is linear; the individuals are not genetically related; there is no hidden confounding factors in the phenotypes.

Likelihood: we start by writing the usual linear regression for one individual

[math]\displaystyle{ \forall i \in \{1,\ldots,N\}, \; y_i = \mu + \beta_1 g_i + \beta_2 \mathbf{1}_{g_i=1} + \epsilon_i \; \text{ with } \; \epsilon_i \; \overset{i.i.d}{\sim} \; \mathcal{N}(0,\tau^{-1}) }[/math]

where [math]\displaystyle{ \beta_1 }[/math] is in fact the additive effect of the SNP, noted [math]\displaystyle{ a }[/math] from now on, and [math]\displaystyle{ \beta_2 }[/math] is the dominance effect of the SNP, [math]\displaystyle{ d = a k }[/math].

Let's now write the model in matrix notation:

[math]\displaystyle{ Y = X B + E \text{ where } B = [ \mu \; a \; d ]^T }[/math]

This gives the following multivariate Normal distribution for the phenotypes:

[math]\displaystyle{ Y | X, \tau, B \sim \mathcal{N}(XB, \tau^{-1} I_N) }[/math]

Even though we can write the likelihood as a multivariate Normal, I still keep the term "univariate" in the title because the regression has a single response, [math]\displaystyle{ Y }[/math]. It is usual to keep the term "multivariate" for the case where there is a matrix of responses (i.e. multiple phenotypes).

The likelihood of the parameters given the data is therefore:

[math]\displaystyle{ \mathcal{L}(\tau, B) = \mathsf{P}(Y | X, \tau, B) }[/math]

[math]\displaystyle{ \mathcal{L}(\tau, B) = \left(\frac{\tau}{2 \pi}\right)^{\frac{N}{2}} exp \left( -\frac{\tau}{2} (Y - XB)^T (Y - XB) \right) }[/math]

Priors: we use the usual conjugate prior

[math]\displaystyle{ \mathsf{P}(\tau, B) = \mathsf{P}(\tau) \mathsf{P}(B | \tau) }[/math]

A Gamma distribution for [math]\displaystyle{ \tau }[/math]:

[math]\displaystyle{ \tau \sim \Gamma(\kappa/2, \, \lambda/2) }[/math]

which means:

[math]\displaystyle{ \mathsf{P}(\tau) = \frac{\frac{\lambda}{2}^{\kappa/2}}{\Gamma(\frac{\kappa}{2})} \tau^{\frac{\kappa}{2}-1} e^{-\frac{\lambda}{2} \tau} }[/math]

And a multivariate Normal distribution for [math]\displaystyle{ B }[/math]:

[math]\displaystyle{ B | \tau \sim \mathcal{N}(\vec{0}, \, \tau^{-1} \Sigma_B) \text{ with } \Sigma_B = diag(\sigma_{\mu}^2, \sigma_a^2, \sigma_d^2) }[/math]

which means:

[math]\displaystyle{ \mathsf{P}(B | \tau) = \left(\frac{\tau}{2 \pi}\right)^{\frac{3}{2}} |\Sigma_B|^{-\frac{1}{2}} exp \left(-\frac{\tau}{2} B^T \Sigma_B^{-1} B \right) }[/math]

Joint posterior (1):

[math]\displaystyle{ \mathsf{P}(\tau, B | Y, X) = \mathsf{P}(\tau | Y, X) \mathsf{P}(B | Y, X, \tau) }[/math]

Conditional posterior of B:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) = \frac{\mathsf{P}(B, Y | X, \tau)}{\mathsf{P}(Y | X, \tau)} }[/math]

Let's neglect the normalization constant for now:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B) }[/math]

Similarly, let's keep only the terms in [math]\displaystyle{ B }[/math] for the moment:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto exp(B^T \Sigma_B^{-1} B) exp((Y-XB)^T(Y-XB)) }[/math]

We expand:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto exp(B^T \Sigma_B^{-1} B - Y^TXB -B^TX^TY + B^TX^TXB) }[/math]

We factorize some terms:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto exp(B^T (\Sigma_B^{-1} + X^TX) B - Y^TXB -B^TX^TY) }[/math]

Importantly, let's define:

[math]\displaystyle{ \Omega = (\Sigma_B^{-1} + X^TX)^{-1} }[/math]

We can see that [math]\displaystyle{ \Omega^T=\Omega }[/math], which means that [math]\displaystyle{ \Omega }[/math] is a symmetric matrix. This is particularly useful here because we can use the following equality: [math]\displaystyle{ \Omega^{-1}\Omega^T=I }[/math].

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto exp(B^T \Omega^{-1} B - (X^TY)^T\Omega^{-1}\Omega^TB -B^T\Omega^{-1}\Omega^TX^TY) }[/math]

This now becomes easy to factorizes totally:

[math]\displaystyle{ \mathsf{P}(B | Y, X, \tau) \propto exp((B^T - \Omega X^TY)^T\Omega^{-1}(B - \Omega X^TY)) }[/math]

We recognize the kernel of a Normal distribution, allowing us to write the conditional posterior as:

[math]\displaystyle{ B | Y, X, \tau \sim \mathcal{N}(\Omega X^TY, \tau^{-1} \Omega) }[/math]

Posterior of [math]\displaystyle{ \tau }[/math]:

Similarly to the equations above:

[math]\displaystyle{ \mathsf{P}(\tau | Y, X) \propto \mathsf{P}(\tau) \mathsf{P}(Y | X, \tau) }[/math]

But now, to handle the second term, we need to integrate over [math]\displaystyle{ B }[/math], thus effectively taking into account the uncertainty in [math]\displaystyle{ B }[/math]:

[math]\displaystyle{ \mathsf{P}(\tau | Y, X) \propto \mathsf{P}(\tau) \int \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B) \mathsf{d}B }[/math]

Again, we use the priors and likelihoods specified above (but everything inside the integral is kept inside it, even if it doesn't depend on [math]\displaystyle{ B }[/math]!):

[math]\displaystyle{ \mathsf{P}(\tau | Y, X) \propto \tau^{\frac{\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{3/2} \tau^{N/2} exp(-\frac{\tau}{2} B^T \Sigma_B^{-1} B) exp(-\frac{\tau}{2} (Y - XB)^T (Y - XB)) \mathsf{d}B }[/math]

As we used a conjugate prior for [math]\displaystyle{ \tau }[/math], we know that we expect a Gamma distribution for the posterior. Therefore, we can take [math]\displaystyle{ \tau^{N/2} }[/math] out of the integral and start guessing what looks like a Gamma distribution. We also factorize inside the exponential:

[math]\displaystyle{ \mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{3/2} exp \left[-\frac{\tau}{2} \left( (B - \Omega X^T Y)^T \Omega^{-1} (B - \Omega X^T Y) - Y^T X \Omega X^T Y + Y^T Y \right) \right] \mathsf{d}B }[/math]

We recognize the conditional posterior of [math]\displaystyle{ B }[/math]. This allows us to use the fact that the pdf of the Normal distribution integrates to one:

[math]\displaystyle{ \mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} exp\left[-\frac{\tau}{2} (Y^T Y - Y^T X \Omega X^T Y) \right] }[/math]

We finally recognize a Gamma distribution, allowing us to write the posterior as:

[math]\displaystyle{ \tau | Y, X \sim \Gamma \left( \frac{N+\kappa}{2}, \; \frac{1}{2} (Y^T Y - Y^T X \Omega X^T Y + \lambda) \right) }[/math]

Joint posterior (2): sometimes it is said that the joint posterior follows a Normal Inverse Gamma distribution:

[math]\displaystyle{ B, \tau | Y, X \sim \mathcal{N}IG(\Omega X^TY, \; \tau^{-1}\Omega, \; \frac{N+\kappa}{2}, \; \frac{\lambda^\ast}{2}) }[/math]

where [math]\displaystyle{ \lambda^\ast = Y^T Y - Y^T X \Omega X^T Y + \lambda }[/math]

Marginal posterior of B: we can now integrate out [math]\displaystyle{ \tau }[/math]:

[math]\displaystyle{ \mathsf{P}(B | Y, X) = \int \mathsf{P}(\tau) \mathsf{P}(B | Y, X, \tau) \mathsf{d}\tau }[/math]

[math]\displaystyle{ \mathsf{P}(B | Y, X) = \frac{\frac{\lambda^\ast}{2}^{\frac{N+\kappa}{2}}}{(2\pi)^\frac{3}{2} |\Omega|^{\frac{1}{2}} \Gamma(\frac{N+\kappa}{2})} \int \tau^{\frac{N+\kappa+3}{2}-1} exp \left[-\tau \left( \frac{\lambda^\ast}{2} + (B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY) \right) \right] \mathsf{d}\tau }[/math]

Here we recognize the formula to integrate the Gamma function:

[math]\displaystyle{ \mathsf{P}(B | Y, X) = \frac{\frac{\lambda^\ast}{2}^{\frac{N+\kappa}{2}} \Gamma(\frac{N+\kappa+3}{2})}{(2\pi)^\frac{3}{2} |\Omega|^{\frac{1}{2}} \Gamma(\frac{N+\kappa}{2})} \left( \frac{\lambda^\ast}{2} + (B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY) \right)^{-\frac{N+\kappa+3}{2}} }[/math]

And we now recognize a multivariate Student's t-distribution:

[math]\displaystyle{ \mathsf{P}(B | Y, X) = \frac{\Gamma(\frac{N+\kappa+3}{2})}{\Gamma(\frac{N+\kappa}{2}) \pi^\frac{3}{2} |\lambda^\ast \Omega|^{\frac{1}{2}} } \left( 1 + \frac{(B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY)}{\lambda^\ast} \right)^{-\frac{N+\kappa+3}{2}} }[/math]

We hence can write:

[math]\displaystyle{ B | Y, X \sim \mathcal{S}_{N+\kappa}(\Omega X^TY, \; (Y^T Y - Y^T X \Omega X^T Y + \lambda) \Omega) }[/math]

Bayes Factor: one way to answer our goal above ("is there an effect of the genotype on the phenotype?") is to do hypothesis testing.

We want to test the following null hypothesis:

[math]\displaystyle{ H_0: \; a = d = 0 }[/math]

In Bayesian modeling, hypothesis testing is performed with a Bayes factor, which in our case can be written as:

[math]\displaystyle{ \mathrm{BF} = \frac{\mathsf{P}(Y | X, a \neq 0, d \neq 0)}{\mathsf{P}(Y | X, a = 0, d = 0)} }[/math]

We can shorten this into:

[math]\displaystyle{ \mathrm{BF} = \frac{\mathsf{P}(Y | X)}{\mathsf{P}_0(Y)} }[/math]

Note that, compare to frequentist hypothesis testing which focuses on the null, the Bayes factor requires to explicitly model the data under the alternative. This makes a big difference when interpreting the results (see below).

Let's start with the numerator:

[math]\displaystyle{ \mathsf{P}(Y | X) = \int \mathsf{P}(\tau) \mathsf{P}(Y | X, \tau) \mathsf{d}\tau }[/math]

First, let's calculate what is inside the integral:

[math]\displaystyle{ \mathsf{P}(Y | X, \tau) = \frac{\mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B)}{\mathsf{P}(B | Y, X, \tau)} }[/math]

Using the formula obtained previously and doing some algebra gives:

[math]\displaystyle{ \mathsf{P}(Y | X, \tau) = \left( \frac{\tau}{2 \pi} \right)^{\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} exp\left( -\frac{\tau}{2} (Y^TY - Y^TX\Omega X^TY) \right) }[/math]

Now we can integrate out [math]\displaystyle{ \tau }[/math] (note the small typo in equation 9 of supplementary text S1 of Servin & Stephens):

[math]\displaystyle{ \mathsf{P}(Y | X) = (2\pi)^{-\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} \frac{\frac{\lambda}{2}^{\frac{\kappa}{2}}}{\Gamma(\frac{\kappa}{2})} \int \tau^{\frac{N+\kappa}{2}-1} exp \left( -\frac{\tau}{2} (Y^TY - Y^TX\Omega X^TY + \lambda) \right) }[/math]

Inside the integral, we recognize the almost-complete pdf of a Gamma distribution. As it has to integrate to one, we get:

[math]\displaystyle{ \mathsf{P}(Y | X) = (2\pi)^{-\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} \left( \frac{\lambda}{2} \right)^{\frac{\kappa}{2}} \frac{\Gamma(\frac{N+\kappa}{2})}{\Gamma(\frac{\kappa}{2})} \left( \frac{Y^TY - Y^TX\Omega X^TY + \lambda}{2} \right)^{-\frac{N+\kappa}{2}} }[/math]

We can use this expression also under the null. In this case, as we need neither [math]\displaystyle{ a }[/math] nor [math]\displaystyle{ d }[/math], [math]\displaystyle{ B }[/math] is simply [math]\displaystyle{ \mu }[/math], [math]\displaystyle{ \Sigma_B }[/math] is [math]\displaystyle{ \sigma_{\mu}^2 }[/math] and [math]\displaystyle{ X }[/math] is a vector of 1's. We can also defines [math]\displaystyle{ \Omega_0 = ((\sigma_{\mu}^2)^{-1} + N)^{-1} }[/math]. In the end, this gives:

[math]\displaystyle{ \mathsf{P}_0(Y) = (2\pi)^{-\frac{N}{2}} \frac{|\Omega_0|^{\frac{1}{2}}}{\sigma_{\mu}} \left( \frac{\lambda}{2} \right)^{\frac{\kappa}{2}} \frac{\Gamma(\frac{N+\kappa}{2})}{\Gamma(\frac{\kappa}{2})} \left( \frac{Y^TY - \Omega_0 N^2 \bar{Y}^2 + \lambda}{2} \right)^{-\frac{N+\kappa}{2}} }[/math]

We can therefore write the Bayes factor:

[math]\displaystyle{ \mathrm{BF} = \left( \frac{|\Omega|}{\Omega_0} \right)^{\frac{1}{2}} \frac{1}{\sigma_a \sigma_d} \left( \frac{Y^TY - Y^TX\Omega X^TY + \lambda}{Y^TY - \Omega_0 N^2 \bar{Y}^2 + \lambda} \right)^{-\frac{N+\kappa}{2}} }[/math]

When the Bayes factor is large, we say that there is enough evidence in the data to support the alternative. Indeed, the Bayesian testing procedure corresponds to measuring support for the specific alternative hypothesis compared to the null hypothesis. Importantly, note that, for a frequentist testing procedure, we would say that there is enough evidence in the data to reject the null. However we wouldn't say anything about the alternative as we don't model it.

The threshold to say that a Bayes factor is large depends on the field. It is possible to use the Bayes factor as a test statistic when doing permutation testing, and then control the false discovery rate. This can give an idea of a reasonable threshold.

Hyperparameters: the model has 5 hyperparameters, [math]\displaystyle{ \{\kappa, \, \lambda, \, \sigma_{\mu}, \, \sigma_a, \, \sigma_d\} }[/math]. How should we choose them?

Such a question is never easy to answer. But note that all hyperparameters are not that important, especially in typical quantitative genetics applications. For instance, we are mostly interested in those that determine the magnitude of the effects, [math]\displaystyle{ \sigma_a }[/math] and [math]\displaystyle{ \sigma_d }[/math], so let's deal with the others first.

As explained in Servin & Stephens, the posteriors for [math]\displaystyle{ \tau }[/math] and [math]\displaystyle{ B }[/math] change appropriately with shifts ([math]\displaystyle{ y+c }[/math]) and scaling ([math]\displaystyle{ y \times c }[/math]) in the phenotype when taking their limits. This also gives us a new Bayes factor, the one used in practice (see Guan & Stephens, 2008):

[math]\displaystyle{ \mathrm{lim}_{\sigma_{\mu} \rightarrow \infty \; ; \; \lambda \rightarrow 0 \; ; \; \kappa \rightarrow 0 } \; \mathrm{BF} = \left( \frac{N}{|\Sigma_B^{-1} + X^TX|} \right)^{\frac{1}{2}} \frac{1}{\sigma_a \sigma_d} \left( \frac{Y^TY - Y^TX (\Sigma_B^{-1} + X^TX)^{-1} X^TY}{Y^TY - N \bar{Y}^2} \right)^{-\frac{N}{2}} }[/math]

Now, for the important hyperparameters, [math]\displaystyle{ \sigma_a }[/math] and [math]\displaystyle{ \sigma_d }[/math], it is usual to specify a grid of values, i.e. [math]\displaystyle{ M }[/math] pairs [math]\displaystyle{ (\sigma_a, \sigma_d) }[/math]. For instance, Guan & Stephens used the following grid:

[math]\displaystyle{ M=4 \; ; \; \sigma_a \in \{0.05, 0.1, 0.2, 0.4\} \; ; \; \sigma_d = \frac{\sigma_a}{4} }[/math]

Then, we can average the Bayes factors obtained over the grid using, as a first approximation, equal weights:

[math]\displaystyle{ \mathrm{BF} = \sum_{m \, \in \, \text{grid}} \frac{1}{M} \, \mathrm{BF}(\sigma_a^{(m)}, \sigma_d^{(m)}) }[/math]

In eQTL studies, the weights can be estimated from the data using a hierarchical model (see below), by pooling all genes together as in Veyrieras et al (PLoS Genetics, 2010).

Implementation: the following R function is adapted from Servin & Stephens supplementary text 1.

BF <- function(G=NULL, Y=NULL, sigma.a=NULL, sigma.d=NULL, get.log10=TRUE){
  stopifnot(! is.null(G), ! is.null(Y), ! is.null(sigma.a), ! is.null(sigma.d))
  subset <- complete.cases(Y) & complete.cases(G)
  Y <- Y[subset]
  G <- G[subset]
  stopifnot(length(Y) == length(G))
  N <- length(G)
  X <- cbind(rep(1,N), G, G == 1)
  inv.Sigma.B <- diag(c(0, 1/sigma.a^2, 1/sigma.d^2))
  inv.Omega <- inv.Sigma.B + t(X) %*% X
  inv.Omega0 <- N
  tY.Y <- t(Y) %*% Y
  log10.BF <- as.numeric(0.5 * log10(inv.Omega0) -
                         0.5 * log10(det(inv.Omega)) -
                         log10(sigma.a) - log10(sigma.d) -
                         (N/2) * (log10(tY.Y - t(Y) %*% X %*% solve(inv.Omega)
                                        %*% t(X) %*% cbind(Y)) -
                                  log10(tY.Y - N*mean(Y)^2)))
  if(get.log10)
    return(log10.BF)
  else
    return(10^log10.BF)
}

In the same vein as what is explained here, we can simulate data under different scenarios and check the BFs:

N <- 300    # play with it
PVE <- 0.1  # play with it
grid <- c(0.05, 0.1, 0.2, 0.4, 0.8, 1.6, 3.2)
MAF <- 0.3
G <- rbinom(n=N, size=2, prob=MAF)
tau <- 1
a <- sqrt((2/5) * (PVE / (tau * MAF * (1-MAF) * (1-PVE))))
d <- a / 2
mu <- rnorm(n=1, mean=0, sd=10)
Y <- mu + a * G + d * (G == 1) + rnorm(n=N, mean=0, sd=tau)
for(m in 1:length(grid))
  print(BF(G, Y, grid[m], grid[m]/4))

Binary phenotype: using a similar notation, we model case-control studies with a logistic regression where the probability to be a case is [math]\displaystyle{ \mathsf{P}(y_i = 1) = p_i }[/math].

There are many equivalent ways to write the likelihood, the usual one being:

[math]\displaystyle{ y_i | p_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i) }[/math] with the log-odds (logit function) being [math]\displaystyle{ \mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1} }[/math]

Let's use [math]\displaystyle{ X_i^T=[1 \; g_i \; \mathbf{1}_{g_i=1}] }[/math] to denote the [math]\displaystyle{ i }[/math]-th row of the design matrix [math]\displaystyle{ X }[/math]. We can also keep the same definition as above for [math]\displaystyle{ B=[\mu \; a \; d]^T }[/math]. Thus we have:

[math]\displaystyle{ p_i = \frac{e^{X_i^TB}}{1 + e^{X_i^TB}} }[/math]

As the [math]\displaystyle{ y_i }[/math]'s can only take [math]\displaystyle{ 0 }[/math] and [math]\displaystyle{ 1 }[/math] as values, the likelihood can be written as:

[math]\displaystyle{ \mathcal{L}(B) = \mathsf{P}(Y | X, B) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i} }[/math]

We still use the same prior as above for [math]\displaystyle{ B }[/math] (but there is no [math]\displaystyle{ \tau }[/math] anymore), so that:

[math]\displaystyle{ B | \Sigma_B \sim \mathcal{N}_3(0, \Sigma_B) }[/math]

where [math]\displaystyle{ \Sigma_B }[/math] is a 3 x 3 matrix with [math]\displaystyle{ [\sigma_\mu^2 \; \sigma_a^2 \; \sigma_d^2] }[/math] on the diagonal and 0 elsewhere.

As above, the Bayes factor is used to compare the two models:

[math]\displaystyle{ \mathrm{BF} = \frac{\mathsf{P}(Y | X, M1)}{\mathsf{P}(Y | X, M0)} = \frac{\mathsf{P}(Y | X, a \neq 0, d \neq 0)}{\mathsf{P}(Y | X, a=0, d=0)} = \frac{\int \mathsf{P}(B) \mathsf{P}(Y | X, B) \mathrm{d}B}{\int \mathsf{P}(\mu) \mathsf{P}(Y | X, \mu) \mathrm{d}\mu} }[/math]

The interesting point here is that there is no way to analytically calculate these integrals (marginal likelihoods). Therefore, we will use Laplace's method to approximate them, as in Guan & Stephens (2008).

Starting with the numerator:

[math]\displaystyle{ \mathsf{P}(Y|X,M1) = \int \exp \left[ N \left( \frac{1}{N} \mathrm{ln} \, \mathsf{P}(B) + \frac{1}{N} \mathrm{ln} \, \mathsf{P}(Y | X, B) \right) \right] \mathsf{d}B }[/math]

[math]\displaystyle{ \mathsf{P}(Y|X,M1) = \int \exp \left\{ N \left[ \frac{1}{N} \left( \mathrm{ln} \left( (2 \pi)^{-\frac{3}{2}} \, \frac{1}{\sigma_\mu \sigma_a \sigma_d} \, \exp\left( -\frac{1}{2} (\frac{\mu^2}{\sigma_\mu^2} + \frac{a^2}{\sigma_a^2} + \frac{d^2}{\sigma_d^2}) \right) \right) \right) + \frac{1}{N} \left( \sum_{i=1}^N \left( y_i \, \mathrm{ln} (p_i) + (1-y_i) \, \mathrm{ln} (1-p_i) \right) \right) \right] \right\} \mathsf{d}B }[/math]

Let's use [math]\displaystyle{ f }[/math] to denote the function inside the exponential:

[math]\displaystyle{ \mathsf{P}(Y|X,M1) = \int \exp \left( N \; f(B) \right) \mathsf{d}B }[/math]

The function [math]\displaystyle{ f }[/math] is defined by:

[math]\displaystyle{ f: \mathbb{R}^3 \rightarrow \mathbb{R} }[/math]

[math]\displaystyle{ f(B) = \frac{1}{N} \left( -\frac{3}{2} \mathrm{ln}(2 \pi) - \frac{1}{2} \mathrm{ln}(|\Sigma_B|) - \frac{1}{2}(B^T \Sigma_B^{-1} B) \right) + \frac{1}{N} \sum_{i=1}^N \left( y_i \, X_i^T B - \mathrm{ln}(1 + e^{X_i^TB}) \right) }[/math]

This function will then be used to approximate the integral, like this:

[math]\displaystyle{ \mathsf{P}(Y|X,M1) \approx N^{-3/2} (2 \pi)^{3/2} |H(B^\star)|^{-1/2} e^{N f(B^\star)} }[/math]

where [math]\displaystyle{ H }[/math] is the Hessian of [math]\displaystyle{ f }[/math] and [math]\displaystyle{ B^\star = [\mu^\star a^\star d^\star]^T }[/math] is the point at which [math]\displaystyle{ f }[/math] is maximized.

We therefore need to find [math]\displaystyle{ B^\star }[/math]. As it maximizes [math]\displaystyle{ f }[/math], we need to calculate the first derivatives of [math]\displaystyle{ f }[/math]. Let's do this the univariate way:

[math]\displaystyle{ \frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial \beta} }[/math]

where [math]\displaystyle{ \beta }[/math] is [math]\displaystyle{ \mu }[/math], [math]\displaystyle{ a }[/math] or [math]\displaystyle{ d }[/math].

A simple form for the first derivatives of [math]\displaystyle{ p_i }[/math] also exists when writing [math]\displaystyle{ p_i = e^{X_i^tB} (1 + e^{X_i^tB})^{-1} }[/math]:

[math]\displaystyle{ \frac{\partial p_i}{\partial \beta} = \left[ e^{X_i^tB} (1 + e^{X_i^tB})^{-1} + e^{X_i^tB} \left( -e^{X_i^tB} (1 + e^{X_i^tB})^{-2} \right) \right] \frac{\partial X_i^TB}{\partial \beta} }[/math]

[math]\displaystyle{ \frac{\partial p_i}{\partial \beta} = \left[ \frac{e^{X_i^tB} (1 + e^{X_i^tB}) - (e^{X_i^tB})^2}{(1 + e^{X_i^tB})^2} \right] \frac{\partial X_i^TB}{\partial \beta} }[/math]

[math]\displaystyle{ \frac{\partial p_i}{\partial \beta} = \left[ p_i (1 - p_i) \right] \frac{\partial X_i^TB}{\partial \beta} }[/math]

where [math]\displaystyle{ \frac{\partial X_i^TB}{\partial \beta} }[/math] is equal to [math]\displaystyle{ 1, \, g_i, \, \mathbf{1}_{g_i=1} }[/math] when [math]\displaystyle{ \beta }[/math] corresponds respectively to [math]\displaystyle{ \mu, \, a, \, d }[/math].

This simplifies the first derivatives of [math]\displaystyle{ f }[/math] into:

[math]\displaystyle{ \frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial \beta} }[/math]

When setting [math]\displaystyle{ \frac{\partial f}{\partial \beta}(\beta^\star) = 0 }[/math], we observe that [math]\displaystyle{ \beta^\star }[/math] is present not only alone but also inside the sum, in the [math]\displaystyle{ p_i }[/math]'s: indeed [math]\displaystyle{ p_i }[/math] is a non-linear function of [math]\displaystyle{ B }[/math]. This means that an iterative procedure is required, typically Newton's method.

To use it, we need the second derivatives of [math]\displaystyle{ f }[/math]:

[math]\displaystyle{ \frac{\partial^2 f}{\partial \beta^2} = - \frac{1}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N \left[ (-p_i(1-p_i)\frac{\partial X_i^TB}{\partial \beta}) + (y_i-p_i)\frac{\partial^2 X_i^TB}{\partial \beta^2} \right] }[/math]

The second derivatives of [math]\displaystyle{ X_i^TB }[/math] are all equal to 0:

[math]\displaystyle{ \frac{\partial^2 f}{\partial \beta^2} = - \frac{1}{N \, \sigma_\beta^2} - \frac{1}{N} \sum_{i=1}^N p_i(1-p_i)\frac{\partial X_i^TB}{\partial \beta} }[/math]

Note that the second derivatives of [math]\displaystyle{ f }[/math] are strictly negative. Therefore, [math]\displaystyle{ f }[/math] is globally convex, which means that it has a unique global maximum, at [math]\displaystyle{ B^\star }[/math]. As a consequence, we have the right to use Laplace's method to approximate the integral around its maximum.

finding the maximums: iterative procedure, update equations or generic solver -> to do

implementation: in R -> to do

finding the effect sizes and their std error: to do

Link between Bayes factor and P-value: see Wakefield (2008)

to do

Hierarchical model: pooling genes, learn weights for grid and genomic annotations, see Veyrieras et al (PLoS Genetics, 2010)

to do

Multiple SNPs with LD: joint analysis of multiple SNPs, handle correlation between them, see Guan & Stephens (Annals of Applied Statistics, 2011) for MCMC, see Carbonetto & Stephens (Bayesian Analysis, 2012) for Variational Bayes

to do

Confounding factors in phenotype: factor analysis, see Stegle et al (PLoS Computational Biology, 2010)

to do

Genetic relatedness: linear mixed model, see Zhou & Stephens (Nature Genetics, 2012)

to do

Discrete phenotype: count data as from RNA-seq, Poisson-like likelihood, see Sun (Biometrics, 2012)

to do

Multiple phenotypes: matrix-variate distributions, tensors

to do

Non-independent genes: enrichment in known pathways, learn "modules"

to do

References:
- Servin & Stephens (PLoS Genetics, 2007)
- Guan & Stephens (PLoS Genetics, 2008)
- Stephens & Balding (Nature Reviews Genetics, 2009)

User:Timothee Flutre/Notebook/Postdoc/2011/11/10: Difference between revisions

Revision as of 15:41, 3 February 2013

Bayesian model of univariate linear regression for QTL detection

Navigation menu

Page actions

Page actions

Personal tools

Navigation

Search

research

Tools

@@ Line 9: / Line 9: @@
-''See Servin & Stephens (PLoS Genetics, 2007).''
+''This page aims at helping people like me, interested in quantitative genetics, to get a better understanding of some Bayesian models, most importantly the impact of the modeling assumptions as well as the underlying maths. It starts with a simple model, and gradually increases the scope to relax assumptions. See references to scientific articles at the end.''
-* '''Data''': let's assume that we obtained data from N individuals. We note <math>y_1,\ldots,y_N</math> the (quantitative) phenotypes (e.g. expression level at a given gene), and <math>g_1,\ldots,g_N</math> the genotypes at a given SNP (as allele dose, 0, 1 or 2).
+* '''Data''': let's assume that we obtained data from N individuals. We note <math>y_1,\ldots,y_N</math> the (quantitative) phenotypes (e.g. expression levels at a given gene), and <math>g_1,\ldots,g_N</math> the genotypes at a given SNP (encoded as allele dose: 0, 1 or 2).
@@ Line 21: / Line 21: @@
-* '''Likelihood''': <math>\forall i \in \{1,\ldots,N\}, \; y_i = \mu + \beta_1 g_i + \beta_2 \mathbf{1}_{g_i=1} + \epsilon_i \text{ with } \epsilon_i \overset{i.i.d}{\sim} \mathcal{N}(0,\tau^{-1})</math>
+* '''Likelihood''': we start by writing the usual linear regression for one individual
+<math>\forall i \in \{1,\ldots,N\}, \; y_i = \mu + \beta_1 g_i + \beta_2 \mathbf{1}_{g_i=1} + \epsilon_i \; \text{ with } \; \epsilon_i \; \overset{i.i.d}{\sim} \; \mathcal{N}(0,\tau^{-1})</math>
 where <math>\beta_1</math> is in fact the additive effect of the SNP, noted <math>a</math> from now on, and <math>\beta_2</math> is the dominance effect of the SNP, <math>d = a k</math>.
-Let's now write in matrix notation:
+Let's now write the model in matrix notation:
 <math>Y = X B + E \text{ where } B = [ \mu \; a \; d ]^T</math>
-which gives the following conditional distribution for the phenotypes:
+This gives the following [http://en.wikipedia.org/wiki/Multivariate_normal_distribution multivariate Normal distribution] for the phenotypes:
-<math>Y | X, B, \tau \sim \mathcal{N}(XB, \tau^{-1} I_N)</math>
+<math>Y | X, \tau, B \sim \mathcal{N}(XB, \tau^{-1} I_N)</math>
+Even though we can write the likelihood as a multivariate Normal, I still keep the term "univariate" in the title because the regression has a single response, <math>Y</math>.
+It is usual to keep the term "multivariate" for the case where there is a matrix of responses (i.e. multiple phenotypes).
 The likelihood of the parameters given the data is therefore:
@@ Line 37: / Line 42: @@
 <math>\mathcal{L}(\tau, B) = \mathsf{P}(Y | X, \tau, B)</math>
-<math>\mathcal{L}(\tau, B) = \left(\frac{\tau}{2 \pi}\right)^{N/2} exp \left( -\frac{\tau}{2} (Y - XB)^T (Y - XB) \right)</math>
+<math>\mathcal{L}(\tau, B) = \left(\frac{\tau}{2 \pi}\right)^{\frac{N}{2}} exp \left( -\frac{\tau}{2} (Y - XB)^T (Y - XB) \right)</math>
-* '''Priors''': we use the usual conjugate prior
+* '''Priors''': we use the usual [http://en.wikipedia.org/wiki/Conjugate_prior conjugate prior]
 <math>\mathsf{P}(\tau, B) = \mathsf{P}(\tau) \mathsf{P}(B | \tau)</math>
+A [http://en.wikipedia.org/wiki/Gamma_distribution Gamma distribution] for <math>\tau</math>:
 <math>\tau \sim \Gamma(\kappa/2, \, \lambda/2)</math>
+which means:
+<math>\mathsf{P}(\tau) = \frac{\frac{\lambda}{2}^{\kappa/2}}{\Gamma(\frac{\kappa}{2})} \tau^{\frac{\kappa}{2}-1} e^{-\frac{\lambda}{2} \tau}</math>
+And a multivariate Normal distribution for <math>B</math>:
 <math>B | \tau \sim \mathcal{N}(\vec{0}, \, \tau^{-1} \Sigma_B) \text{ with } \Sigma_B = diag(\sigma_{\mu}^2, \sigma_a^2, \sigma_d^2)</math>
+which means:
+<math>\mathsf{P}(B | \tau) = \left(\frac{\tau}{2 \pi}\right)^{\frac{3}{2}} |\Sigma_B|^{-\frac{1}{2}} exp \left(-\frac{\tau}{2} B^T \Sigma_B^{-1} B \right)</math>
-* '''Joint posterior''':
+* '''Joint posterior (1)''':
 <math>\mathsf{P}(\tau, B | Y, X) = \mathsf{P}(\tau | Y, X) \mathsf{P}(B | Y, X, \tau)</math>
@@ Line 55: / Line 72: @@
 * '''Conditional posterior of B''':
-<math>\mathsf{P}(B | Y, X, \tau) = \mathsf{P}(B, Y | X, \tau)</math>
 <math>\mathsf{P}(B | Y, X, \tau) = \frac{\mathsf{P}(B, Y | X, \tau)}{\mathsf{P}(Y | X, \tau)}</math>
-<math>\mathsf{P}(B | Y, X, \tau) = \frac{\mathsf{P}(B | \tau) \mathsf{P}(Y | X, B, \tau)}{\int \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B) \mathsf{d}B}</math>
+Let's neglect the normalization constant for now:
-Here and in the following, we neglect all constants (e.g. normalization constant, <math>Y^TY</math>, etc):
 <math>\mathsf{P}(B | Y, X, \tau) \propto \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B)</math>
-We use the prior and likelihood and keep only the terms in <math>B</math>:
+Similarly, let's keep only the terms in <math>B</math> for the moment:
 <math>\mathsf{P}(B | Y, X, \tau) \propto exp(B^T \Sigma_B^{-1} B) exp((Y-XB)^T(Y-XB))</math>
@@ Line 78: / Line 91: @@
 <math>\mathsf{P}(B | Y, X, \tau) \propto exp(B^T (\Sigma_B^{-1} + X^TX) B - Y^TXB -B^TX^TY)</math>
-Let's define <math>\Omega = (\Sigma_B^{-1} + X^TX)^{-1}</math>. We can see that <math>\Omega^T=\Omega</math>, which means that <math>\Omega</math> is a [http://en.wikipedia.org/wiki/Symmetric_matrix symmetric matrix].
+Importantly, let's define:
+<math>\Omega = (\Sigma_B^{-1} + X^TX)^{-1}</math>
+We can see that <math>\Omega^T=\Omega</math>, which means that <math>\Omega</math> is a [http://en.wikipedia.org/wiki/Symmetric_matrix symmetric matrix].
 This is particularly useful here because we can use the following equality: <math>\Omega^{-1}\Omega^T=I</math>.
@@ Line 104: / Line 121: @@
 Again, we use the priors and likelihoods specified above (but everything inside the integral is kept inside it, even if it doesn't depend on <math>B</math>!):
-<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{1/2} \tau^{N/2} exp(-\frac{\tau}{2} B^T \Sigma_B^{-1} B) exp(-\frac{\tau}{2} (Y - XB)^T (Y - XB)) \mathsf{d}B</math>
+<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{3/2} \tau^{N/2} exp(-\frac{\tau}{2} B^T \Sigma_B^{-1} B) exp(-\frac{\tau}{2} (Y - XB)^T (Y - XB)) \mathsf{d}B</math>
 As we used a conjugate prior for <math>\tau</math>, we know that we expect a Gamma distribution for the posterior.
@@ Line 110: / Line 127: @@
 We also factorize inside the exponential:
-<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{1/2} exp \left[-\frac{\tau}{2} \left( (B - \Omega X^T Y)^T \Omega^{-1} (B - \Omega X^T Y) - Y^T X \Omega X^T Y + Y^T Y \right) \right] \mathsf{d}B</math>
+<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} \int \tau^{3/2} exp \left[-\frac{\tau}{2} \left( (B - \Omega X^T Y)^T \Omega^{-1} (B - \Omega X^T Y) - Y^T X \Omega X^T Y + Y^T Y \right) \right] \mathsf{d}B</math>
 We recognize the conditional posterior of <math>B</math>.
 This allows us to use the fact that the pdf of the Normal distribution integrates to one:
-<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} exp\left[-\frac{\tau}{2} (Y^T X \Omega X^T Y + Y^T Y) \right]</math>
+<math>\mathsf{P}(\tau | Y, X) \propto \tau^{\frac{N+\kappa}{2} - 1} e^{-\frac{\lambda}{2} \tau} exp\left[-\frac{\tau}{2} (Y^T Y - Y^T X \Omega X^T Y) \right]</math>
+We finally recognize a Gamma distribution, allowing us to write the posterior as:
+<math>\tau | Y, X \sim \Gamma \left( \frac{N+\kappa}{2}, \; \frac{1}{2} (Y^T Y - Y^T X \Omega X^T Y + \lambda) \right)</math>
+* '''Joint posterior (2)''': sometimes it is said that the joint posterior follows a Normal Inverse Gamma distribution:
+<math>B, \tau | Y, X \sim \mathcal{N}IG(\Omega X^TY, \; \tau^{-1}\Omega, \; \frac{N+\kappa}{2}, \; \frac{\lambda^\ast}{2})</math>
+where <math>\lambda^\ast = Y^T Y - Y^T X \Omega X^T Y + \lambda</math>
+* '''Marginal posterior of B''': we can now integrate out <math>\tau</math>:
+<math>\mathsf{P}(B | Y, X) = \int \mathsf{P}(\tau) \mathsf{P}(B | Y, X, \tau) \mathsf{d}\tau</math>
+<math>\mathsf{P}(B | Y, X) = \frac{\frac{\lambda^\ast}{2}^{\frac{N+\kappa}{2}}}{(2\pi)^\frac{3}{2} |\Omega|^{\frac{1}{2}} \Gamma(\frac{N+\kappa}{2})} \int \tau^{\frac{N+\kappa+3}{2}-1} exp \left[-\tau \left( \frac{\lambda^\ast}{2} + (B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY) \right) \right] \mathsf{d}\tau</math>
+Here we recognize the formula to [http://en.wikipedia.org/wiki/Gamma_function#Integration_problems integrate the Gamma function]:
+<math>\mathsf{P}(B | Y, X) = \frac{\frac{\lambda^\ast}{2}^{\frac{N+\kappa}{2}} \Gamma(\frac{N+\kappa+3}{2})}{(2\pi)^\frac{3}{2} |\Omega|^{\frac{1}{2}} \Gamma(\frac{N+\kappa}{2})} \left( \frac{\lambda^\ast}{2} + (B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY) \right)^{-\frac{N+\kappa+3}{2}}</math>
+And we now recognize a [http://en.wikipedia.org/wiki/Multivariate_t-distribution multivariate Student's t-distribution]:
+<math>\mathsf{P}(B | Y, X) = \frac{\Gamma(\frac{N+\kappa+3}{2})}{\Gamma(\frac{N+\kappa}{2}) \pi^\frac{3}{2} |\lambda^\ast \Omega|^{\frac{1}{2}} } \left( 1 + \frac{(B - \Omega X^TY)^T \Omega^{-1} (B - \Omega X^TY)}{\lambda^\ast} \right)^{-\frac{N+\kappa+3}{2}}</math>
+We hence can write:
+<math>B | Y, X \sim \mathcal{S}_{N+\kappa}(\Omega X^TY, \; (Y^T Y - Y^T X \Omega X^T Y + \lambda) \Omega)</math>
+* '''Bayes Factor''': one way to answer our goal above ("is there an effect of the genotype on the phenotype?") is to do [http://en.wikipedia.org/wiki/Hypothesis_testing hypothesis testing].
+We want to test the following [http://en.wikipedia.org/wiki/Null_hypothesis null hypothesis]:
+<math>H_0: \; a = d = 0</math>
+In Bayesian modeling, hypothesis testing is performed with a [http://en.wikipedia.org/wiki/Bayes_factor Bayes factor], which in our case can be written as:
+<math>\mathrm{BF} = \frac{\mathsf{P}(Y | X, a \neq 0, d \neq 0)}{\mathsf{P}(Y | X, a = 0, d = 0)}</math>
+We can shorten this into:
+<math>\mathrm{BF} = \frac{\mathsf{P}(Y | X)}{\mathsf{P}_0(Y)}</math>
+Note that, compare to frequentist hypothesis testing which focuses on the null, the Bayes factor requires to explicitly model the data under the alternative.
+This makes a big difference when interpreting the results (see below).
+Let's start with the numerator:
+<math>\mathsf{P}(Y | X) = \int \mathsf{P}(\tau) \mathsf{P}(Y | X, \tau) \mathsf{d}\tau</math>
+First, let's calculate what is inside the integral:
+<math>\mathsf{P}(Y | X, \tau) = \frac{\mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B)}{\mathsf{P}(B | Y, X, \tau)}</math>
+Using the formula obtained previously and doing some algebra gives:
+<math>\mathsf{P}(Y | X, \tau) = \left( \frac{\tau}{2 \pi} \right)^{\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} exp\left( -\frac{\tau}{2} (Y^TY - Y^TX\Omega X^TY) \right)</math>
+Now we can integrate out <math>\tau</math> (note the small typo in equation 9 of supplementary text S1 of Servin & Stephens):
+<math>\mathsf{P}(Y | X) = (2\pi)^{-\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} \frac{\frac{\lambda}{2}^{\frac{\kappa}{2}}}{\Gamma(\frac{\kappa}{2})} \int \tau^{\frac{N+\kappa}{2}-1} exp \left( -\frac{\tau}{2} (Y^TY - Y^TX\Omega X^TY + \lambda) \right)</math>
+Inside the integral, we recognize the almost-complete pdf of a Gamma distribution.
+As it has to integrate to one, we get:
+<math>\mathsf{P}(Y | X) = (2\pi)^{-\frac{N}{2}} \left( \frac{|\Omega|}{|\Sigma_B|} \right)^{\frac{1}{2}} \left( \frac{\lambda}{2} \right)^{\frac{\kappa}{2}} \frac{\Gamma(\frac{N+\kappa}{2})}{\Gamma(\frac{\kappa}{2})} \left( \frac{Y^TY - Y^TX\Omega X^TY + \lambda}{2} \right)^{-\frac{N+\kappa}{2}}</math>
+We can use this expression also under the null.
+In this case, as we need neither <math>a</math> nor <math>d</math>, <math>B</math> is simply <math>\mu</math>, <math>\Sigma_B</math> is <math>\sigma_{\mu}^2</math> and <math>X</math> is a vector of 1's.
+We can also defines <math>\Omega_0 = ((\sigma_{\mu}^2)^{-1} + N)^{-1}</math>.
+In the end, this gives:
+<math>\mathsf{P}_0(Y) = (2\pi)^{-\frac{N}{2}} \frac{|\Omega_0|^{\frac{1}{2}}}{\sigma_{\mu}} \left( \frac{\lambda}{2} \right)^{\frac{\kappa}{2}} \frac{\Gamma(\frac{N+\kappa}{2})}{\Gamma(\frac{\kappa}{2})} \left( \frac{Y^TY - \Omega_0 N^2 \bar{Y}^2 + \lambda}{2} \right)^{-\frac{N+\kappa}{2}}</math>
+We can therefore write the Bayes factor:
+<math>\mathrm{BF} = \left( \frac{|\Omega|}{\Omega_0} \right)^{\frac{1}{2}} \frac{1}{\sigma_a \sigma_d} \left( \frac{Y^TY - Y^TX\Omega X^TY + \lambda}{Y^TY - \Omega_0 N^2 \bar{Y}^2 + \lambda} \right)^{-\frac{N+\kappa}{2}}</math>
+When the Bayes factor is large, we say that there is enough evidence in the data to ''support the alternative''.
+Indeed, the Bayesian testing procedure corresponds to measuring support for the specific alternative hypothesis compared to the null hypothesis.
+Importantly, note that, for a frequentist testing procedure, we would say that there is enough evidence in the data to ''reject the null''.
+However we wouldn't say anything about the alternative as we don't model it.
+The threshold to say that a Bayes factor is large depends on the field. It is possible to use the Bayes factor as a test statistic when doing permutation testing, and then control the false discovery rate. This can give an idea of a reasonable threshold.
+* '''Hyperparameters''': the model has 5 hyperparameters, <math>\{\kappa, \, \lambda, \, \sigma_{\mu}, \, \sigma_a, \, \sigma_d\}</math>. How should we choose them?
+Such a question is never easy to answer. But note that all hyperparameters are not that important, especially in typical quantitative genetics applications. For instance, we are mostly interested in those that determine the magnitude of the effects, <math>\sigma_a</math> and <math>\sigma_d</math>, so let's deal with the others first.
+As explained in Servin & Stephens, the posteriors for <math>\tau</math> and <math>B</math> change appropriately with shifts (<math>y+c</math>) and scaling (<math>y \times c</math>) in the phenotype when taking their limits.
+This also gives us a new Bayes factor, the one used in practice (see Guan & Stephens, 2008):
+<math>\mathrm{lim}_{\sigma_{\mu} \rightarrow \infty \; ; \; \lambda \rightarrow 0 \; ; \; \kappa \rightarrow 0 } \; \mathrm{BF} = \left( \frac{N}{|\Sigma_B^{-1} + X^TX|} \right)^{\frac{1}{2}} \frac{1}{\sigma_a \sigma_d} \left( \frac{Y^TY - Y^TX (\Sigma_B^{-1} + X^TX)^{-1} X^TY}{Y^TY - N \bar{Y}^2} \right)^{-\frac{N}{2}}</math>
+Now, for the important hyperparameters, <math>\sigma_a</math> and <math>\sigma_d</math>, it is usual to specify a grid of values, i.e. <math>M</math> pairs <math>(\sigma_a, \sigma_d)</math>. For instance, Guan & Stephens used the following grid:
+<math>M=4 \; ; \; \sigma_a \in \{0.05, 0.1, 0.2, 0.4\} \; ; \; \sigma_d = \frac{\sigma_a}{4}</math>
+Then, we can average the Bayes factors obtained over the grid using, as a first approximation,  equal weights:
+<math>\mathrm{BF} = \sum_{m \, \in \, \text{grid}} \frac{1}{M} \, \mathrm{BF}(\sigma_a^{(m)}, \sigma_d^{(m)})</math>
+In eQTL studies, the weights can be estimated from the data using a hierarchical model (see below), by pooling all genes together as in Veyrieras ''et al'' (PLoS Genetics, 2010).
+* '''Implementation''': the following R function is adapted from Servin & Stephens supplementary text 1.
+ <nowiki>
+BF <- function(G=NULL, Y=NULL, sigma.a=NULL, sigma.d=NULL, get.log10=TRUE){
+  stopifnot(! is.null(G), ! is.null(Y), ! is.null(sigma.a), ! is.null(sigma.d))
+  subset <- complete.cases(Y) & complete.cases(G)
+  Y <- Y[subset]
+  G <- G[subset]
+  stopifnot(length(Y) == length(G))
+  N <- length(G)
+  X <- cbind(rep(1,N), G, G == 1)
+  inv.Sigma.B <- diag(c(0, 1/sigma.a^2, 1/sigma.d^2))
+  inv.Omega <- inv.Sigma.B + t(X) %*% X
+  inv.Omega0 <- N
+  tY.Y <- t(Y) %*% Y
+  log10.BF <- as.numeric(0.5 * log10(inv.Omega0) -
+.5 * log10(det(inv.Omega)) -
+                         log10(sigma.a) - log10(sigma.d) -
+                         (N/2) * (log10(tY.Y - t(Y) %*% X %*% solve(inv.Omega)
+                                        %*% t(X) %*% cbind(Y)) -
+                                  log10(tY.Y - N*mean(Y)^2)))
+  if(get.log10)
+    return(log10.BF)
+  else
+    return(10^log10.BF)
+}
+</nowiki>
+In the same vein as what is explained [http://openwetware.org/wiki/User:Timothee_Flutre/Notebook/Postdoc/2011/06/28 here], we can simulate data under different scenarios and check the BFs:
+ <nowiki>
+N <- 300    # play with it
+PVE <- 0.1  # play with it
+grid <- c(0.05, 0.1, 0.2, 0.4, 0.8, 1.6, 3.2)
+MAF <- 0.3
+G <- rbinom(n=N, size=2, prob=MAF)
+tau <- 1
+a <- sqrt((2/5) * (PVE / (tau * MAF * (1-MAF) * (1-PVE))))
+d <- a / 2
+mu <- rnorm(n=1, mean=0, sd=10)
+Y <- mu + a * G + d * (G == 1) + rnorm(n=N, mean=0, sd=tau)
+for(m in 1:length(grid))
+  print(BF(G, Y, grid[m], grid[m]/4))
+</nowiki>
+* '''Binary phenotype''': using a similar notation, we model case-control studies with a [http://en.wikipedia.org/wiki/Logistic_regression logistic regression] where the probability to be a case is <math>\mathsf{P}(y_i = 1) = p_i</math>.
+There are many equivalent ways to write the likelihood, the usual one being:
+<math>y_i | p_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i)</math> with the [http://en.wikipedia.org/wiki/Log-odds log-odds] (logit function) being <math>\mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1}</math>
+Let's use <math>X_i^T=[1 \; g_i \; \mathbf{1}_{g_i=1}]</math> to denote the <math>i</math>-th row of the design matrix <math>X</math>. We can also keep the same definition as above for <math>B=[\mu \; a \; d]^T</math>. Thus we have:
+<math>p_i = \frac{e^{X_i^TB}}{1 + e^{X_i^TB}}</math>
+As the <math>y_i</math>'s can only take <math>0</math> and <math>1</math> as values, the likelihood can be written as:
+<math>\mathcal{L}(B) = \mathsf{P}(Y | X, B) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i}</math>
+We still use the same prior as above for <math>B</math> (but there is no <math>\tau</math> anymore), so that:
+<math>B | \Sigma_B \sim \mathcal{N}_3(0, \Sigma_B)</math>
+where <math>\Sigma_B</math> is a 3 x 3 matrix with <math>[\sigma_\mu^2 \; \sigma_a^2 \; \sigma_d^2]</math> on the diagonal and 0 elsewhere.
+As above, the Bayes factor is used to compare the two models:
+<math>\mathrm{BF} = \frac{\mathsf{P}(Y | X, M1)}{\mathsf{P}(Y | X, M0)} = \frac{\mathsf{P}(Y | X, a \neq 0, d \neq 0)}{\mathsf{P}(Y | X, a=0, d=0)} = \frac{\int \mathsf{P}(B) \mathsf{P}(Y | X, B) \mathrm{d}B}{\int \mathsf{P}(\mu) \mathsf{P}(Y | X, \mu) \mathrm{d}\mu}</math>
+The interesting point here is that there is no way to analytically calculate these integrals (marginal likelihoods). Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).
+Starting with the numerator:
+<math>\mathsf{P}(Y|X,M1) = \int \exp \left[ N \left( \frac{1}{N} \mathrm{ln} \, \mathsf{P}(B) + \frac{1}{N} \mathrm{ln} \, \mathsf{P}(Y | X, B) \right) \right] \mathsf{d}B</math>
+<math>\mathsf{P}(Y|X,M1) = \int \exp \left\{ N \left[ \frac{1}{N} \left(  \mathrm{ln} \left( (2 \pi)^{-\frac{3}{2}} \, \frac{1}{\sigma_\mu \sigma_a \sigma_d} \, \exp\left( -\frac{1}{2} (\frac{\mu^2}{\sigma_\mu^2} + \frac{a^2}{\sigma_a^2} + \frac{d^2}{\sigma_d^2}) \right) \right)  \right) + \frac{1}{N} \left(  \sum_{i=1}^N \left( y_i \, \mathrm{ln} (p_i) + (1-y_i) \, \mathrm{ln} (1-p_i) \right)  \right) \right] \right\} \mathsf{d}B</math>
+Let's use <math>f</math> to denote the function inside the exponential:
+<math>\mathsf{P}(Y|X,M1) = \int \exp \left( N \; f(B) \right) \mathsf{d}B</math>
+The function <math>f</math> is defined by:
+<math>f: \mathbb{R}^3 \rightarrow \mathbb{R}</math>
+<math>f(B) = \frac{1}{N} \left( -\frac{3}{2} \mathrm{ln}(2 \pi) - \frac{1}{2} \mathrm{ln}(|\Sigma_B|) - \frac{1}{2}(B^T \Sigma_B^{-1} B) \right) + \frac{1}{N} \sum_{i=1}^N \left( y_i \, X_i^T B - \mathrm{ln}(1 + e^{X_i^TB}) \right)</math>
+This function will then be used to approximate the integral, like this:
+<math>\mathsf{P}(Y|X,M1) \approx N^{-3/2} (2 \pi)^{3/2} |H(B^\star)|^{-1/2} e^{N f(B^\star)}</math>
+where <math>H</math> is the [http://en.wikipedia.org/wiki/Hessian_matrix Hessian] of <math>f</math> and <math>B^\star = [\mu^\star a^\star d^\star]^T</math> is the point at which <math>f</math> is maximized.
+We therefore need to find <math>B^\star</math>. As it maximizes <math>f</math>, we need to calculate the first derivatives of <math>f</math>. Let's do this the univariate way:
+<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial \beta}</math>
+where <math>\beta</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.
+A simple form for the first derivatives of <math>p_i</math> also exists when writing <math>p_i = e^{X_i^tB} (1 + e^{X_i^tB})^{-1}</math>:
+<math>\frac{\partial p_i}{\partial \beta} = \left[ e^{X_i^tB} (1 + e^{X_i^tB})^{-1} + e^{X_i^tB} \left( -e^{X_i^tB} (1 + e^{X_i^tB})^{-2} \right) \right] \frac{\partial X_i^TB}{\partial \beta}</math>
+<math>\frac{\partial p_i}{\partial \beta} = \left[ \frac{e^{X_i^tB} (1 + e^{X_i^tB}) - (e^{X_i^tB})^2}{(1 + e^{X_i^tB})^2} \right] \frac{\partial X_i^TB}{\partial \beta}</math>
+<math>\frac{\partial p_i}{\partial \beta} = \left[ p_i (1 - p_i) \right] \frac{\partial X_i^TB}{\partial \beta}</math>
+where <math>\frac{\partial X_i^TB}{\partial \beta}</math> is equal to <math>1, \, g_i, \, \mathbf{1}_{g_i=1}</math> when <math>\beta</math> corresponds respectively to <math>\mu, \, a, \, d</math>.
+This simplifies the first derivatives of <math>f</math> into:
+<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial \beta}</math>
+When setting <math>\frac{\partial f}{\partial \beta}(\beta^\star) = 0</math>, we observe that <math>\beta^\star</math> is present not only alone but also inside the sum, in the <math>p_i</math>'s: indeed <math>p_i</math> is a non-linear function of <math>B</math>. This means that an iterative procedure is required, typically [http://en.wikipedia.org/wiki/Newton_method_in_optimization Newton's method].
+To use it, we need the second derivatives of <math>f</math>:
+<math>\frac{\partial^2 f}{\partial \beta^2} = - \frac{1}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N \left[ (-p_i(1-p_i)\frac{\partial X_i^TB}{\partial \beta}) + (y_i-p_i)\frac{\partial^2 X_i^TB}{\partial \beta^2} \right]</math>
+The second derivatives of <math>X_i^TB</math> are all equal to 0:
+<math>\frac{\partial^2 f}{\partial \beta^2} = - \frac{1}{N \, \sigma_\beta^2} - \frac{1}{N} \sum_{i=1}^N p_i(1-p_i)\frac{\partial X_i^TB}{\partial \beta}</math>
+Note that the second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star</math>. As a consequence, we have the right to use Laplace's method to approximate the integral around its maximum.
+finding the maximums: iterative procedure, update equations or generic solver -> to do
+implementation: in R -> to do
+finding the effect sizes and their std error: to do
+* '''Link between Bayes factor and P-value''': see Wakefield (2008)
+to do
+* '''Hierarchical model''': pooling genes, learn weights for grid and genomic annotations, see Veyrieras ''et al'' (PLoS Genetics, 2010)
+to do
+* '''Multiple SNPs with LD''': joint analysis of multiple SNPs, handle correlation between them, see Guan & Stephens (Annals of Applied Statistics, 2011) for MCMC, see Carbonetto & Stephens (Bayesian Analysis, 2012) for Variational Bayes
+to do
+* '''Confounding factors in phenotype''': factor analysis, see Stegle ''et al'' (PLoS Computational Biology, 2010)
+to do
+* '''Genetic relatedness''': linear mixed model, see Zhou & Stephens (Nature Genetics, 2012)
+to do
+* '''Discrete phenotype''': count data as from RNA-seq, Poisson-like likelihood, see Sun (Biometrics, 2012)
+to do
+* '''Multiple phenotypes''': matrix-variate distributions, tensors
+to do
+* '''Non-independent genes''': enrichment in known pathways, learn "modules"
+to do
-We finally recognize the following Gamma distribution:
-<math>\tau | Y, X \sim \Gamma \left( \frac{N+\kappa}{2}, \; \frac{1}{2} (Y^T X \Omega X^T Y + Y^T Y + \lambda) \right)</math>
+* '''References''':
+** Servin & Stephens (PLoS Genetics, 2007)
+** Guan & Stephens (PLoS Genetics, 2008)
+** Stephens & Balding (Nature Reviews Genetics, 2009)
 <!-- ##### DO NOT edit below this line unless you know what you are doing. ##### -->