Chapter 5 Hierarchical Model

Consider the following model: \[ \begin{aligned} y_i & \stackrel{ind}{\sim} p(y\mid \theta_i)\\ \theta_i & \stackrel{ind}{\sim} p(\theta\mid \phi) \\ \phi &\sim p(\phi) \end{aligned} \] This is a hierarchical or multilevel model.

The joint posterior distribution can be decomposed via \[ p(\theta, \phi\mid y) = p(\theta\mid \phi, y)p(\phi \mid y) \] where

• \(p(\theta\mid \phi, y) \propto p(y\mid \theta)p(\theta\mid \phi)= \prod_{i = 1}^n p(y_i \mid \theta_i)p(\theta_i \mid \phi) = \prod_{i = 1}^n p(\theta_i \mid \phi, y_i)\)
• $p(y) p(y) p() $
• \(p(y \mid \phi) = \int p(y\mid \theta) p(\theta\mid \phi)d\theta = \prod_{i = 1}^n \int p(y_i\mid \theta_i)p(\theta_i\mid \phi) d\theta_i = \prod_{i = 1}^n p(y_i \mid \phi)\)

Example: Beta-Binomial

Consider \[ Y_i \stackrel{ind}{\sim} Bin(n_i, \theta_i) \quad \theta_i \stackrel{ind}{\sim} Be(\alpha, \beta) \quad \alpha, \beta \sim p(\alpha, \beta) \] In this example, \(\phi = (\alpha, \beta)\). \(p(\theta \mid \alpha, \beta, y) = \prod_{i = 1}^n Be(\theta_i\mid \alpha + y_i, \beta + n_i - y_i)\). The marginal posterior for \(\alpha, \beta\) is \(p(\alpha, \beta \mid y) \propto p(y\mid \alpha, \beta) p(\alpha, \beta)\) with \[ \begin{aligned} p(y\mid \alpha, \beta) &= \prod_{i = 1}^n \int p(y_i\mid \theta_i) p(\theta_i\mid \alpha, \beta) d\theta_i \\ & = \prod_{i = 1}^n \int Bin(y_i\mid n_i, \theta_i)Be(\theta_i \mid \alpha, \beta)d\theta_i \\ & = \prod_{i = 1}^n {n_i\choose y_i}\frac{B(\alpha + y_i , \beta + n_i - y_i)}{B(\alpha, \beta)} \end{aligned} \] Thus, \(y_i \mid \alpha, \beta \stackrel{ind}{\sim} \text{Beta-binomial}(n_i, \alpha, \beta)\).

\(\alpha\): prior successes and \(\beta\): prior failures

A more natural parameterization is

• prior expectation: \(\mu = \frac{\alpha}{\alpha + \beta}\)
• prior sample size \(\eta = \alpha + \beta\)

We can assume the mean (\(\mu\)) and the size (\(\eta\)) are independent a priori. For example, suppose \(\mu \sim Be(20, 30)\) and \(\eta \sim LN(0, 3^2)\) where \(LN\) is short for log-normal distribution.

model_informative_prior = "
data {
int<lower=0> N; // data
int<lower=0> n[N] ;
int<lower=0> y[N];
real<lower=0> a; // prior
real<lower=0> b;
real<lower=0> C;
real m;
}
parameters {
real<lower=0, upper=1> mu;
real<lower=0> eta;
real<lower=0, upper=1> theta[N];
}
transformed parameters {
real<lower=0> alpha;
real<lower=0> beta;
alpha = eta* mu;
beta = eta*(1-mu) ;
}
model {
mu ~ beta(a,b);
eta ~ lognormal(m,C);

# implicit joint distributions
theta ~ beta(alpha,beta);
y ~ binomial(n, theta);
}
"

a = 20; b = 30; m = 0; C = 3
n = 1e4
prior_draws = mutate(data.frame(mu = rbeta(n, a, b), eta = rlnorm(n, m, C)),
                     alpha = eta* mu, beta = eta*(1-mu))

dat = list(y=d$made, n=d$attempts, N=nrow(d), a=a, b=b, m=m, C=C)
m = stan_model(model_code=model_informative_prior)
r = sampling(m, dat, c("mu", "eta", "alpha" "beta", "theta"), iter = 10000)

plot(r, pars = c('eta', 'alpha', 'beta'))
plot(r, pars = c('mu', 'theta'))

In the Bayesian Data Analysis page 110, serveral priors are discussed:

• \((\log(\alpha/\beta), \log(\alpha + \beta)) \propto 1\) leads to an improper posterior
• \((\log(\alpha/\beta), \log(\alpha + \beta)) \sim Unif([-10^{10}, 10^{10}]\times [-10^{10}, 10^{10}])\) while proper and seemingly vague is a very informative prior
• \((\log(\alpha, \beta), \log(\alpha+\beta))\propto \alpha\beta(\alpha+\beta)^{-5/2}\) while leads to a proper posterior and is equivalent to \(p(\alpha, \beta) \propto (\alpha + \beta)^{-5/2}\)

# Stan code for default prior 
model_default_prior = "
data {
int<lower=0> N;
int<lower=0> n [N];
int<lower=0> y [N];
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0, upper=1> theta[N];
}
model {
# default prior
target += -5 * log(alpha+beta)/2;

# implicit joint distributions
theta ~ beta(alpha, beta) ;
y ~ binomial(n, theta);
}
"

An alternative to jointly sampling \(\theta, \alpha, \beta\) is to

