User:Timothee Flutre/Notebook/Postdoc/2011/11/10 - Revision history

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add info about confounders in phenotype

2013-02-28T01:46:28Z

Bayesian model of univariate linear regression for QTL detection: add info about confounders in phenotype

←Older revision		Revision as of 01:46, 28 February 2013
Line 388:		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, see Zhou & Stephens (Nature Genetics, 2012)	+	* '''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, see Sun (Biometrics, 2012)	+	* '''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 408:		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

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add multivar eq for logistic reg

2013-02-22T02:44:01Z

Bayesian model of univariate linear regression for QTL detection: add multivar eq for logistic reg

Show changes

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */

2013-02-03T22:41:10Z

Bayesian model of univariate linear regression for QTL detection

←Older revision		Revision as of 22:41, 3 February 2013
Line 354:		Line 354:
	<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial \beta}</math>		<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial \beta}</math>

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

	To use it, we need the second derivatives of <math>f</math>:		To use it, we need the second derivatives of <math>f</math>:

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ use beta instead of x

2013-02-03T22:39:57Z

Bayesian model of univariate linear regression for QTL detection: use beta instead of x

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ use N for Laplace's method

2013-01-27T19:12:40Z

Bayesian model of univariate linear regression for QTL detection: use N for Laplace's method

←Older revision		Revision as of 19:12, 27 January 2013
Line 310:		Line 310:
	<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>		<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 calculate these integrals ~~analytically~~. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).	+	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:		Starting with the numerator:

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

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

	Let's use <math>f</math> to denote the function inside the exponential:		Let's use <math>f</math> to denote the function inside the exponential:

-	<math>\mathsf{P}(Y\|X,M1) = \int \exp \left( f(B) \right) \mathsf{d}B</math>	+	<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:		The function <math>f</math> is defined by:
Line 326:		Line 326:
	<math>f: \mathbb{R}^3 \rightarrow \mathbb{R}</math>		<math>f: \mathbb{R}^3 \rightarrow \mathbb{R}</math>

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

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

-	<math>\mathsf{P}(Y\|X) \approx (2 \pi)^{-3/2} \|H(B^\star)\|^{-1/2} ~~\exp~~^{f(B^\star)}</math>	+	<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.		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.

-	~~First~~ we need to calculate the first derivatives of <math>f</math>~~, noting that they all have a very similar form~~:	+	We therefore need to find <math>B^\star</math>. As it maximizes <math>f</math>, we need to calculate the first derivatives of <math>f</math>. Let's do this the univariate way:

-	<math>\frac{\partial f}{\partial x} = - \frac{x}{\sigma_x^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>	+	<math>\frac{\partial f}{\partial x} = - \frac{x}{N \, \sigma_x^2} + \frac{1}{N} \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>

	where <math>x</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.		where <math>x</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.

-	~~Also, there is a~~ simple form for the first derivatives of <math>p_i</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 x} = \~~frac{ (\frac{\partial X_i^TB}{\partial x}~~ e^{X_i^TB}) (1+e^{X_i^TB}) - e^{X_i^TB} (~~\frac{\partial X_i^TB}{\partial x}~~ e^{X_i^TB}~~)}{~~(1+e^{X_i^TB})^2} = \~~frac{p_i^2}{e^{X_i^TB}}~~ \frac{\partial X_i^TB}{\partial x}</math>	+	<math>\frac{\partial p_i}{\partial x} = \left[ e^{X_i^tB} (1 + e^{X_i^tB})^{-1} + e^{X_i^tB} \left( -e^{X_i^tB} (1 + e^{X_i^tB})^{-2} \right) \right] \frac{\partial X_i^TB}{\partial x}</math>

