User:Timothee Flutre/Notebook/Postdoc/2011/11/10

From OpenWetWare

(Difference between revisions)
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|----------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| ( improve R code for simulations)      |                                                                                                                                                                                            | Current revision (21:46, 27 February 2013) (view source) ( add info about confounders in phenotype) |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| (10 intermediate revisions not shown.) |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| Line 23:                               |                                                                                                                                                                                            | Line 23:                                                                                            |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        | * '''Likelihood''': we start by writing the usual linear regression for one individual                                                                                                     |                                                                                                     | * '''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>    | +                                                                                                   | <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>. |                                                                                                     | 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>.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
| Line 236:                              |                                                                                                                                                                                            | Line 236:                                                                                           |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        | <math>\mathrm{BF} = \sum_{m \, \in \, \text{grid}} \frac{1}{M} \, \mathrm{BF}(\sigma_a^{(m)}, \sigma_d^{(m)})</math>                                                                       |                                                                                                     | <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).                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| Line 284:                              |                                                                                                                                                                                            | Line 286:                                                                                           |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Binary phenotype''': case-control like in GWAS, logistic regression, see Guan & Stephens (2008) for Laplace approximation                                                             | +                                                                                                   | * '''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>.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | to do                                                                                                                                                                                      | +                                                                                                   | There are many equivalent ways to write the likelihood, the usual one being:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>y_i | p_i \; \overset{i.i.d}{\sim} \; \mathrm{Binomial}(1,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>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Link between Bayes factor and P-value''': see Wakeley (2008)                                                                                                                          | +                                                                                                   | 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>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | 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 \mathrm{ln}(p_i) + (1-y_i) \mathrm{ln}(1 - p_i) \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 two things: <math>H</math> and <math>B^\star</math>. Note that for both we need to calculate the first derivatives of <math>f</math>. As <math>f</math> is multi-dimensional (it takes values in <math>\mathbb{R}^3</math>), we need to calculate its [http://en.wikipedia.org/wiki/Gradient gradient].                                                                                                                                                                                                                                                                                                                                                                                                                                                        |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | In the following, some formulas from [http://en.wikipedia.org/wiki/Matrix_calculus matrix calculus] are sometimes required. In such cases, I will use the [http://en.wikipedia.org/wiki/Matrix_calculus#Layout_conventions numerator layout].                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\nabla f = \frac{\partial f}{\partial B} = \left( \frac{\partial f}{\partial \mu} \; \frac{\partial f}{\partial a} \; \frac{\partial f}{\partial d} \right)</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\nabla f = - \frac{1}{2N} \frac{\partial B^T\Sigma_B^{-1}B}{\partial B} + \frac{1}{N} \sum_i \left( y_i \frac{\partial \mathrm{ln}(p_i)}{\partial B} + (1-y_i) \frac{\partial \mathrm{ln}(1-p_i)}{\partial B} \right)</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\nabla f = - \frac{1}{N} B^T\Sigma_B^{-1} + \frac{1}{N} \sum_i \left( \frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial B}</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 B} = \frac{\partial e^{X_i^TB}}{\partial B} (1 + e^{X_i^TB})^{-1} + e^{X_i^TB} \frac{\partial (1 + e^{X_i^TB})^{-1}}{\partial B}</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\frac{\partial p_i}{\partial B} = e^{X_i^TB} \frac{\partial X_i^TB}{\partial B} (1 + e^{X_i^TB})^{-1} - e^{X_i^TB} (1 + e^{X_i^TB})^{-2} \frac{\partial (1 + e^{X_i^TB})}{\partial B}</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\frac{\partial p_i}{\partial B} = p_i X_i^T - p_i (1 + e^{X_i^TB})^{-1} e^{X_i^TB} \frac{\partial X_i^TB}{\partial B}</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\frac{\partial p_i}{\partial B} = p_i (1 - p_i) X_i^T</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | This simplifies the gradient of <math>f</math> into:                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>\nabla f = - \frac{1}{N} B^T\Sigma_B^{-1} + \frac{1}{N} \sum_i (y_i - p_i) X_i^T</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | To find <math>B^\star</math>, we set <math>\nabla f(B^\star) = 0</math>. However, in this equation, <math>B^\star</math> is present not only alone but also in the <math>p_i</math>'s. As <math>p_i</math> is a non-linear function of <math>B</math>, the equation can't be solved directly but an iterative procedure is required, typically a [http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method conjugate gradient method] (as in Guan & Stephens) or [http://en.wikipedia.org/wiki/Newton_method_in_optimization Newton's method]. The former only requires <math>f</math> and <math>\nabla f</math> while the latter also requires <math>H</math>. Remember that, in any case, we need <math>H</math> for the Laplace approximation, so let's calculate it: |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>H = - \frac{1}{N} \Sigma_B^{-1} - \frac{1}{N} \sum_i \frac{\partial p_i}{\partial B} X_i^T</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>H = - \frac{1}{N} \Sigma_B^{-1} - \frac{1}{N} \sum_i X_i^T p_i (1-p_i) X_i</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | <math>H = - \frac{1}{N} (\Sigma_B^{-1} + X^T W X)</math>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | where <math>W</math> is the N x N matrix with <math>p_i(1-p_i)</math> on the diagonal.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | Note that all 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 of <math>f</math> around its maximum.                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | 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                                                                                                                                                                                      |                                                                                                     | to do                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |
| Line 304:                              |                                                                                                                                                                                            | Line 388:                                                                                           |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Confounding factors in phenotype''': factor analysis, see Stegle ''et al'' (PLoS Computational Biology, 2010)                                                                         | +                                                                                                   | * '''Confounders in phenotype''': it is well known in molecular biology that any experiment involving several assays (e.g. measuring gene expression levels with a DNA microarray) suffers from "unknown confounders", the most (in)famous being the so-called "batch effects".                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | to do                                                                                                                                                                                      | +                                                                                                   | For instance, samples from individual 1 and 2 are correlated with each other because they were treated another day than all the other samples. Such a correlations has nothing to do with the genotype at a given SNP (in most cases). However, the core model, <math>y_i = \mu + \beta g_i + \epsilon_i</math> assumes that the errors are uncorrelated between individuals: <math>\epsilon_i \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(0,\tau^{-1})</math>. If this is not the case, i.e. if the <math>y_i</math>'s are correlated but this correlation has nothing to do with the <math>g_i</math>'s, then more variance in the errors won't be accounted for, and we'll loose power when trying to detect weak, yet non-zero <math>\beta</math>.                            |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | An intuitive way of removing these confounders is to realize that we can use all gene expression levels to try to identify them. Indeed, batch effects are very likely to affect all genes in a sample (though possibly at different magnitudes). As the effect of the confounders are, as a first approximation, typically much bigger than the effect of a SNP genotype, we can try to learn the confounders using all gene expression levels, and only them. So let's put all of them into a <math>G \times N</math> matrix <math>Y_1</math> with genes in rows and individuals in columns.                                                                                                                                                                                   |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | For the moment, the data are expressed in the [http://en.wikipedia.org/wiki/Standard_basis standard basis], i.e. the basis of the observations. But some confounders are present in these data, they contribute with noise and redundancy and hence dilute the signal. The idea is, first, to identify a new basis which will correspond to a "mix" of the original samples (e.g. one component of this "mix" may correspond to the day at which the samples were processed), and second, to remove these components from the data in order to only keep the signal.                                                                                                                                                                                                             |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | to be continued                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            | +                                                                                                   | see also factor analysis, see Stegle ''et al'' (PLoS Computational Biology, 2010)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Genetic relatedness''': linear mixed model                                                                                                                                            | +                                                                                                   | * '''Confounders in genotype''': mainly pop structure and genetic relatedness, linear mixed model (LMM), see Zhou & Stephens (Nature Genetics, 2012)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        | to do                                                                                                                                                                                      |                                                                                                     | to do                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Discrete phenotype''': count data as from RNA-seq, Poisson-like likelihood                                                                                                            | +                                                                                                   | * '''Discrete phenotype''': count data (e.g. from RNA-seq), Poisson-like likelihood, generalized linear model (GLM), see Sun (Biometrics, 2012)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        | to do                                                                                                                                                                                      |                                                                                                     | to do                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |
| Line 324:                              |                                                                                                                                                                                            | Line 416:                                                                                           |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| -                                      | * '''Non-independent genes''': enrichment in known pathways, learn "modules"                                                                                                               | +                                                                                                   | * '''Non-independent genes''': enrichment in known pathways, learn "modules", distributions on networks                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
|                                        |                                                                                                                                                                                            |                                                                                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
|                                        | to do                                                                                                                                                                                      |                                                                                                     | to do                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            |

Current revision

Project name Main project page
Previous entry      Next entry

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 $y_1,\ldots,y_N$ the (quantitative) phenotypes (e.g. expression levels at a given gene), and $g_1,\ldots,g_N$ 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

$\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})$

where β1 is in fact the additive effect of the SNP, noted a from now on, and β2 is the dominance effect of the SNP, d = ak.

Let's now write the model in matrix notation:

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

This gives the following multivariate Normal distribution for the phenotypes:

$Y | X, \tau, B \sim \mathcal{N}(XB, \tau^{-1} I_N)$

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, Y. 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:

$\mathcal{L}(\tau, B) = \mathsf{P}(Y | X, \tau, B)$

$\mathcal{L}(\tau, B) = \left(\frac{\tau}{2 \pi}\right)^{\frac{N}{2}} exp \left( -\frac{\tau}{2} (Y - XB)^T (Y - XB) \right)$

• Priors: we use the usual conjugate prior

$\mathsf{P}(\tau, B) = \mathsf{P}(\tau) \mathsf{P}(B | \tau)$

A Gamma distribution for τ:

$\tau \sim \Gamma(\kappa/2, \, \lambda/2)$

which means:

$\mathsf{P}(\tau) = \frac{\frac{\lambda}{2}^{\kappa/2}}{\Gamma(\frac{\kappa}{2})} \tau^{\frac{\kappa}{2}-1} e^{-\frac{\lambda}{2} \tau}$

And a multivariate Normal distribution for B:

$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)$

which means:

$\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)$

• Joint posterior (1):

$\mathsf{P}(\tau, B | Y, X) = \mathsf{P}(\tau | Y, X) \mathsf{P}(B | Y, X, \tau)$

• Conditional posterior of B:

$\mathsf{P}(B | Y, X, \tau) = \frac{\mathsf{P}(B, Y | X, \tau)}{\mathsf{P}(Y | X, \tau)}$

Let's neglect the normalization constant for now:

$\mathsf{P}(B | Y, X, \tau) \propto \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B)$

Similarly, let's keep only the terms in B for the moment:

$\mathsf{P}(B | Y, X, \tau) \propto exp(B^T \Sigma_B^{-1} B) exp((Y-XB)^T(Y-XB))$

We expand:

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

We factorize some terms:

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

Importantly, let's define:

$\Omega = (\Sigma_B^{-1} + X^TX)^{-1}$

We can see that ΩT = Ω, which means that Ω is a symmetric matrix. This is particularly useful here because we can use the following equality: Ω − 1ΩT = I.

$\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)$

This now becomes easy to factorizes totally:

$\mathsf{P}(B | Y, X, \tau) \propto exp((B^T - \Omega X^TY)^T\Omega^{-1}(B - \Omega X^TY))$

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

$B | Y, X, \tau \sim \mathcal{N}(\Omega X^TY, \tau^{-1} \Omega)$

• Posterior of τ:

Similarly to the equations above:

$\mathsf{P}(\tau | Y, X) \propto \mathsf{P}(\tau) \mathsf{P}(Y | X, \tau)$

But now, to handle the second term, we need to integrate over B, thus effectively taking into account the uncertainty in B:

$\mathsf{P}(\tau | Y, X) \propto \mathsf{P}(\tau) \int \mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B) \mathsf{d}B$

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 B!):

$\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$

As we used a conjugate prior for τ, we know that we expect a Gamma distribution for the posterior. Therefore, we can take τN / 2 out of the integral and start guessing what looks like a Gamma distribution. We also factorize inside the exponential:

$\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$

We recognize the conditional posterior of B. This allows us to use the fact that the pdf of the Normal distribution integrates to one:

$\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]$

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

$\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)$

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

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

where $\lambda^\ast = Y^T Y - Y^T X \Omega X^T Y + \lambda$

• Marginal posterior of B: we can now integrate out τ:

$\mathsf{P}(B | Y, X) = \int \mathsf{P}(\tau) \mathsf{P}(B | Y, X, \tau) \mathsf{d}\tau$

$\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$

Here we recognize the formula to integrate the Gamma function:

$\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}}$

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

$\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}}$

We hence can write:

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

• 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:

$H_0: \; a = d = 0$

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

$\mathrm{BF} = \frac{\mathsf{P}(Y | X, a \neq 0, d \neq 0)}{\mathsf{P}(Y | X, a = 0, d = 0)}$

We can shorten this into:

$\mathrm{BF} = \frac{\mathsf{P}(Y | X)}{\mathsf{P}_0(Y)}$

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).

$\mathsf{P}(Y | X) = \int \mathsf{P}(\tau) \mathsf{P}(Y | X, \tau) \mathsf{d}\tau$

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

$\mathsf{P}(Y | X, \tau) = \frac{\mathsf{P}(B | \tau) \mathsf{P}(Y | X, \tau, B)}{\mathsf{P}(B | Y, X, \tau)}$

Using the formula obtained previously and doing some algebra gives:

$\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)$

Now we can integrate out τ (note the small typo in equation 9 of supplementary text S1 of Servin & Stephens):

$\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)$

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

$\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}}$

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

$\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}}$

We can therefore write the Bayes factor:

$\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}}$

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, $\{\kappa, \, \lambda, \, \sigma_{\mu}, \, \sigma_a, \, \sigma_d\}$. 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, σa and σd, so let's deal with the others first.

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

$\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}}$

Now, for the important hyperparameters, σa and σd, it is usual to specify a grid of values, i.e. M pairs (σa,σd). For instance, Guan & Stephens used the following grid:

$M=4 \; ; \; \sigma_a \in \{0.05, 0.1, 0.2, 0.4\} \; ; \; \sigma_d = \frac{\sigma_a}{4}$

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

$\mathrm{BF} = \sum_{m \, \in \, \text{grid}} \frac{1}{M} \, \mathrm{BF}(\sigma_a^{(m)}, \sigma_d^{(m)})$

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 $\mathsf{P}(y_i = 1) = p_i$.

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

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

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

$p_i = \frac{e^{X_i^TB}}{1 + e^{X_i^TB}}$

As the yi's can only take 0 and 1 as values, the likelihood can be written as:

$\mathcal{L}(B) = \mathsf{P}(Y | X, B) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i}$

We still use the same prior as above for B (but there is no τ anymore), so that:

$B | \Sigma_B \sim \mathcal{N}_3(0, \Sigma_B)$

where ΣB is a 3 x 3 matrix with $(\sigma_\mu^2 \; \sigma_a^2 \; \sigma_d^2)$ on the diagonal and 0 elsewhere.

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

$\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}$

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:

$\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$

Let's use f to denote the function inside the exponential:

$\mathsf{P}(Y|X,M1) = \int \exp \left( N \; f(B) \right) \mathsf{d}B$

The function f is defined by:

$f: \mathbb{R}^3 \rightarrow \mathbb{R}$

$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 \mathrm{ln}(p_i) + (1-y_i) \mathrm{ln}(1 - p_i) \right)$

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

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

where H is the Hessian of f and $B^\star = (\mu^\star a^\star d^\star)^T$ is the point at which f is maximized.

We therefore need two things: H and $B^\star$. Note that for both we need to calculate the first derivatives of f. As f is multi-dimensional (it takes values in $\mathbb{R}^3$), we need to calculate its gradient.

In the following, some formulas from matrix calculus are sometimes required. In such cases, I will use the numerator layout.

$\nabla f = \frac{\partial f}{\partial B} = \left( \frac{\partial f}{\partial \mu} \; \frac{\partial f}{\partial a} \; \frac{\partial f}{\partial d} \right)$

$\nabla f = - \frac{1}{2N} \frac{\partial B^T\Sigma_B^{-1}B}{\partial B} + \frac{1}{N} \sum_i \left( y_i \frac{\partial \mathrm{ln}(p_i)}{\partial B} + (1-y_i) \frac{\partial \mathrm{ln}(1-p_i)}{\partial B} \right)$

$\nabla f = - \frac{1}{N} B^T\Sigma_B^{-1} + \frac{1}{N} \sum_i \left( \frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial B}$

A simple form for the first derivatives of pi also exists when writing $p_i = e^{X_i^tB} (1 + e^{X_i^tB})^{-1}$:

$\frac{\partial p_i}{\partial B} = \frac{\partial e^{X_i^TB}}{\partial B} (1 + e^{X_i^TB})^{-1} + e^{X_i^TB} \frac{\partial (1 + e^{X_i^TB})^{-1}}{\partial B}$

$\frac{\partial p_i}{\partial B} = e^{X_i^TB} \frac{\partial X_i^TB}{\partial B} (1 + e^{X_i^TB})^{-1} - e^{X_i^TB} (1 + e^{X_i^TB})^{-2} \frac{\partial (1 + e^{X_i^TB})}{\partial B}$

$\frac{\partial p_i}{\partial B} = p_i X_i^T - p_i (1 + e^{X_i^TB})^{-1} e^{X_i^TB} \frac{\partial X_i^TB}{\partial B}$

$\frac{\partial p_i}{\partial B} = p_i (1 - p_i) X_i^T$

This simplifies the gradient of f into:

$\nabla f = - \frac{1}{N} B^T\Sigma_B^{-1} + \frac{1}{N} \sum_i (y_i - p_i) X_i^T$

To find $B^\star$, we set $\nabla f(B^\star) = 0$. However, in this equation, $B^\star$ is present not only alone but also in the pi's. As pi is a non-linear function of B, the equation can't be solved directly but an iterative procedure is required, typically a conjugate gradient method (as in Guan & Stephens) or Newton's method. The former only requires f and $\nabla f$ while the latter also requires H. Remember that, in any case, we need H for the Laplace approximation, so let's calculate it:

$H = - \frac{1}{N} \Sigma_B^{-1} - \frac{1}{N} \sum_i \frac{\partial p_i}{\partial B} X_i^T$

$H = - \frac{1}{N} \Sigma_B^{-1} - \frac{1}{N} \sum_i X_i^T p_i (1-p_i) X_i$

$H = - \frac{1}{N} (\Sigma_B^{-1} + X^T W X)$

where W is the N x N matrix with pi(1 − pi) on the diagonal.

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

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

• Confounders in phenotype: it is well known in molecular biology that any experiment involving several assays (e.g. measuring gene expression levels with a DNA microarray) suffers from "unknown confounders", the most (in)famous being the so-called "batch effects".

For instance, samples from individual 1 and 2 are correlated with each other because they were treated another day than all the other samples. Such a correlations has nothing to do with the genotype at a given SNP (in most cases). However, the core model, yi = μ + βgi + εi assumes that the errors are uncorrelated between individuals: $\epsilon_i \overset{\mathrm{i.i.d}}{\sim} \mathcal{N}(0,\tau^{-1})$. If this is not the case, i.e. if the yi's are correlated but this correlation has nothing to do with the gi's, then more variance in the errors won't be accounted for, and we'll loose power when trying to detect weak, yet non-zero β.

An intuitive way of removing these confounders is to realize that we can use all gene expression levels to try to identify them. Indeed, batch effects are very likely to affect all genes in a sample (though possibly at different magnitudes). As the effect of the confounders are, as a first approximation, typically much bigger than the effect of a SNP genotype, we can try to learn the confounders using all gene expression levels, and only them. So let's put all of them into a $G \times N$ matrix Y1 with genes in rows and individuals in columns.

For the moment, the data are expressed in the standard basis, i.e. the basis of the observations. But some confounders are present in these data, they contribute with noise and redundancy and hence dilute the signal. The idea is, first, to identify a new basis which will correspond to a "mix" of the original samples (e.g. one component of this "mix" may correspond to the day at which the samples were processed), and second, to remove these components from the data in order to only keep the signal.

to be continued

see also factor analysis, see Stegle et al (PLoS Computational Biology, 2010)

• Confounders in genotype: mainly pop structure and genetic relatedness, linear mixed model (LMM), see Zhou & Stephens (Nature Genetics, 2012)

to do

• Discrete phenotype: count data (e.g. from RNA-seq), Poisson-like likelihood, generalized linear model (GLM), see Sun (Biometrics, 2012)

to do

• Multiple phenotypes: matrix-variate distributions, tensors

to do

• Non-independent genes: enrichment in known pathways, learn "modules", distributions on networks

to do

• References:
• Servin & Stephens (PLoS Genetics, 2007)
• Guan & Stephens (PLoS Genetics, 2008)
• Stephens & Balding (Nature Reviews Genetics, 2009)
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 89, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": false, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.3, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "github"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.8363431096076965, "perplexity_flag": "middle"}