1. sample \(\alpha, \beta \sim p(\alpha, \beta\mid y)\), and then
2. sample \(\theta \stackrel{ind}{\sim} p(\theta_i \mid \alpha, \beta, y_i) \stackrel{d}{=} Be(\alpha + y_i, \beta + n_i - y_i)\).

The marginal posterior for \(\alpha, \beta\) is \[ p(\alpha, \beta\mid y) \propto p(y\mid \alpha, \beta)p(\alpha, \beta) = \left[\prod_{i = 1}^n \text{Beta-Binomial}(y_i \mid n_i, \alpha, \beta) \right]p(\alpha, \beta) \]

# Marginalized (integrated) theta out of the model
model_marginalized = "
data {
int<lower=0> N;
int<lower=0> n[N];
int<lower=0> y[N];
}
parameters {
real<lower=0> alpha;
real<lower=0> beta;
}
model {
target += -5 * log(alpha+beta)/2;
y ~ beta_binomial(n, alpha, beta);
}
"

The forth approach:

# Stan Beta-Binomial
model_seasons = "
data {
int<lower=0> N; int<lower=0> n[N]; int<lower=0> y[N] ;
real<lower=0> a; real<lower=0> b; real<lower=0> C; real m;
}
parameters {
real<lower=0, upper=1> mu;
real<lower=0> eta;
}
transformed parameters {
real<lower=0> alpha;
real<lower=0> beta;
alpha = eta * mu;
beta = eta * (1-mu) ;
}
model {
mu ~ beta(a,b);
eta ~ lognormal(m,C);
y ~ beta_binomial(n, alpha, beta);
}
generated quantities {
real<lower=0, upper=1> theta[N] ;
for (i in 1:N) 
    theta [i] = beta_rng(alpha+y[i], beta + n[i]-y[i]);
}
"

Exchangeability: The set \(Y_1, Y_2, \ldots, Y_n\) is exchangeable if the joint probabilty \(p(y_1, \ldots, y_n)\) is invariant to permutation of the indices. That is, for any permutation \(\pi\), \(p(y_1, \ldots, y_n) = p(y_{\pi_1}, \ldots, y_{\pi_n})\).

Theorem 1: All independent and identically distributed random variables are exchangeable.

Theorem 2 (de Finetti’s Theorem): A sequence of random variables \((y_1, y_2, \ldots)\) is infinitely exchangeable iff, for all \(n\), \[ p(y_1, \ldots, y_n) = \int \prod_{i = 1}^n p(y_i\mid \theta)P(d\theta) \] for some measure \(P\) on \(\theta\).

Although hierarchical models are typically written using the conditional independence notation above, the assumptions underlying the model are exchangeability and functional forms for the priors. \[ y_i \stackrel{ind}{\sim} p(y\mid \theta_i) \quad \theta_i \stackrel{ind}{\sim} p(\theta\mid \phi) \quad \phi \sim p(\phi) \]

Example 2: Normal hierarchical model

Suppose we have the following model: \[ \begin{aligned} y_{ij} & \sim N(\theta_i, \sigma^2)\\ \theta_i & \sim N(\mu, \tau^2) \end{aligned} \] with \(j = 1,\ldots, n_i\), \(i = 1, \ldots, I\) and \(n = \sum_{i = 1}^I n_i\). We assume \(\sigma^2 = s^2\) is known and \(p(\mu, \tau^2) \propto p(\tau)\), i.e. assume an improper uniform prior on \(\mu\).

(Section 5.4 in Bayesian Data Analysis) Let \(\bar y_{i.} = \frac{1}{n_i} \sum_{i = 1}^{n_i} y_{ij}\) and \(s^2_i = s^2/n_{i}\), then \[ \left\{\begin{aligned} p(\tau \mid y) & \propto p(\tau) V_{\mu}^{1 / 2} \prod_{i=1}^{\mathrm{I}}\left(s_{i}^{2}+\tau^{2}\right)^{-1 / 2} \exp \left(-\frac{\left(\bar{y}_{i}-\hat{\mu}\right)^{2}}{2\left(s_{i}^{2}+\tau^{2}\right)}\right) \\ \mu \mid \tau, y & \sim N\left(\hat{\mu}, V_{\mu}\right) & \\ \theta_{i} \mid \mu, \tau, y & \sim N\left(\hat{\theta}_{i}, V_{i}\right) & \\ \end{aligned}\right. \] with \[ V_{\mu}^{-1} =\sum_{j=1}^{J} \frac{1}{s_{i}^{2}+\tau^{2}} \quad \hat{\mu}=V_{\mu}\left(\sum_{i=1}^{\mathrm{I}} \frac{\bar{y}_{i}}{s_{i}^{2}+\tau^{2}}\right) \\ V_{i}^{-1} =\frac{1}{s_{i}^{2}}+\frac{1}{\tau^{2}} \quad \hat{\theta}_{i} =V_{i}\left(\frac{\bar{y}_{i}}{s_{i}^{2}}+\frac{\mu}{\tau^{2}}\right) \] One choice for \(p(\tau)\) is \(Ca^+(0, 1)\)

There are some alternative distributions for \(\theta_i\):

• Heavy-tailed: \(\theta_i \sim t_{\nu}(\mu, \tau^2)\)
• Peak at zero: \(\theta_i \sim Laplace(\mu, \tau^2)\)
• Point mass at zero: \(\theta_i \sim \pi\delta_0 + (1 - \pi)N(\mu, \tau^2)\)

Hierarchical models

• allow the data to inform us about similarities across groups
• provide data driven shrinkage toward a grand mean
  • lots of shrinkage when means are similar
  • little shrinkage when means are different