---
title: "Hierarchical 1-way Analysis of Variance"
author: "Paul J. Northrop"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Hierarchical 1-way Analysis of Variance}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: bang.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(comment = "#>", collapse = TRUE)
set.seed(47)
```

The *bang* package simulates from the posterior distributions involved in certain Bayesian models.  See the vignette [Introducing bang: Bayesian Analysis, No Gibbs](bang-a-vignette.html) for an introduction.  In this vignette we consider the hierarchical 1-way Analysis of variance (ANOVA) model:
\begin{equation*}
Y_{ij} = \mu + \alpha_i + \epsilon_{ij} = \theta_i + \epsilon_{ij}, \quad\mbox{for} \,\, i = 1, \ldots, I, \, j = 1, \ldots, n_i, 
\end{equation*}
where  $\alpha_i \sim N(0, \sigma_\alpha^2)$, so that $\theta_i = \mu + \alpha_i \sim N(\mu, \sigma_\alpha^2)$, and $\epsilon_{ij} \sim N(0, \sigma^2)$ and all random variables are independent.  Variability of the response variable $Y$ about an overall level $\mu$ is decomposed into contributions $\alpha_i, i = 1, \ldots, I$ from an explanatory factor indicating membership of group $i$ and a random error term $\epsilon_{ij}$.

This model has $I + 3$ parameters: $\boldsymbol{\mathbf{\phi}} = (\mu, \sigma_\alpha, \sigma)$ and $\boldsymbol{\mathbf{\alpha}} = (\alpha_1, \ldots, \alpha_I)$.  Equivalently, we could replace $\boldsymbol{\mathbf{\alpha}}$ by $\boldsymbol{\mathbf{\theta}} = (\theta_1, \ldots, \theta_I)$.  The full posterior $\pi(\boldsymbol{\mathbf{\alpha}},\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}}) = \pi(\boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{y}}) \, \pi(\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}})$ can be factorized into products of lower-dimensional densities.  See the [Appendix](#appendix) for details.  If it is assumed that $\mu$ and $(\sigma_\alpha, \sigma)$ are independent *a priori* and a normal prior is set for $\mu$ then the only challenging part of sampling from $\pi(\boldsymbol{\mathbf{\alpha}},\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}})$ is to simulate from the two-dimensional density $\pi(\sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})$.  Otherwise, we need to simulate from $\pi(\mu, \sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})$.  The *bang* function `hanova1` uses the *rust* package [@rust] to simulate from these densities.  In sampling from the marginal posterior density $\pi(\phi \mid \boldsymbol{\mathbf{y}})$ we use a parameterization designed to improve the efficiency of sampling.  In this case we work with $(\log\sigma_\alpha, \log\sigma)$. 
We illustrate the use of `hanova1` using two example datasets.  Unless stated otherwise we use the default hyperprior $\pi(\mu, \sigma_\alpha, \sigma) \propto 1/\sigma$, that is, a uniform prior for $(\mu, \sigma_\alpha, \log\sigma)$ for $\sigma_\alpha>0, \sigma>0$ (see Sections 5.7 and 11.6 of @BDA2014).  A user-defined prior can be set using `set_user_prior`.

## Late 21st Century Global Temperature Projection Data

The data frame `temp2` contains indices of global temperature change from late 20th century (1970-1999) to late 21st century (2069-2098) based on data produced by the Fifth Coupled Model Intercomparison Project (CMIP5).  The dataset is the union of four subsets, each based on a different greenhouse gas emissions scenario called a   Representative Concentration Pathway (RCP). Here we analyse only data for RCP2.6.  Of course, inferences about the overall temperature change parameter $\mu$ are important, but it is also interesting to compare the magnitudes of $\sigma_\alpha$ and $\sigma$.  If, for example, if $\sigma_\alpha$ is much greater than $\sigma$ then uncertainty about temperature projection associated with choice of GCM is greater than that associated with the choice of simulation run from a particular GCM. See @NC2014 for more information.

```{r, include = FALSE}
set.seed(7052009)
```

```{r}
library(bang)
# Extract RCP2.6 data
RCP26_2 <- temp2[temp2$RCP == "rcp26", ]
```

There are 61 observations in total distributed rather unevenly across the 38 GCMs.  Only 28 of the GCMs have at least one run available for RCP2.6.

```{r}
# Number of observations
length(RCP26_2[, "index"])
# Numbers of runs for each GCM
table(RCP26_2[, "GCM"])
# Number of GCMs with at least one run
sum(table(RCP26_2[, "GCM"]) > 0)
```

We use `hanova1` to sample from the posterior distribution of the parameters based on the (default) improper uniform prior for $(\mu, \log\sigma_\alpha, \sigma)$, described in Section 11.6 of @BDA2014.  This prior fits in to the special case considered in the [Appendix](#appendix), with an infinite prior variance $\sigma_0^2$ for $\mu$.

```{r, fig.show='hold', fig.width = 3.45, fig.height = 3.45}
# The response is the index, the explanatory factor is the GCM
temp_res <- hanova1(resp = RCP26_2[, "index"], fac = RCP26_2[, "GCM"])
# Plots relating to the posterior sample of the variance parameters
plot(temp_res, params = "ru")
plot(temp_res, ru_scale = TRUE)
```

The plot on the left shows the values sampled from the posterior distribution of $(\sigma_\alpha, \sigma)$ with superimposed density contours.  We see that the posterior distribution is located at values for which $\sigma_\alpha$ is much greater than $\sigma$, by a factor of nearly 10.  On the right is a similar plot displayed on the scale used for sampling, $(\rho_1, \rho_2)=(\log \sigma_\alpha, \log \sigma)$, after relocation of the posterior mode to the origin and rotation and scaling to near circularity of the density contours.

We summarize the marginal posterior distribution of $\mu$ using a histogram with a superimposed kernel density estimate.

```{r, fig.width = 5, fig.align='center'}
hist(temp_res$sim_vals[, "mu"],  main = "", xlab = expression(mu), prob = TRUE)
lines(density(temp_res$sim_vals[, "mu"]))
```

The following plot summarizes the estimated marginal posterior densities of the mean index for each GCM.
```{r}
plot(temp_res, params = "pop", which_pop = "all", one_plot = TRUE)
```

## Coagulation time data

In the temperature projection example sampling was conducted on the scale $(\log\sigma_\alpha, \log\alpha)$ and was unproblematic.  It would also have been unproblematic had we sampled on the original $(\sigma_\alpha, \sigma)$ scale.  To show that there are examples where the latter is not the case we consider a small dataset presented in Section 11.6 of @BDA2014.  The response variable is the coagulation time of blood drawn from 24 animals allocated to different diets.  The crucial aspect of this dataset is that the explanatory factor (diet) has only 4 levels.  This means that there is little information about $\sigma_\alpha$ in the data.  Unless some prior information about $\sigma_\alpha$ is provided the posterior distribution for $\sigma_\alpha$ will tend to have a heavy upper tail [@Gelman2006]. 

The generalized ratio-of-uniforms method used by *rust* can fail for heavy-tailed densities and this is indeed the case for these data if we try to sample directly from $\pi(\sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})$ using the *rust* package's default settings for the generalized ratio-of-uniforms algorithm.  One solution is to reparameterize to $(\log\sigma_\alpha, \log\sigma)$, which `hanova1` implements by default.  Another possibility is to increase the generalized ratio-of-uniforms tuning parameter `r` from the default value of 1/2 used in *rust*.  These approaches are illustrated below.

```{r, fig.show='hold', fig.width = 3.45, fig.height = 3.45}
coag1 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2], n = 10000)
coag2 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2], n = 10000,
                 param = "original", r = 1)
