Bayesian Regression Models for Whole-Genome Analyses

Meuwissen et al. (2001) introduced three regression models for whole-genome prediction of breeding value of the form

$$y_{i}=\mu+\sum_{j=1}^{k}X_{ij}\alpha_{j}+e_{i},$$

where $y_{i}$ is the phenotypic value, $\mu$ is the intercept, $X_{ij}$ is $j^{th}$ marker covariate of animal $i$, $\alpha_{j}$ is the partial regression coefficient of $X_{ij}$, and $e_{i}$ are identically and independently distributed residuals with mean zero and variance $\sigma_{e}^{2}.$ In most current analyses, $X_{ij}$ are SNP genotype covariates that can be coded as 0, 1 and 2, depending on the number of B alleles at SNP locus $j$.

In all three of their models, a flat prior was used for the intercept and a scaled inverted chi-square distribution for $\sigma_{e}^{2}.$ The three models introduced by Meuwissen et al. (2001) differ only in the prior used for $\alpha_{j}$.

BLUP

In their first model, which they called BLUP, a normal distribution with mean zero and known variance, $\sigma_{\alpha}^{2}$, is used as the prior for $\alpha_{j}$.

The meaning of $\sigma_{\alpha}^{2}$

Assume the QTL are in the marker panel. Then, the genotypic value $g_{i}$ for a randomly sampled animal $i$ can be written as $$g_{i}=\mu+\mathbf{x}'_{i}\boldsymbol{\alpha},$$ where $\mathbf{x}'_{i}$ is the vector of SNP genotype covariates and $\boldsymbol{\alpha}$ is the vector of regression coefficients. Note that randomly sampled animals differ only in $\mathbf{x}'_{i}$ and have $\boldsymbol{\alpha}$ in common. Thus, genotypic variability is entirely due to variability in the genotypes of animals. So, $\sigma_{\alpha}^{2}$ is not the genetic variance at a locus (Fernando et al., 2007; Gianola et al., 2009).

Relationship of $\sigma_{\alpha}^{2}$ to genetic variance

Assume loci with effect on trait are in linkage equilibrium. Then, the additive genetic variance is $$V_{A}=\sum_{j}^{k}2p_{j}q_{j}\alpha_{j}^{2},$$ where $p_{j}=1-q_{j}$ is gene frequency at SNP locus $j$. Letting $U_{j}=2p_{j}q_{j}$ and $V_{j}=\alpha_{j}^{2}$, $$V_{A}=\sum_{j}^{k}U_{j}V_{j}.$$ For a randomly

sampled locus, covariance between $U_{j}$ and $V_{j}$ is $$C_{UV}=\frac{\sum_{j}U_{j}V_{j}}{k}-(\frac{\sum_{j}U_{j}}{k})(\frac{\sum_{j}V_{j}}{k})$$ Rearranging this expression for $C_{UV}$ gives $$\sum_{j}U_{j}V_{j}=kC_{UV}+(\sum_{j}U_{j})(\frac{\sum_{j}V_{j}}{k})$$ So,

$$V_{A}=kC_{UV}+(\sum_{j}2p_{j}q_{j})(\frac{\sum_{j}\alpha_{j}^{2}}{k}).$$ Letting $\sigma_{\alpha}^{2}=\frac{\sum_{j}\alpha_{j}^{2}}{k}$ gives $$V_{A}=kC_{UV}+(\sum_{j}2p_{j}q_{j})\sigma_{\alpha}^{2}$$ and

$$\sigma_{\alpha}^{2}=\frac{V_{A}-kC_{UV}}{\sum_{j}2p_{j}q_{j}},$$ which gives $$\sigma_{\alpha}^{2}=\frac{V_{A}}{\sum_{j}2p_{j}q_{j}},$$ if gene frequency is independent of the effect of the gene.

Full-conditionals:

The joint posterior for all the parameters is proportional to

