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

−

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

+

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

−

* '''Goal''': we want (i) to assess the evidence in the data for an effect of the genotype on the phenotype, and (ii) estimate the posterior distribution of this effect.

+

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

Line 21:

Line 21:

−

* '''Likelihood''':

+

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

This now becomes easy to factorizes totally (remember to add the term in <math>Y</math> to [https://learnbayes.org/index.php?option=com_content&view=article&id=77:completesquare&catid=83:reference&Itemid=479 complete the square]):

In this case, as we need neither <math>a</math> nor <math>d</math>, <math>B</math> is simply <math>\mu</math>, <math>\Sigma_B</math> is <math>\sigma_{\mu}^2</math> and <math>X</math> is a vector of 1's.

+

We can also defines <math>\Omega_0 = ((\sigma_{\mu}^2)^{-1} + N)^{-1}</math>.

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.

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

+

+

As explained in Servin & Stephens, the posteriors for <math>\tau</math> and <math>B</math> change appropriately with shifts (<math>y+c</math>) and scaling (<math>y \times c</math>) in the phenotype when taking their limits.

+

This also gives us a new Bayes factor, the one used in practice (see Guan & Stephens, 2008):

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

In the same vein as what is explained [http://openwetware.org/wiki/User:Timothee_Flutre/Notebook/Postdoc/2011/06/28 here], we can simulate data under different scenarios and check the BFs:

+

+

<nowiki>

+

N <- 300 # play with it

+

PVE <- 0.1 # play with it

+

grid <- c(0.05, 0.1, 0.2, 0.4, 0.8, 1.6, 3.2)

+

MAF <- 0.3

+

G <- rbinom(n=N, size=2, prob=MAF)

+

tau <- 1

+

a <- sqrt((2/5) * (PVE / (tau * MAF * (1-MAF) * (1-PVE))))

+

d <- a / 2

+

mu <- rnorm(n=1, mean=0, sd=10)

+

Y <- mu + a * G + d * (G == 1) + rnorm(n=N, mean=0, sd=tau)

+

for(m in 1:length(grid))

+

print(BF(G, Y, grid[m], grid[m]/4))

+

</nowiki>

+

+

+

* '''Binary phenotype''': using a similar notation, we model case-control studies with a [http://en.wikipedia.org/wiki/Logistic_regression logistic regression] where the probability to be a case is <math>\mathsf{P}(y_i = 1) = p_i</math>.

+

+

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

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:

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

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

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:

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.

* '''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, <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.

Latest revision as of 20:05, 26 September 2017

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle y_1,\ldots,y_N}
the (quantitative) phenotypes (e.g. expression levels at a given gene), and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle 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

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

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{\displaystyle Y}.
It is usual to keep the term "multivariate" for the case where there is a matrix of responses (i.e. multiple phenotypes).

We can see that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \Omega^T=\Omega}
, which means that Ω{\displaystyle \Omega } is a symmetric matrix.
This is particularly useful here because we can use the following equality: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \Omega^{-1}\Omega^T=I}
.

This now becomes easy to factorizes totally (remember to add the term in Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle Y}
to complete the square):

But now, to handle the second term, we need to integrate over Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle B}
, thus effectively taking into account the uncertainty in Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle 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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle B}
!):

As we used a conjugate prior for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \tau}
, we know that we expect a Gamma distribution for the posterior.
Therefore, we can take Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \tau^{N/2}}
out of the integral and start guessing what looks like a Gamma distribution.
We also factorize inside the exponential:

We recognize the conditional posterior of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle B}
.
This allows us to use the fact that the pdf of the Normal distribution integrates to one:

Marginal posterior of B: we can now integrate out Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://api.formulasearchengine.com/v1/":): {\displaystyle \tau}
:

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

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

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.

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{\displaystyle \sigma _{a}} and σd{\displaystyle \sigma _{d}}, so let's deal with the others first.

As explained in Servin & Stephens, the posteriors for τ{\displaystyle \tau } and B{\displaystyle B} change appropriately with shifts (y+c{\displaystyle y+c}) and scaling (y×c{\displaystyle 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):

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

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

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

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

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

To find B⋆{\displaystyle B^{\star }}, we set ∇f(B⋆)=0{\displaystyle \nabla f(B^{\star })=0}. However, in this equation, B⋆{\displaystyle B^{\star }} is present not only alone but also in the pi{\displaystyle p_{i}}'s. As pi{\displaystyle p_{i}} is a non-linear function of B{\displaystyle 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{\displaystyle f} and ∇f{\displaystyle \nabla f} while the latter also requires H{\displaystyle H}. Remember that, in any case, we need H{\displaystyle H} for the Laplace approximation, so let's calculate it:

where W{\displaystyle W} is the N x N matrix with pi(1−pi){\displaystyle p_{i}(1-p_{i})} on the diagonal.

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

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{\displaystyle y_{i}=\mu +\beta g_{i}+\epsilon _{i}} assumes that the errors are uncorrelated between individuals: ϵi∼i.i.dN(0,τ−1){\displaystyle \epsilon _{i}{\overset {\mathrm {i.i.d} }{\sim }}{\mathcal {N}}(0,\tau ^{-1})}. If this is not the case, i.e. if the yi{\displaystyle y_{i}}'s are correlated but this correlation has nothing to do with the gi{\displaystyle g_{i}}'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 β{\displaystyle \beta }.

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×N{\displaystyle G\times N} matrix Y1{\displaystyle Y_{1}} 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.