plot(coag1, params = "ru")
plot(coag1, ru_scale = TRUE)
```

The heaviness of the upper tail of the marginal posterior density of $\sigma_\alpha$ is evident in the plot on the left.  Parameter transformation to $(\rho_1, \rho_2)=(\log\sigma_\alpha, \log\sigma)$ results in a density (in the plot on the right) from which it is easier to simulate.

We produce some summaries of the posterior sample stored in `coag1`.  The summaries calculated from `coag2` are very similar.

```{r}
probs <- c(2.5, 25, 50, 75, 97.5) / 100
all1 <- cbind(coag1$theta_sim_vals, coag1$sim_vals)
round(t(apply(all1, 2, quantile, probs = probs)), 1)
```

These posterior summaries are similar to those presented in Table 11.3 of @BDA2014 (where $\sigma_\alpha$ is denoted $\tau$), which were obtained using Gibbs sampling. 

When the number of groups is small @Gelman2006 advocates the use of a half-Cauchy prior for $\sigma_\alpha$.  The code below implements this using independent half-Cauchy priors for $\sigma_\alpha$ and $\sigma$, that is,
\[ 
\pi(\sigma_\alpha, \sigma) \propto 
\left(1 + \frac{\sigma_\alpha^2}{A_\alpha^2}\right)^{-1}
\left(1 + \frac{\sigma^2}{A^2}\right)^{-1}, \quad \sigma_\alpha>0, \, \sigma>0.
\]
We choose, somewhat arbitrarily, $\sigma_\alpha = 10$: in practice $\sigma_\alpha$ should be set by considering the problem in hand.  See @Gelman2006 for details. We set $A$ to be large enough to result in an effectively flat prior for $\sigma$.

```{r}
coag3 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2],
                 param = "original", prior = "cauchy", hpars = c(10, 1e6))
