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}$.
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}$.
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).
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,
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
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.
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.
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).
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})}$.
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}$.
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.
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.
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.
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).
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})}$.
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\}. $$
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}$.
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.
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.
D Gianola, G de los Campos, W G Hill, E Manfredi, and R Fernando. Additive genetic variability and the bayesian alphabet. Genetics, 183(1):347– 363, Sep 2009.
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.
Here, the scale parameter, $S^2_{\beta}$, in scaled inverse chi-squarred prior for $\sigma^2_j$ is treated as an unknown.
The density function for this prior distribution is: $$ f(\sigma_{j}^{2}|S_{\beta}^{2}, \nu_{\beta}) = \frac{(S_{\beta}^{2}\nu_{\beta}/2)^{\nu_{\beta}/2}} {\Gamma(\nu_{\beta}/2)}(\sigma_{j}^{2})^{-(2+\nu_{\beta})/2}\exp\left\{ -\frac{\nu_{\beta}S_{\beta}^{2}}{2\sigma_{j}^{2}}\right\}, $$
where $S_{\beta}^{2}$ and $\nu_{\beta}$ are the scale and the degrees of freedom parameters for this distribution. The scale parameter $S_{\beta}^{2}$ does not appear anywhere else in $f(\boldsymbol{\theta}|\mathbf{y})$.
Thus, the Gamma$(S^2_{\beta}| a,b)$ distribution with density function: $$ f(S^2_{\beta}|a,b) = \frac {b^a(S^2_{\beta})^{a-1}} {\Gamma(a)} \exp\{-bS^2_{\beta}\} $$ is a conjugate prior for $S^2_{\beta}$.
So, using Gamma$(S^2_{\beta}| a,b)$ as the prior for $S^2_{\beta}$, its full conditional can be written as
$$ \begin{aligned} f(S^2_{\beta}|\text{ELSE}) & \propto \prod_{j=1}^{k}(S^2_{\beta})^{\nu_{\beta}/2}\exp\left\{ -\frac{\nu_{\beta}S_{\beta}^{2}}{2\sigma_{j}^{2}}\right\} \\ & \times b^a(S^2_{\beta})^{a-1}\exp\{-bS^2_{\beta}\}\\ &\times (S^2_{\beta})^{(k\nu_{\beta})/2 + a - 1} \exp\left\{ -S_{\beta}^{2}(\sum_{j=1}^k\frac{\nu_{\beta}}{2\sigma_{j}^{2}} + b)\right\}, \end{aligned} $$which can be recognized as a Gamma distribution with shape $a^*= (k\nu_{\beta})/2 + a$ and rate $b^*=(\sum_{j=1}^k\frac{\nu_{\beta}}{2\sigma_{j}^{2}} + b)$ parameters.