-	~~where~~ <math>\frac{\partial ~~X_i^TB~~}{\partial x}~~</math> is <math>1</math> for <math>~~\~~mu</math>, <math>g_i</math> for <math>a</math>, and <math>~~\~~mathbf~~{1}_{~~g_i=~~1}</math> ~~for <math>d</math>.~~	+	<math>\frac{\partial p_i}{\partial x} = \left[ \frac{e^{X_i^tB} (1 + e^{X_i^tB}) - (e^{X_i^tB})^2}{(1 + e^{X_i^tB})^2} \right] \frac{\partial X_i^TB}{\partial x}</math>

-	~~Second, we need to calculate the second derivatives of~~ <math>f</math>:	+	<math>\frac{\partial p_i}{\partial x} = \left[ p_i (1 - p_i) \right] \frac{\partial X_i^TB}{\partial x}</math>

-	<math>\frac{\partial^~~2 f~~}{\partial x^2} ~~= - \frac{~~1}{\~~sigma_~~\~~mu^2} +~~ \~~sum_~~{i=1}~~^N \left(-\frac~~{~~y_i}{e^{X_i^TB}} - \frac{(~~1~~-y_i)p_i^2}{e^{X_i^TB}(1-p_i)^2~~} \~~right)~~ \~~frac{~~\~~partial X_i^TB}{\partial x}~~</math>	+	where <math>\frac{\partial X_i^TB}{\partial x}</math> is equal to <math>1, \, g_i, \, \mathbf{1}_{g_i=1}</math> when <math>x</math> is equal respectively to <math>\mu, \, a, \, d</math>.

-	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star</math>. As a consequence, we have the right to use ~~the~~ Laplace's method to approximate the integral around its maximum.	+	This simplifies the first derivatives of <math>f</math> into:
		+
		+	<math>\frac{\partial f}{\partial x} = - \frac{x}{N \, \sigma_x^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial x}</math>
		+
		+	When setting <math>\frac{\partial f}{\partial x}(x^\star) = 0</math>, we observe that <math>x^\star</math> is present not only alone but also inside the sum, in the <math>p_i</math>'s. This means that an iterative procedure is required, typically [http://en.wikipedia.org/wiki/Newton_method_in_optimization Newton's method].
		+
		+	To use it, we need the second derivatives of <math>f</math>:
		+
		+	<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{N \, \sigma_x^2} + \frac{1}{N} \sum_{i=1}^N \left[ (-p_i(1-p_i)\frac{\partial X_i^TB}{\partial x}) + (y_i-p_i)\frac{\partial^2 X_i^TB}{\partial x^2} \right]</math>
		+
		+	The second derivatives of <math>X_i^TB</math> are all equal to 0:
		+
		+	<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{N \, \sigma_x^2} - \frac{1}{N} \sum_{i=1}^N p_i(1-p_i)\frac{\partial X_i^TB}{\partial x}</math>
		+
		+	Note that the second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star</math>. As a consequence, we have the right to use Laplace's method to approximate the integral around its maximum.

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

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add details on logistic regression

2013-01-27T16:47:28Z

Bayesian model of univariate linear regression for QTL detection: add details on logistic regression

←Older revision		Revision as of 16:47, 27 January 2013
Line 290:		Line 290:
	There are many equivalent ways to write the likelihood, the usual one being:		There are many equivalent ways to write the likelihood, the usual one being:

-	<math>y_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i) ~~\; \text{~~ with the log-odds being ~~} \;~~ \mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1}</math>	+	<math>y_i \| p_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i)</math> with the [http://en.wikipedia.org/wiki/Log-odds log-odds] (logit function) being <math>\mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1}</math>

-	~~Using~~ <math>X_i</math> to denote the <math>i</math>-th row of the design matrix <math>X</math> ~~and keeping~~ the same definition as above for <math>B</math>, we have:	+	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>		<math>p_i = \frac{e^{X_i^TB}}{1 + e^{X_i^TB}}</math>
Line 298:		Line 298:
	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:		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) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i}</math>	+	<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) ~~and the Bayes factor now is~~:	+	We still use the same prior as above for <math>B</math> (but there is no <math>\tau</math> anymore), so that:

-	<math>\mathrm{BF} = \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>	+	<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 calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).		The interesting point here is that there is no way to calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).
Line 308:		Line 314:
	Starting with the numerator:		Starting with the numerator:

-	<math>\mathsf{P}(Y\|X,B) = \int \exp \left( \mathrm{ln} \, \mathsf{P}(B) + \mathrm{ln} \, \mathsf{P}(Y \| X, B) \right) \mathsf{d}B</math>	+	<math>\mathsf{P}(Y\|X,M1) = \int \exp \left( \mathrm{ln} \, \mathsf{P}(B) + \mathrm{ln} \, \mathsf{P}(Y \| X, B) \right) \mathsf{d}B</math>

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

	Let's use <math>f</math> to denote the function inside the exponential:		Let's use <math>f</math> to denote the function inside the exponential:

-	<math>\mathsf{P}(Y\|X,B) = \int \exp \left( f(B) \right) \mathsf{d}B</math>	+	<math>\mathsf{P}(Y\|X,M1) = \int \exp \left( 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{3}{2} \mathrm{ln}(2 \pi) - \frac{1}{2} \mathrm{ln}(\|\Sigma_B\|) - \frac{1}{2}(B^T \Sigma_B^{-1} B) + \sum_{i=1}^N \left( y_i \, X_i^T B - \mathrm{ln}(1 + e^{X_i^TB}) \right)</math>
		+
		+	This function will then be used to approximate the integral, like this:
		+
		+	<math>\mathsf{P}(Y\|X) \approx (2 \pi)^{-3/2} \|H(B^\star)\|^{-1/2} \exp^{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.

	First we need to calculate the first derivatives of <math>f</math>, noting that they all have a very similar form:		First we need to calculate the first derivatives of <math>f</math>, noting that they all have a very similar form:
Line 332:		Line 350:
	<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{\sigma_\mu^2} + \sum_{i=1}^N \left(-\frac{y_i}{e^{X_i^TB}} - \frac{(1-y_i)p_i^2}{e^{X_i^TB}(1-p_i)^2} \right) \frac{\partial X_i^TB}{\partial x}</math>		<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{\sigma_\mu^2} + \sum_{i=1}^N \left(-\frac{y_i}{e^{X_i^TB}} - \frac{(1-y_i)p_i^2}{e^{X_i^TB}(1-p_i)^2} \right) \frac{\partial X_i^TB}{\partial x}</math>

-	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star ~~= [\mu^\star a^\star d^\star]^T~~</math>. As a consequence, we have the right to use the Laplace's method to approximate the integral around its maximum.	+	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star</math>. As a consequence, we have the right to use the Laplace's method to approximate the integral around its maximum.

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

	implementation: in R -> to do		implementation: in R -> to do

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add 2nd deriv for binary phenotypes

2013-01-20T21:00:29Z

Bayesian model of univariate linear regression for QTL detection: add 2nd deriv for binary phenotypes

←Older revision		Revision as of 21:00, 20 January 2013
Line 327:		Line 327:

	where <math>\frac{\partial X_i^TB}{\partial x}</math> is <math>1</math> for <math>\mu</math>, <math>g_i</math> for <math>a</math>, and <math>\mathbf{1}_{g_i=1}</math> for <math>d</math>.		where <math>\frac{\partial X_i^TB}{\partial x}</math> is <math>1</math> for <math>\mu</math>, <math>g_i</math> for <math>a</math>, and <math>\mathbf{1}_{g_i=1}</math> for <math>d</math>.
-
-	~~This gives us one equation per parameter:~~
-
-	~~<math>\frac{\partial f}{\partial \mu} = - \frac{\mu}{\sigma_\mu^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right)</math>~~
-
-	~~<math>\frac{\partial f}{\partial a} = - \frac{a}{\sigma_a^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) g_i</math>~~
-
-	~~<math>\frac{\partial f}{\partial d} = - \frac{d}{\sigma_d^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \mathbf{1}_{g_i=1}</math>~~

	Second, we need to calculate the second derivatives of <math>f</math>:		Second, we need to calculate the second derivatives of <math>f</math>:

-	~~...~~	+	<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{\sigma_\mu^2} + \sum_{i=1}^N \left(-\frac{y_i}{e^{X_i^TB}} - \frac{(1-y_i)p_i^2}{e^{X_i^TB}(1-p_i)^2} \right) \frac{\partial X_i^TB}{\partial x}</math>

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

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ minor changes

2013-01-11T22:04:40Z

Bayesian model of univariate linear regression for QTL detection: minor changes

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add info for binary phenotypes

2013-01-03T23:01:40Z

Bayesian model of univariate linear regression for QTL detection: add info for binary phenotypes

←Older revision		Revision as of 23:01, 3 January 2013
Line 300:		Line 300:
	<math>\mathcal{L}(B) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i}</math>		<math>\mathcal{L}(B) = \prod_{i=1}^N p_i^{y_i} (1-p_i)^{1-y_i}</math>

-	We still use the same ~~priors~~ as above for <math>B</math> and the Bayes factor now is:	+	We still use the same prior as above for <math>B</math> (but there is no <math>\tau</math> anymore) and the Bayes factor now is:

	<math>\mathrm{BF} = \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>		<math>\mathrm{BF} = \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 calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).		The interesting point here is that there is no way to calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).
		+
		+	Let's start with the numerator:
		+
		+	<math>\mathsf{P} (Y\|X) = \int \exp \left( \mathrm{ln} \, \mathsf{P}(B) + \mathrm{ln} \, \mathsf{P}(Y \| X, B) \right) \mathsf{d}B</math>
		+
		+	<math>\mathsf{P} (Y\|X) = \int \exp \left[ \mathrm{ln} \left( (2 \pi)^{-\frac{3}{2}} \, \frac{1}{\sigma_\mu \sigma_a \sigma_d} \, \exp\left( -\frac{1}{2} (\frac{\mu^2}{\sigma_\mu^2} + \frac{a^2}{\sigma_a^2} + \frac{d^2}{\sigma_d^2}) \right) \right) + \sum_{i=1}^N y_i \, \mathrm{ln} (p_i) + \sum_{i=1}^N (1-y_i) \, \mathrm{ln} (1-p_i) \right] \mathsf{d}B</math>
		+
		+	using <math>f</math> to denote the function inside the exponential:
		+
		+	<math>\mathsf{P} (Y\|X) = \int \exp \left( f(B) \right) \mathsf{d}B</math>
		+
		+	First we need to calculate the first derivatives of <math>f</math>, by noting that they all have a very similar form:
		+
		+	<math>\frac{\partial f}{\partial x} = - \frac{x}{\sigma_x^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>
		+
		+	where <math>x</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.
		+
		+	Also, there is a simple form for the first derivatives of <math>p_i</math>:
		+
		+	<math>\frac{\partial p_i}{\partial x} = \frac{ (\frac{\partial X_i^TB}{\partial x} e^{X_i^TB}) (1+e^{X_i^TB}) - e^{X_i^TB} (\frac{\partial X_i^TB}{\partial x} e^{X_i^TB})}{(1+e^{X_i^TB})^2} = \frac{p_i^2}{e^{X_i^TB}} \frac{\partial X_i^TB}{\partial x}</math>
		+
		+	where <math>\frac{\partial X_i^TB}{\partial x}</math> is <math>1</math> for <math>\mu</math>, <math>g_i</math> for <math>a</math>, and <math>\mathbf{1}_{g_i=1}</math> for <math>d</math>.
		+
		+	This gives us one equation per parameter:
		+
		+	<math>\frac{\partial f}{\partial \mu} = - \frac{\mu}{\sigma_\mu^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right)</math>
		+
		+	<math>\frac{\partial f}{\partial a} = - \frac{a}{\sigma_a^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) g_i</math>
		+
		+	<math>\frac{\partial f}{\partial d} = - \frac{d}{\sigma_d^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \mathbf{1}_{g_i=1}</math>
		+
		+	Second, we need to calculate the second derivatives of <math>f</math>:
		+
		+	...
		+
		+	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex and hence has a unique global maximum, at <math>B^\star = [\mu^\star a^\star d^\star]^T</math>. As a consequence, we have the right to use the Laplace's method to approximate each integral of the Bayes factor around their respective maximum.
		+
		+	finding the maximums: iterative procedure or generic solver -> to do
		+
		+	implementation: in R -> to do
		+
		+	finding the effect sizes and their std error: to do

Timothee Flutre: /* Bayesian model of univariate linear regression for QTL detection */ add refs

2013-01-03T20:37:26Z

Bayesian model of univariate linear regression for QTL detection: add refs

←Older revision		Revision as of 20:37, 3 January 2013
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 305:		Line 307:


-	* '''Link between Bayes factor and P-value''': see ~~Wakeley~~ (2008)	+	* '''Link between Bayes factor and P-value''': see Wakefield (2008)

	to do		to do
Line 325:		Line 327:


-	* '''Genetic relatedness''': linear mixed model	+	* '''Genetic relatedness''': linear mixed model, 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 as from RNA-seq, Poisson-like likelihood, see Sun (Biometrics, 2012)

	to do		to do

←Older revision		Revision as of 22:39, 3 February 2013
Line 336:		Line 336:
	We therefore need to find <math>B^\star</math>. As it maximizes <math>f</math>, we need to calculate the first derivatives of <math>f</math>. Let's do this the univariate way:		We therefore need to find <math>B^\star</math>. As it maximizes <math>f</math>, we need to calculate the first derivatives of <math>f</math>. Let's do this the univariate way:

-	<math>\frac{\partial f}{\partial x} = - \frac{x}{N \, \~~sigma_x~~^2} + \frac{1}{N} \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>	+	<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial \beta}</math>

-	where <math>x</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.	+	where <math>\beta</math> is <math>\mu</math>, <math>a</math> or <math>d</math>.

	A simple form for the first derivatives of <math>p_i</math> also exists when writing <math>p_i = e^{X_i^tB} (1 + e^{X_i^tB})^{-1}</math>:		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 x} = \left[ e^{X_i^tB} (1 + e^{X_i^tB})^{-1} + e^{X_i^tB} \left( -e^{X_i^tB} (1 + e^{X_i^tB})^{-2} \right) \right] \frac{\partial X_i^TB}{\partial x}</math>	+	<math>\frac{\partial p_i}{\partial \beta} = \left[ e^{X_i^tB} (1 + e^{X_i^tB})^{-1} + e^{X_i^tB} \left( -e^{X_i^tB} (1 + e^{X_i^tB})^{-2} \right) \right] \frac{\partial X_i^TB}{\partial \beta}</math>

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

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