```

## Appendix {#appendix}

Consider the hierarchical 1-way ANOVA model
\begin{equation*}
Y_{ij} = \mu + \alpha_i + \epsilon_{ij}, \quad\mbox{for} \,\, i = 1, \ldots, I, \, 
j = 1, \ldots, n_i,  \label{eqn:1way}
\end{equation*}
where  $\alpha_i \sim N(0, \sigma_\alpha^2)$ and $\epsilon_{ij} \sim N(0, \sigma^2)$
and all random variables are independent.
We specify a prior density $\pi(\boldsymbol{\mathbf{\phi}})$ for 
$\boldsymbol{\mathbf{\phi}} = (\mu, \sigma_\alpha, \sigma)$.
Let 
$\boldsymbol{\mathbf{y}} = \{y_{ij}, i=1, \ldots, I, j = 1, \ldots, n_i\}$
and 
$\boldsymbol{\mathbf{\alpha}} = (\alpha_1, \ldots, \alpha_I)$.

### Marginal posterior density of $\boldsymbol{\mathbf{\phi}}$

The joint posterior density of 
$(\boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{\alpha}})$ satisfies
\begin{align*}
\pi(\boldsymbol{\mathbf{\alpha}},\boldsymbol{\mathbf{\phi}} 
\mid \boldsymbol{\mathbf{y}}) &\propto 
\pi(\boldsymbol{\mathbf{\phi}})
\prod\limits_{i=1}^{I}\prod\limits_{j=1}^{n_i} 
\sigma^{-1} \exp \left[-\frac{1}{2\sigma^2}(y_{ij} - \mu - \alpha_i)^2\right]
\prod_{i=1}^{I} 
\sigma_\alpha^{-1} \exp\left[-\frac{1}{2\sigma^2_\alpha}\alpha_i^2\right]  \nonumber \\
&= \pi(\boldsymbol{\mathbf{\phi}}) \, \sigma^{-n_\cdot} \, \sigma_\alpha^{-I} 
\prod_{i=1}^I \exp \left\{ -\frac12 \left[ \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu - \alpha_i)^2 
+\sigma_\alpha^{-2} \alpha_i^2 \right] \right\},
\label{eqn:joint}
\end{align*}

where $n_\cdot = \sum_{i=1}^{I} n_i$.
Let $\bar{y}_{i\cdot} = (1/n_i) \sum_{j=1}^{n_i} y_{ij}$.
Completing the\ square in $\alpha_i$ in the term in square brackets gives
\begin{align*}
&\sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu - \alpha_i)^2 + \sigma_\alpha^{-2} \alpha_i^2 \\
&=
\sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu)^2 
- 2 n_i \sigma^{-2} (\bar{y}_{i\cdot} - \mu) \alpha_i + (n_i \sigma^{-2} + \sigma_\alpha^{-2})\alpha_i^2 \\
&= (\sigma^{-2} n_i + \sigma_{\alpha}^{-2}) 
\left[ \alpha_i - \frac{n_i \sigma^{-2} (\bar{y}_{i\cdot} - \mu)}{n_i \sigma^{-2} + \sigma_{\alpha}^{-2}} \right]^2
- \frac{n_i^2 \sigma^{-4} (\bar{y}_{i\cdot} - \mu)^2}{n_i \sigma^{-2} + \sigma_{\alpha}^{-2}}
+ \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu)^2 \\
&= a_i (\alpha_i - b_i / a_i)^2 + c_i - b_i^2 / a_i,
\end{align*}
where
$a_i = n_i \sigma^{-2} + \sigma_\alpha^{-2}, 
b_i = n_i \sigma^{-2} (\bar{y}_{i\cdot} - \mu)$
and
$c_i = \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu)^2$.
Therefore,
\begin{align*}
\pi(\boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{y}}) 
&\propto 
\pi(\boldsymbol{\mathbf{\phi}}) \, \sigma^{-n_\cdot} \, \sigma_\alpha^{-I} 
\prod_{i=1}^I \exp\left[ -\frac12 \left( c_i - b_i ^ 2 / a_i \right) \right]
\exp\left[ -\frac{1}{2a_i^{-1}} \left( \alpha_i - b_i / a_i \right)^2 \right]. \\
\end{align*}

The marginal posterior density of $\boldsymbol{\mathbf{\phi}}$ is given by 
\begin{align*}
\pi(\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}}) 
&=
\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}
\pi(\boldsymbol{\mathbf{\alpha}}, \boldsymbol{\mathbf{\phi}} 
\mid \boldsymbol{\mathbf{y}}) {\rm ~d}\alpha_1 \cdots {\rm ~d}\alpha_I \\
&\propto 
\pi(\boldsymbol{\mathbf{\phi}}) \, \sigma^{-n_\cdot} \, \sigma_\alpha^{-I} 
\prod_{i=1}^I \exp\left[ -\frac12 \left( c_i - b_i ^ 2 / a_i \right) \right]
\int_{-\infty}^{\infty} \exp\left[ -\frac{1}{2a_i^{-1}} \left( \alpha_i - b_i / a_i \right)^2 \right] {\rm ~d}\alpha_i \\
&\propto 
\pi(\boldsymbol{\mathbf{\phi}}) \, \sigma^{-n_\cdot} \, \sigma_\alpha^{-I} 
\prod_{i=1}^I a_i^{-1/2} \exp\left[ -\frac12 \left( c_i - b_i ^ 2 / a_i \right) \right] \\
&\propto 
\pi(\boldsymbol{\mathbf{\phi}}) \, \sigma^{I-n_\cdot} \, 
\prod_{i=1}^I (\sigma_\alpha^2+\sigma^2/n_i)^{-1/2} \exp\left[ -\frac12 \left( c_i - b_i ^2 / a_i \right) \right].
\end{align*}
Manipulating $c_i - b_i ^2 / a_i$ gives
\begin{align*}
c_i - b_i ^ 2 / a_i 
  &= \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu)^2 
   - \frac{n_i^2 \sigma^{-4} (\bar{y}_{i\cdot} - \mu)^2}
  {n_i \sigma^{-2} + \sigma_\alpha^{-2}} \\
  &= \sigma^{-2} \left[ \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2 
   + n_i (\bar{y}_{i\cdot} - \mu)^2 \right]
   - \frac{n_i^2 \sigma^{-4} (\bar{y}_{i\cdot} - \mu)^2}
  {n_i \sigma^{-2} + \sigma_\alpha^{-2}} \\
  &= \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2 +
   \frac{\left[ n_i \sigma^{-2} \left(n_i \sigma^{-2}+\sigma_\alpha^{-2}\right)
    - n_i^2 \sigma^{-4} \right]}{n_i \sigma^{-2} + \sigma_\alpha^{-2}}
    (\bar{y}_{i\cdot} - \mu)^2 \\
  &= \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2 
   + \frac{1}{\sigma_\alpha^2 + \sigma^2 / n_i} (\mu - \bar{y}_{i\cdot})^2.
\end{align*}
Therefore, the joint posterior density of $\boldsymbol{\mathbf{\phi}}$ satisfies
\begin{align*}
\pi(\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}}) 
&\propto 
\pi(\boldsymbol{\mathbf{\phi}}) \, 
\sigma^{I-n_\cdot} \, 
\exp\left( -\frac{1}{2\sigma^{2}} S \right)
\prod_{i=1}^I \sigma_i^{-1} 
\exp\left[ -\frac{1}{2 \sigma_i^2} (\mu - \mu_i)^2 \right],
\end{align*}
where $S = \sum_{i=1}^I \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2$,
$\mu_i = \bar{y}_{i\cdot}$ and $\sigma_i^2 = \sigma_\alpha^2 + \sigma^2/n_i$.

#### A special case

Suppose that a $N(\mu_0, \sigma_0^2$) prior distribution is specified for $\mu$, 
so that $\pi(\mu) \propto \exp[-\sigma_0^{-2}(\mu-\mu_0)^2) / 2]$ and that $\mu$ and
$(\sigma_\alpha, \sigma)$ are independent *a priori*.
We derive the marginal posterior density for $(\sigma_\alpha, \sigma)$ in this case
and the conditional posterior density of $\mu$ given $(\sigma_\alpha, \sigma)$.
The joint posterior density for $(\mu, \sigma_\alpha, \sigma)$ satisfies
\begin{align*}
\pi(\mu, \sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})
&\propto \pi(\mu) \, \pi(\sigma_\alpha, \sigma) \, 
\sigma^{I-n_\cdot} \, 
\exp\left( -\frac{1}{2\sigma^{2}} S \right)
\prod_{i=1}^I \sigma_i^{-1} 
\exp\left[ -\frac{1}{2 \sigma_i^2} (\mu - \mu_i)^2 \right], \\
&=
\pi(\sigma_\alpha, \sigma) \, 
\sigma^{I-n_\cdot} \, 
\exp\left( -\frac{1}{2\sigma^{2}} S \right)
\left(\prod_{i=1}^I \sigma_i^{-1} \right)
\exp\left[
-\frac{1}{2} \sum_{i=0}^{I} \sigma_i^{-2} (\mu - \mu_i)^2 
\right],
\end{align*}
Completing the square in $\mu$ in $\sum_{i=0}^{I} \sigma_i^{-2} (\mu - \mu_i)^2$
gives
\begin{align*}
\sum_{i=0}^{I} \sigma_i^{-2} (\mu - \mu_i)^2 
&= 
\sum_{i=0}^{I} \sigma_i^{-2} (\mu^2 - 2 \mu \mu_i + \mu_i^2) \\
&=
\mu^2 \sum_{i=0}^{I} \sigma_i^{-2} 
- 2 \mu \sum_{i=0}^{I} \mu_i \sigma_i^{-2} 
+ \sum_{i=0}^{I} \mu_i ^2 \sigma_i^{-2} \\ 
&=
\left(\sum_{i=0}^{I} \sigma_i^{-2} \right)
\left( \mu - \frac{\sum_{i=0}^{I} \mu_i \sigma_i^{-2}}
{\sum_{i=0}^{I} \sigma_i^{-2}} \right)^{\!\!2}
- \frac{\left( \sum_{i=0}^{I} \mu_i \sigma_i^{-2} \right)^2}
{\sum_{i=0}^{I} \sigma_i^{-2}}
+ \sum_{i=0}^{I} \mu_i^{2} \sigma_i^{-2} \\
&= 
S_0 (\mu - S_1 S_0^{-1})^2
- S_1^2 S_0^{-1} + S_2, \\
\end{align*}
where $S_j = \sum_{i=0}^{I} \mu_i^{j} \sigma_i^{-2}, \, j = 0, 1, 2$.
Therefore,
\begin{align*}
\pi(\sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})
&=
\int_{-\infty}^{\infty} 
\pi(\mu, \sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})
{\rm ~d}\mu \\
&\propto
\pi(\sigma_\alpha, \sigma) \, 
\sigma^{I-n_\cdot} \, 
\exp\left\{ -\frac{1}{2\sigma^{2}} S \right\}
\left(\prod_{i=1}^I \sigma_i^{-1} \right)
\exp\left( \frac12 S_1^2 S_0^{-1} - \frac12 S_2 \right) \\
&\phantom{\propto} \times 
\int_{-\infty}^{\infty} 
\exp\left\{ -\frac{1}{2 S_0^{-1}} (\mu - S_1 S_0^{-1})^2
\right\}
{\rm ~d}\mu \\
&\propto
\pi(\sigma_\alpha, \sigma) \, 
\sigma^{I-n_\cdot} \, 
\exp\left\{ -\frac{1}{2\sigma^{2}} S \right\}
\left(\prod_{i=1}^I \sigma_i^{-1} \right)
\exp\left( \frac12 S_1^2 S_0^{-1} - \frac12 S_2 \right) 
S_0^{-1/2}.
\end{align*}

Retaining only the terms in $\pi(\mu, \sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})$ that involve $\mu$ gives
\begin{align*}
\pi(\mu \mid \sigma_\alpha, \sigma, \boldsymbol{\mathbf{y}}) 
&\propto 
\exp\left[ -\frac{1}{2 S_0^{-1}} (\mu - S_1 / S_0)^2 \right].
\end{align*}
Therefore, $\mu \mid \sigma_\alpha, \sigma, \boldsymbol{\mathbf{y}}
\sim N(S_1 S_0^{-1}, S_0^{-1})$.

### The conditional posterior density of $\boldsymbol{\mathbf{\theta}}$ given $\boldsymbol{\mathbf{\phi}}$

The conditional posterior density of $\boldsymbol{\mathbf{\alpha}}$ given 
$\boldsymbol{\mathbf{\phi}}$ satisfies 
\begin{align*}
\pi(\boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{\phi}}, 
\boldsymbol{\boldsymbol{\mathbf{y}}}) 
&\propto \prod_{i=1}^I \exp\left[ -\frac{1}{2a_i^{-1}} \left( \alpha_i - b_i / a_i \right)^2 \right],
\end{align*}
because these are the only terms in $\pi (\boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{y}})$ that involve $\alpha_i, i = 1, \ldots, I$.
Therefore, conditional on $\boldsymbol{\mathbf{\phi}}$, $\alpha_1, \ldots, \alpha_I$ are independent *a posteriori*
and $\alpha_i \mid \boldsymbol{\mathbf{\phi}}, 
\boldsymbol{\mathbf{y}} \sim N(b_i a_i^{-1}, a_i^{-1})$,
where
\begin{align*}
b_i a_i^{-1} &= 
\frac{\sigma_\alpha^2 \, (\bar{y}_{i\cdot} - \mu)} 
{\sigma_\alpha^{2} + \sigma^2 / n_i} \\
a_i^{-1} &= \frac{\sigma_\alpha^2 \, \sigma^2 / n_i}
{\sigma_\alpha^{2} + \sigma^2 / n_i}. \\
\end{align*}
Noting that $\theta_i = \mu + \alpha_i$ gives $\theta_i \mid \boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{y}} \sim N(\hat{\theta}_i, a_i^{-1})$, where
\begin{align*}
\hat{\theta}_i &= 
\frac{\bar{y}_{i\cdot} \sigma_\alpha^2 + \mu \, \sigma ^ 2 / n_i} 
{\sigma_\alpha^{2} + \sigma^2 / n_i}. \\
\end{align*}

### Factorisations for simulation

In the most general case considered here the factorisation 
\begin{align*}
\pi(\boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{y}})
&= \pi(\boldsymbol{\mathbf{\phi}} \mid \boldsymbol{\mathbf{y}})
\, \prod_{i=1}^{I} \pi(\alpha_i \mid \boldsymbol{\mathbf{\phi}}, 
\boldsymbol{\boldsymbol{\mathbf{y}}})
\end{align*}
means that we can simulate first from the three-dimensional $\pi(\boldsymbol{\mathbf{\phi}}, \mid \boldsymbol{\mathbf{y}})$ and then, conditional on the value of $\boldsymbol{\mathbf{\phi}}$, simulate independently from each of the normal distributions of $\alpha_i \mid \boldsymbol{\mathbf{\phi}}, \boldsymbol{\boldsymbol{\mathbf{y}}}$.

In the special case detailed above the factorisation becomes
\begin{align*}
\pi(\boldsymbol{\mathbf{\phi}}, \boldsymbol{\mathbf{\alpha}} \mid \boldsymbol{\mathbf{y}})
&= \pi(\sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})
\, \pi(\mu \mid \sigma_\alpha, \sigma, \boldsymbol{\mathbf{y}})
\, \prod_{i=1}^{I} \pi(\alpha_i \mid \boldsymbol{\mathbf{\phi}}, 
\boldsymbol{\boldsymbol{\mathbf{y}}}).
\end{align*}
Therefore, the first stage can be performed by simulating from the two-dimensional $\pi(\sigma_\alpha, \sigma \mid \boldsymbol{\mathbf{y}})$ and then, conditional on the value of $(\sigma_\alpha, \sigma)$, from the normal distribution of $\mu \mid \sigma_\alpha, \sigma, \boldsymbol{\mathbf{y}}$.

## References

<script type="text/x-mathjax-config">
   MathJax.Hub.Config({  "HTML-CSS": { minScaleAdjust: 125, availableFonts: [] }  });
</script>

