hierarchical-models

Hierarchical Models

When to Use

• Nested/grouped data (students in schools, patients in hospitals)
• Repeated measurements on subjects
• Meta-analysis with study-level variation
• Partial pooling between complete pooling and no pooling

Core Concept: Partial Pooling

Group means shrink toward overall mean based on:
- Within-group sample size
- Within-group variance
- Between-group variance

Stan Implementation

Centered Parameterization (Default)

data {
  int<lower=0> N;           // Total observations
  int<lower=0> J;           // Number of groups
  array[N] int<lower=1,upper=J> group;
  vector[N] y;
}
parameters {
  real mu;                  // Population mean
  real<lower=0> tau;        // Between-group SD
  real<lower=0> sigma;      // Within-group SD
  vector[J] theta;          // Group means
}
model {
  // Hyperpriors
  mu ~ normal(0, 10);
  tau ~ cauchy(0, 2.5);
  sigma ~ exponential(1);

  // Group effects
  theta ~ normal(mu, tau);

  // Likelihood
  y ~ normal(theta[group], sigma);
}

Non-Centered Parameterization (Better for weak data/small tau)

parameters {
  real mu;
  real<lower=0> tau;
  real<lower=0> sigma;
  vector[J] theta_raw;      // Standard normal
}
transformed parameters {
  vector[J] theta = mu + tau * theta_raw;
}
model {
  theta_raw ~ std_normal();
  // ... rest same
}

When to use non-centered: Divergences, small tau, few observations per group.

JAGS Implementation

model {
  for (i in 1:N) {
    y[i] ~ dnorm(theta[group[i]], tau.y)
  }

  for (j in 1:J) {
    theta[j] ~ dnorm(mu, tau.theta)
  }

  # Hyperpriors
  mu ~ dnorm(0, 0.0001)
  tau.theta <- pow(sigma.theta, -2)
  sigma.theta ~ dunif(0, 100)
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif(0, 100)
}

Classic Example: Eight Schools

data {
  int<lower=0> J;
  array[J] real y;          // Observed effects
  array[J] real<lower=0> sigma;  // Known SEs
}
parameters {
  real mu;
  real<lower=0> tau;
  vector[J] theta_raw;
}
transformed parameters {
  vector[J] theta = mu + tau * theta_raw;
}
model {
  mu ~ normal(0, 5);
  tau ~ cauchy(0, 5);
  theta_raw ~ std_normal();
  y ~ normal(theta, sigma);
}