-	where <math>\frac{\partial X_i^TB}{\partial x}</math> is equal to <math>1, \, g_i, \, \mathbf{1}_{g_i=1}</math> when <math>x</math> ~~is equal~~ respectively to <math>\mu, \, a, \, d</math>.	+	where <math>\frac{\partial X_i^TB}{\partial \beta}</math> is equal to <math>1, \, g_i, \, \mathbf{1}_{g_i=1}</math> when <math>\beta</math> corresponds respectively to <math>\mu, \, a, \, d</math>.

	This simplifies the first derivatives of <math>f</math> into:		This simplifies the first derivatives of <math>f</math> into:

-	<math>\frac{\partial f}{\partial x} = - \frac{x}{N \, \~~sigma_x~~^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial x}</math>	+	<math>\frac{\partial f}{\partial \beta} = - \frac{\beta}{N \, \sigma_\beta^2} + \frac{1}{N} \sum_{i=1}^N (y_i - p_i ) \frac{\partial X_i^TB}{\partial \beta}</math>

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

	To use it, we need the second derivatives of <math>f</math>:		To use it, we need the second derivatives of <math>f</math>:

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

	The second derivatives of <math>X_i^TB</math> are all equal to 0:		The second derivatives of <math>X_i^TB</math> are all equal to 0:

-	<math>\frac{\partial^2 f}{\partial x^2} = - \frac{1}{N \, \~~sigma_x~~^2} - \frac{1}{N} \sum_{i=1}^N p_i(1-p_i)\frac{\partial X_i^TB}{\partial x}</math>	+	<math>\frac{\partial^2 f}{\partial \beta^2} = - \frac{1}{N \, \sigma_\beta^2} - \frac{1}{N} \sum_{i=1}^N p_i(1-p_i)\frac{\partial X_i^TB}{\partial \beta}</math>

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

←Older revision		Revision as of 22:04, 11 January 2013
Line 290:		Line 290:
	There are many equivalent ways to write the likelihood, the usual one being:		There are many equivalent ways to write the likelihood, the usual one being:

-	<math>y_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i) \; \text{ with } \; \mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1}</math>	+	<math>y_i \; \overset{i.i.d}{\sim} \; Bernoulli(p_i) \; \text{ with the log-odds being } \; \mathrm{ln} \frac{p_i}{1 - p_i} = \mu + a \, g_i + d \, \mathbf{1}_{g_i=1}</math>

	Using <math>X_i</math> to denote the <math>i</math>-th row of the design matrix <math>X</math> and keeping the same definition as above for <math>B</math>, we have:		Using <math>X_i</math> to denote the <math>i</math>-th row of the design matrix <math>X</math> and keeping the same definition as above for <math>B</math>, we have:
