Rats: A Normal Hierarchical Model

An example from OpenBUGS [38] and section 6 of Gelfand et al. [25] concerning 30 rats whose weights were measured at each of five consecutive weeks.

Model

Weights are modeled as

y_{i,j} &\sim \text{Normal}\left(\alpha_i + \beta_i (x_j - \bar{x}), \sigma_c\right) \quad\quad i=1,\ldots,30; j=1,\ldots,5 \\
\alpha_i &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\
\beta_i &\sim \text{Normal}(\mu_\beta, \sigma_\beta) \\
\mu_\alpha, \mu_\beta &\sim \text{Normal}(0, 1000) \\
\sigma^2_\alpha, \sigma^2_\beta, \sigma^2_c &\sim \text{InverseGamma}(0.001, 0.001),

where y_{i,j} is repeated weight measurement j on rat i, and x_j is the day on which the measurement was taken.

Analysis Program

using Mamba

## Data
rats = (Symbol => Any)[
  :y =>
    [151, 199, 246, 283, 320,
     145, 199, 249, 293, 354,
     147, 214, 263, 312, 328,
     155, 200, 237, 272, 297,
     135, 188, 230, 280, 323,
     159, 210, 252, 298, 331,
     141, 189, 231, 275, 305,
     159, 201, 248, 297, 338,
     177, 236, 285, 350, 376,
     134, 182, 220, 260, 296,
     160, 208, 261, 313, 352,
     143, 188, 220, 273, 314,
     154, 200, 244, 289, 325,
     171, 221, 270, 326, 358,
     163, 216, 242, 281, 312,
     160, 207, 248, 288, 324,
     142, 187, 234, 280, 316,
     156, 203, 243, 283, 317,
     157, 212, 259, 307, 336,
     152, 203, 246, 286, 321,
     154, 205, 253, 298, 334,
     139, 190, 225, 267, 302,
     146, 191, 229, 272, 302,
     157, 211, 250, 285, 323,
     132, 185, 237, 286, 331,
     160, 207, 257, 303, 345,
     169, 216, 261, 295, 333,
     157, 205, 248, 289, 316,
     137, 180, 219, 258, 291,
     153, 200, 244, 286, 324],
  :x => [8.0, 15.0, 22.0, 29.0, 36.0]
]
rats[:xbar] = mean(rats[:x])
rats[:N] = size(rats[:y], 1)
rats[:T] = size(rats[:y], 2)

rats[:rat] = Integer[div(i - 1, 5) + 1 for i in 1:150]
rats[:week] = Integer[(i - 1) % 5 + 1 for i in 1:150]
rats[:X] = rats[:x][rats[:week]]
rats[:Xm] = rats[:X] - rats[:xbar]


## Model Specification

model = Model(

  y = Stochastic(1,
    @modelexpr(alpha, beta, rat, Xm, s2_c,
      begin
        mu = alpha[rat] + beta[rat] .* Xm
        MvNormal(mu, sqrt(s2_c))
      end
    ),
    false
  ),

  alpha = Stochastic(1,
    @modelexpr(mu_alpha, s2_alpha,
      Normal(mu_alpha, sqrt(s2_alpha))
    ),
    false
  ),

  alpha0 = Logical(
    @modelexpr(mu_alpha, xbar, mu_beta,
      mu_alpha - xbar * mu_beta
    )
  ),

  mu_alpha = Stochastic(
    :(Normal(0.0, 1000)),
    false
  ),

  s2_alpha = Stochastic(
    :(InverseGamma(0.001, 0.001)),
    false
  ),

  beta = Stochastic(1,
    @modelexpr(mu_beta, s2_beta,
      Normal(mu_beta, sqrt(s2_beta))
    ),
    false
  ),

  mu_beta = Stochastic(
    :(Normal(0.0, 1000))
  ),

  s2_beta = Stochastic(
    :(InverseGamma(0.001, 0.001)),
    false
  ),

  s2_c = Stochastic(
    :(InverseGamma(0.001, 0.001))
  )

)


## Initial Values
inits = [
  [:y => rats[:y], :alpha => fill(250, 30), :beta => fill(6, 30),
   :mu_alpha => 150, :mu_beta => 10, :s2_c => 1, :s2_alpha => 1,
   :s2_beta => 1],
  [:y => rats[:y], :alpha => fill(20, 30), :beta => fill(0.6, 30),
   :mu_alpha => 15, :mu_beta => 1, :s2_c => 10, :s2_alpha => 10,
   :s2_beta => 10]
]


## Sampling Scheme
scheme = [Slice([:s2_c], [10.0]),
          AMWG([:alpha], fill(100.0, 30)),
          Slice([:mu_alpha, :s2_alpha], [100.0, 10.0], :univar),
          AMWG([:beta], ones(30)),
          Slice([:mu_beta, :s2_beta], [1.0, 1.0], :univar)]
setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, rats, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

Results

Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750

Empirical Posterior Estimates:
           Mean       SD      Naive SE     MCSE       ESS
mu_beta   6.18623 0.1072236 0.001238112 0.00215008 2486.9858
 alpha0 106.57458 3.6532621 0.042184237 0.05738100 4053.4550
   s2_c  37.01996 5.5173157 0.063708474 0.18973116  845.6260

Quantiles:
          2.5%     25.0%     50.0%     75.0%     97.5%
mu_beta  5.97845   6.11535   6.18663   6.25654   6.40124
 alpha0 99.36560 104.13484 106.54939 109.01257 113.78295
   s2_c 27.69329  33.09522  36.58565  40.29530  49.37235