$$\begin{aligned} f(\boldsymbol{\theta}|\mathbf{y}) & \propto f(\mathbf{y}|\boldsymbol{\theta})f(\boldsymbol{\theta})\\ & \propto \left(\sigma_{e}^{2}\right)^{-n/2}\exp\left\{ -\frac{(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})'(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})}{2\sigma_{e}^{2}}\right\} \\ & \times \prod_{j=1}^{k}\left(\sigma_{\alpha}^{2}\right)^{-1/2}\exp\left\{ -\frac{\alpha_{j}^{2}}{2\sigma_{\alpha}^{2}}\right\} \\ & \times (\sigma_{\alpha}^{2})^{-(\nu_{\alpha}+2)/2}\exp\left\{ -\frac{\nu_{\alpha}S_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}\right\} \\ & \times (\sigma_{e}^{2})^{-(2+\nu_{e})/2}\exp\left\{ -\frac{\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} ,\end{aligned}$$

where $\boldsymbol{\theta}$ denotes all the unknowns.

Full-conditional for $\mu$

The full-conditional for $\mu$ is a normal distribution with mean $\hat{\mu}$ and variance $\frac{\sigma_{e}^{2}}{n}$, where $\hat{\mu}$ is the least-squares estimate of $\mu$ in the model $$\mathbf{y-\sum_{j=1}^{k}}\mathbf{X}_{j}\alpha_{j}=\mathbf{1}\mu+\mathbf{e},$$ and $\frac{\sigma_{e}^{2}}{n}$ is the variance of this estimator ($n$ is the number of observations).

Full-conditional for $\alpha_{j}$

$$\begin{aligned} f(\alpha_{j}|\text{ELSE}) & \propto \exp\left\{ -\frac{(\mathbf{w}_{j}-\mathbf{X}_{j}\alpha_{j})'(\mathbf{w}_{j}-\mathbf{X}_{j}\alpha_{j})}{2\sigma_{e}^{2}}\right\} \\ & \times \exp\left\{ -\frac{\alpha_{j}^{2}}{2\sigma_{\alpha}^{2}}\right\} \\ & \propto \exp\left\{ -\frac{[\mathbf{w}'_{j}\mathbf{w}_{j}-2\mathbf{w}'_{j}\mathbf{X}_{j}\alpha_{j}+\alpha_{j}^{2}(\mathbf{x}'_{j}\mathbf{x}_{j}+\sigma_{e}^{2}/\sigma_{\alpha}^{2})]}{2\sigma_{e}^{2}}\right\} \\ & \propto \exp\left\{ -\frac{(\alpha_{j}-\hat{\alpha_{j}})^{2}}{\frac{2\sigma_{e}^{2}}{(\mathbf{x}'_{j}\mathbf{x}_{j}+\sigma_{e}^{2}/\sigma_{\alpha}^{2})}}\right\} ,\end{aligned}$$

where $$\mathbf{w}_{j}=\mathbf{y}-\mathbf{1}\mu-\sum_{l\neq j}\mathbf{X}_{l}\alpha_{l}.$$ So, the full-conditional for $\alpha_{j}$ is a normal distribution with mean $$\hat{\alpha}_{j}=\frac{\mathbf{X}'_{j}\mathbf{w}_{j}}{(\mathbf{x}'_{j}\mathbf{x}_{j}+\sigma_{e}^{2}/\sigma_{\alpha}^{2})}$$ and variance $\frac{\sigma_{e}^{2}}{(\mathbf{x}'_{j}\mathbf{x}_{j}+\sigma_{e}^{2}/\sigma_{\alpha}^{2})}$.

Full-conditional for $\sigma_{\alpha}^{2}$

$$\begin{aligned} f(\sigma_{\alpha}^{2}|\text{ELSE}) & \propto \prod_{j=1}^{k}\left(\sigma_{\alpha}^{2}\right)^{-1/2}\exp\left\{ -\frac{\alpha_{j}^{2}}{2\sigma_{\alpha}^{2}}\right\} \\ & \times (\sigma_{\alpha}^{2})^{-(\nu_{\alpha}+2)/2}\exp\left\{ -\frac{\nu_{\alpha}S_{\alpha}^{2}}{2\sigma_{\alpha}^{2}}\right\} \\ & \propto (\sigma_{\alpha}^{2})^{-(k+\nu_{\alpha}+2)/2}\exp\left\{ -\frac{\sum_{j=1}^{k}\alpha_{j}^{2}+\nu_{\alpha}S_{\beta\alpha}^{2}}{2\sigma_{\alpha}^{2}}\right\} ,\end{aligned}$$

and this is proportional to a scaled inverted chi-square distribution with $\tilde{\nu}_{\alpha}=\nu_{\alpha}+k$ and scale parameter $\tilde{S}_{\alpha}^{2}=(\sum_{k}\alpha_{j}^{2}+\nu_{\alpha}S_{\alpha}^{2})/\tilde{\nu}_{\alpha}$.

Full-conditional for $\sigma_{e}^{2}$

$$\begin{aligned} f(\sigma_{e}^{2}|\text{ELSE}) & \propto \left(\sigma_{e}^{2}\right)^{-n/2} \exp\left\{ -\frac{ (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})' (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j}) }{2\sigma_{e}^{2}}\right\} \\ & \times (\sigma_{e}^{2})^{-(2+\nu_{e})/2}\exp\left\{ -\frac{\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} \\ & \propto (\sigma_{e}^{2})^{-(n+2+\nu_{e})/2}\exp\left\{ -\frac{(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})'(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})+\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} ,\end{aligned}$$

which is proportional to a scaled inverted chi-square density with $\tilde{\nu}_{e}=n+\nu_{e}$ degrees of freedom and $\tilde{S_{e}^{2}}=\frac{(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})'(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\alpha_{j})+\nu_{e}S_{e}^{2}}{\tilde{\nu_{e}}}$ scale parameter.

BayesB

Model

The usual model for BayesB is:

$$y_{i}=\mu+\sum_{j=1}^{k}X_{ij}\alpha_{j}+e_{i},$$ where the prior $\mu$ is flat and the prior for $\alpha_{j}$ is a mixture distribution:

$$\alpha_{j}=\begin{cases} 0 & \text{probability}\,\pi\\ \sim N(0,\sigma_{j}^{2}) & \text{probability}\,(1-\pi) \end{cases}, $$

where $\sigma_{j}^{2}$ has a scaled inverted chi-square prior with scale parameter $S_{\alpha}^{2}$ and $\nu_{\alpha}$ degrees of freedom. The residual is normally distributed with mean zero and variance $\sigma_{e}^{2}$, which has a scaled inverted chi-square prior with scale parameter $S_{e}^{2}$ and $\nu_{e}$ degrees of freedom. Meuwissen et al. @Meuwissen.THE.ea.2001a gave a Metropolis-Hastings sampler to jointly sample $\sigma_{j}^{2}$ and $\alpha_{j}$. Here, we will show how the Gibbs sampler can be used in BayesB.

In order to use the Gibbs sampler, the model is written as

$$y_{i}=\mu+\sum_{j=1}^{k}X_{ij}\beta_{j}\delta_{j}+e_{i},$$ where $\beta_{j}\sim N(0,\sigma_{j}^{2})$ and $\delta_{j}$ is Bernoulli($1-\pi$):

$$\begin{aligned} \delta_{j} & = & \begin{cases} 0 & \text{probability}\,\pi\\ 1 & \text{probability}\,(1-\pi) \end{cases}.\\\end{aligned}$$

Other priors are the same as in the usual model. Note that in this model, $\alpha_{j}=\beta_{j}\delta_{j}$ has a mixture distribution as in the usual BayesB model.

Full-conditionals:

The joint posterior for all the parameters is proportional to

$$\begin{aligned} f(\boldsymbol{\theta}|\mathbf{y}) & \propto f(\mathbf{y}|\boldsymbol{\theta})f(\boldsymbol{\theta})\\ & \propto \left(\sigma_{e}^{2}\right)^{-n/2}\exp\left\{ -\frac{ (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})' (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j}) }{2\sigma_{e}^{2}}\right\} \\ & \times \prod_{j=1}^{k}\left(\sigma_{j}^{2}\right)^{-1/2}\exp\left\{ -\frac{\beta_{j}^{2}}{2\sigma_{j}^{2}}\right\} \\ & \times \prod_{j=1}^{k}\pi^{(1-\delta_{j})}(1-\pi)^{\delta_{j}}\\ & \times \prod_{j=1}^{k}(\sigma_{j}^{2})^{-(\nu_{\beta}+2)/2}\exp\left\{ -\frac{\nu_{\beta}S_{\beta}^{2}}{2\sigma_{j}^{2}}\right\} \\ & \times (\sigma_{e}^{2})^{-(2+\nu_{e})/2}\exp\left\{ -\frac{\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} ,\end{aligned}$$

where $\boldsymbol{\theta}$ denotes all the unknowns.

Full-conditional for $\mu$

The full-conditional for $\mu$ is a normal distribution with mean $\hat{\mu}$ and variance $\frac{\sigma_{e}^{2}}{n}$, where $\hat{\mu}$ is the least-squares estimate of $\mu$ in the model $$\mathbf{y-\sum_{j=1}^{k}}\mathbf{X}_{j}\beta_{j}\delta_{j}=\mathbf{1}\mu+\mathbf{e},$$ and $\frac{\sigma_{e}^{2}}{n}$ is the variance of this estimator ($n$ is the number of observations).

Full-conditional for $\beta_{j}$

$$\begin{aligned} f(\beta_{j}|\text{ELSE}) & \propto \exp\left\{ -\frac{(\mathbf{w}_{j}-\mathbf{X}_{j}\beta_{j}\delta_{j})'(\mathbf{w}_{j}-\mathbf{X}_{j}\beta_{j}\delta_{j})}{2\sigma_{e}^{2}}\right\} \\ & \times \exp\left\{ -\frac{\beta_{j}^{2}}{2\sigma_{j}^{2}}\right\} \\ & \propto \exp\left\{ -\frac{[\mathbf{w}'_{j}\mathbf{w}_{j}-2\mathbf{w}'_{j}\mathbf{X}_{j}\beta_{j}\delta_{j}+\beta_{j}^{2}(\mathbf{x}'_{j}\mathbf{x}_{j}\delta_{j}+\sigma_{e}^{2}/\sigma_{j}^{2})]}{2\sigma_{e}^{2}}\right\} \\ & \propto \exp\left\{ -\frac{(\beta_{j}-\hat{\beta_{j}})^{2}}{\frac{2\sigma_{e}^{2}}{(\mathbf{x}'_{j}\mathbf{x}_{j}\delta_{j}+\sigma_{e}^{2}/\sigma_{j}^{2})}}\right\} ,\end{aligned}$$

where $$\mathbf{w}_{j}=\mathbf{y}-\mathbf{1}\mu-\sum_{l\neq j}\mathbf{X}_{l}\beta_{l}\delta_{l}.$$ So, the full-conditional for $\beta_{j}$ is a normal distribution with mean $$\hat{\beta}_{j}=\frac{\mathbf{X}'_{j}\mathbf{w}_{j}\delta_{j}}{(\mathbf{x}'_{j}\mathbf{x}_{j}\delta_{j}+\sigma_{e}^{2}/\sigma_{j}^{2})}$$ and variance $\frac{\sigma_{e}^{2}}{(\mathbf{x}'_{j}\mathbf{x}_{j}\delta_{j}+\sigma_{e}^{2}/\sigma_{j}^{2})}$.

Full-conditional for $\delta_{j}$

$$ \Pr(\delta_{j}=1|\text{ELSE})\propto\frac{h(\delta_{j}=1)}{h(\delta_{j}=1)+h(\delta_{j}=0)}, $$

where $$ h(\delta_{j})=\pi^{(1-\delta_{j})}(1-\pi)^{\delta_{j}}\exp\left\{ -\frac{ (\mathbf{w}_{j}-\mathbf{X}_{j}\beta_{j}\delta_{j})' (\mathbf{w}_{j}-\mathbf{X}_{j}\beta_{j}\delta_{j}) }{2\sigma_{e}^{2}}\right\} . $$

Full-conditional for $\sigma_{j}^{2}$

$$\begin{aligned} f(\sigma_{j}^{2}|\text{ELSE}) &\propto \left(\sigma_{j}^{2}\right)^{-1/2}\exp\left\{ -\frac{\beta_{j}^{2}}{2\sigma_{j}^{2}}\right\} \\ &\times (\sigma_{j}^{2})^{-(\nu_{\beta}+2)/2}\exp\left\{ -\frac{\nu_{\beta}S_{\beta}^{2}}{2\sigma_{j}^{2}}\right\} \\ &\propto (\sigma_{j}^{2})^{-(1+\nu_{\beta}+2)/2}\exp\left\{ -\frac{\beta_{j}^{2}+\nu_{\beta}S_{\beta}^{2}}{2\sigma_{j}^{2}}\right\} ,\end{aligned}$$

and this is proportional to a scaled inverted chi-square distribution with $\tilde{\nu}_{j}=\nu_{\beta}+1$ and scale parameter $\tilde{S}_{j}^{2}=(\beta_{j}^{2}+\nu_{\beta}S_{\beta}^{2})/\tilde{\nu}_{j}$.

Full-conditional for $\sigma_{e}^{2}$

$$\begin{aligned} f(\sigma_{e}^{2}|\text{ELSE}) & \propto \left(\sigma_{e}^{2}\right)^{-n/2}\exp\left\{ -\frac{ (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})' (\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j}) }{2\sigma_{e}^{2}}\right\} \\ & \times (\sigma_{e}^{2})^{-(2+\nu_{e})/2}\exp\left\{ -\frac{\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} \\ & \propto (\sigma_{e}^{2})^{-(n+2+\nu_{e})/2}\exp\left\{ -\frac{(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})'(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})+\nu_{e}S_{e}^{2}}{2\sigma_{e}^{2}}\right\} ,\end{aligned}$$

which is proportional to a scaled inverted chi-square density with $\tilde{\nu}_{e}=n+\nu_{e}$ degrees of freedom and $\tilde{S_{e}^{2}}=\frac{(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})'(\mathbf{y}-\mathbf{1}\mu-\sum\mathbf{X}_{j}\beta_{j}\delta_{j})+\nu_{e}S_{e}^{2}}{\tilde{\nu_{e}}}$ scale parameter.

References

  1. R. L. Fernando, D. Habier, C. Stricker, J. C. M. Dekkers, and L. R. Totir. Genomic selection. Acta Agriculturae Scandinavica, Section A - Animal Science, 57(4):192–195, 2007.

  2. D Gianola, G de los Campos, W G Hill, E Manfredi, and R Fernando. Ad- ditive genetic variability and the bayesian alphabet. Genetics, 183(1):347– 363, Sep 2009.

  3. T. H. E. Meuwissen, B. J. Hayes, and M. E. Goddard. Prediction of total genetic value using genome-wide dense marker maps. Genetics, 157:1819– 1829, 2001.