Line 306:		Line 306:
	The interesting point here is that there is no way to calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).		The interesting point here is that there is no way to calculate these integrals analytically. Therefore, we will use [http://en.wikipedia.org/wiki/Laplace_approximation Laplace's method] to approximate them, as in Guan & Stephens (2008).

-	~~Let's start~~ with the numerator:	+	Starting with the numerator:

-	<math>\mathsf{P} (Y\|X) = \int \exp \left( \mathrm{ln} \, \mathsf{P}(B) + \mathrm{ln} \, \mathsf{P}(Y \| X, B) \right) \mathsf{d}B</math>	+	<math>\mathsf{P}(Y\|X,B) = \int \exp \left( \mathrm{ln} \, \mathsf{P}(B) + \mathrm{ln} \, \mathsf{P}(Y \| X, B) \right) \mathsf{d}B</math>

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

-	~~using~~ <math>f</math> to denote the function inside the exponential:	+	Let's use <math>f</math> to denote the function inside the exponential:

-	<math>\mathsf{P} (Y\|X) = \int \exp \left( f(B) \right) \mathsf{d}B</math>	+	<math>\mathsf{P}(Y\|X,B) = \int \exp \left( f(B) \right) \mathsf{d}B</math>

-	First we need to calculate the first derivatives of <math>f</math>, by noting that they all have a very similar form:	+	First we need to calculate the first derivatives of <math>f</math>, noting that they all have a very similar form:

	<math>\frac{\partial f}{\partial x} = - \frac{x}{\sigma_x^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>		<math>\frac{\partial f}{\partial x} = - \frac{x}{\sigma_x^2} + \sum_{i=1}^N \left(\frac{y_i}{p_i} - \frac{1-y_i}{1-p_i} \right) \frac{\partial p_i}{\partial x}</math>
Line 340:		Line 340:
	...		...

-	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex ~~and hence~~ has a unique global maximum, at <math>B^\star = [\mu^\star a^\star d^\star]^T</math>. As a consequence, we have the right to use the Laplace's method to approximate ~~each integral of~~ the ~~Bayes factor~~ around ~~their respective~~ maximum.	+	The second derivatives of <math>f</math> are strictly negative. Therefore, <math>f</math> is globally convex, which means that it has a unique global maximum, at <math>B^\star = [\mu^\star a^\star d^\star]^T</math>. As a consequence, we have the right to use the Laplace's method to approximate the integral around its maximum.

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