Hierarchical 1-way Analysis of Variance
The bang package simulates from the posterior distributions involved in certain Bayesian models. See the vignette Introducing bang: Bayesian Analysis, No Gibbs for an introduction. In this vignette we consider the hierarchical 1-way Analysis of variance (ANOVA) model: where αi ∼ N(0, σα2), so that θi = μ + αi ∼ N(μ, σα2), and ϵij ∼ N(0, σ2) and all random variables are independent. Variability of the response variable Y about an overall level μ is decomposed into contributions αi, i = 1, …, I from an explanatory factor indicating membership of group i and a random error term ϵij.
This model has I + 3
parameters: ϕ = (μ, σα, σ)
and α = (α1, …, αI).
Equivalently, we could replace α by θ = (θ1, …, θI).
The full posterior π(α, ϕ ∣ y) = π(α ∣ ϕ, y) π(ϕ ∣ y)
can be factorized into products of lower-dimensional densities. See the
Appendix for details. If it is assumed that
μ and (σα, σ)
are independent a priori and a normal prior is set for μ then the only challenging part of
sampling from π(α, ϕ ∣ y)
is to simulate from the two-dimensional density π(σα, σ ∣ y).
Otherwise, we need to simulate from π(μ, σα, σ ∣ y).
The bang function hanova1 uses the rust
package (Northrop 2017) to simulate from
these densities. In sampling from the marginal posterior density π(ϕ ∣ y)
we use a parameterization designed to improve the efficiency of
sampling. In this case we work with (log σα, log σ).
We illustrate the use of hanova1 using two example
datasets. Unless stated otherwise we use the default hyperprior π(μ, σα, σ) ∝ 1/σ,
that is, a uniform prior for (μ, σα, log σ)
for σα > 0, σ > 0
(see Sections 5.7 and 11.6 of Gelman et al.
(2014)). 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 μ
are important, but it is also interesting to compare the magnitudes of
σα and
σ. If, for example, if σα is much
greater than σ 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 Northrop and Chandler
(2014) for more information.
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.
# Number of observations
length(RCP26_2[, "index"])
#> [1] 61
# Numbers of runs for each GCM
table(RCP26_2[, "GCM"])
#>
#> ACCESS1-0 ACCESS1-3 bcc-csm1-1 bcc-csm1-1-m BNU-ESM
#> 0 0 1 1 1
#> CanESM2 CCSM4 CESM1-BGC CESM1-CAM5 CMCC-CM
#> 5 6 0 3 0
#> CMCC-CMS CNRM-CM5 CSIRO-Mk3-6-0 EC-EARTH FGOALS-g2
#> 0 1 10 2 1
#> FIO-ESM GFDL-CM3 GFDL-ESM2G GFDL-ESM2M GISS-E2-H
#> 3 1 1 1 1
#> GISS-E2-H-CC GISS-E2-R GISS-E2-R-CC HadGEM2-AO HadGEM2-CC
#> 0 1 0 1 0
#> HadGEM2-ES inmcm4 IPSL-CM5A-LR IPSL-CM5A-MR IPSL-CM5B-LR
#> 4 0 4 1 0
#> MIROC-ESM MIROC-ESM-CHEM MIROC5 MPI-ESM-LR MPI-ESM-MR
#> 1 1 3 3 1
#> MRI-CGCM3 NorESM1-M NorESM1-ME
#> 1 1 1
# Number of GCMs with at least one run
sum(table(RCP26_2[, "GCM"]) > 0)
#> [1] 28We use hanova1 to sample from the posterior distribution
of the parameters based on the (default) improper uniform prior for
(μ, log σα, σ),
described in Section 11.6 of Gelman et al.
(2014). This prior fits in to the special case considered in the
Appendix, with an infinite prior variance σ02 for μ.
# 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 (σα, σ) with superimposed density contours. We see that the posterior distribution is located at values for which σα is much greater than σ, by a factor of nearly 10. On the right is a similar plot displayed on the scale used for sampling, (ρ1, ρ2) = (log σα, log σ), 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 μ using a histogram with a superimposed kernel density estimate.
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.
Coagulation time data
In the temperature projection example sampling was conducted on the scale (log σα, log α) and was unproblematic. It would also have been unproblematic had we sampled on the original (σα, σ) 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 Gelman et al. (2014). 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 σα in the data. Unless some prior information about σα is provided the posterior distribution for σα will tend to have a heavy upper tail (Gelman 2006).
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 π(σα, σ ∣ y)
using the rust package’s default settings for the generalized
ratio-of-uniforms algorithm. One solution is to reparameterize to (log σα, log σ),
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.
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 σα is evident in the plot on the left. Parameter transformation to (ρ1, ρ2) = (log σα, log σ) 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.
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)
#> 2.5% 25% 50% 75% 97.5%
#> theta[1] 58.8 60.4 61.2 62.0 63.8
#> theta[2] 64.0 65.2 65.9 66.5 67.9
#> theta[3] 65.7 67.1 67.8 68.4 69.8
#> theta[4] 59.4 60.5 61.1 61.7 62.9
#> mu 54.7 62.2 64.0 65.7 73.2
#> sigma[alpha] 2.0 3.5 5.0 7.9 27.0
#> sigma 1.8 2.2 2.4 2.7 3.4These posterior summaries are similar to those presented in Table 11.3 of Gelman et al. (2014) (where σα is denoted τ), which were obtained using Gibbs sampling.
When the number of groups is small Gelman (2006) advocates the use of a half-Cauchy prior for σα. The code below implements this using independent half-Cauchy priors for σα and σ, 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, σα = 10: in practice σα should be set by considering the problem in hand. See Gelman (2006) for details. We set A to be large enough to result in an effectively flat prior for σ.
Appendix
Consider the hierarchical 1-way ANOVA model where αi ∼ N(0, σα2) and ϵij ∼ N(0, σ2) and all random variables are independent. We specify a prior density π(ϕ) for ϕ = (μ, σα, σ). Let y = {yij, i = 1, …, I, j = 1, …, ni} and α = (α1, …, αI).
Marginal posterior density of ϕ
The joint posterior density of (ϕ, α) satisfies
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 αi in the term in square brackets gives where ai = niσ−2 + σα−2, bi = niσ−2(ȳi⋅ − μ) and $c_i = \sigma^{-2} \sum_{j=1}^{n_i} (y_{ij} - \mu)^2$. Therefore,
The marginal posterior density of ϕ is given by Manipulating ci − bi2/ai gives Therefore, the joint posterior density of ϕ satisfies where $S = \sum_{i=1}^I \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_{i\cdot})^2$, μi = ȳi⋅ and σi2 = σα2 + σ2/ni.
A special case
Suppose that a N(μ0, σ02) prior distribution is specified for μ, so that π(μ) ∝ exp [−σ0−2(μ − μ0)2)/2] and that μ and (σα, σ) are independent a priori. We derive the marginal posterior density for (σα, σ) in this case and the conditional posterior density of μ given (σα, σ). The joint posterior density for (μ, σα, σ) satisfies Completing the square in μ in $\sum_{i=0}^{I} \sigma_i^{-2} (\mu - \mu_i)^2$ gives where $S_j = \sum_{i=0}^{I} \mu_i^{j} \sigma_i^{-2}, \, j = 0, 1, 2$. Therefore,
Retaining only the terms in π(μ, σα, σ ∣ y) that involve μ gives Therefore, μ ∣ σα, σ, y ∼ N(S1S0−1, S0−1).
The conditional posterior density of θ given ϕ
The conditional posterior density of α given ϕ satisfies because these are the only terms in π(ϕ, α ∣ y) that involve αi, i = 1, …, I. Therefore, conditional on ϕ, α1, …, αI are independent a posteriori and αi ∣ ϕ, y ∼ N(biai−1, ai−1), where Noting that θi = μ + αi gives θi ∣ ϕ, y ∼ N(θ̂i, ai−1), where
Factorisations for simulation
In the most general case considered here the factorisation means that we can simulate first from the three-dimensional π(ϕ, ∣y) and then, conditional on the value of ϕ, simulate independently from each of the normal distributions of αi ∣ ϕ, y.
In the special case detailed above the factorisation becomes Therefore, the first stage can be performed by simulating from the two-dimensional π(σα, σ ∣ y) and then, conditional on the value of (σα, σ), from the normal distribution of μ ∣ σα, σ, y.