Epilepsy: Repeated Measures on Poisson Counts

An example from OpenBUGS [44], Thall and Vail [92] Breslow and Clayton [10] concerning the effects of treatment, baseline seizure counts, and age on follow-up seizure counts at four visits in 59 patients.

Model

Counts are modelled as

y_{i,j} &\sim \text{Poisson}(\mu_{i,j}) \quad\quad i=1,\ldots,59; j=1,\ldots,4 \\
\log(\mu_{i,j}) &= \alpha_0 + \alpha_\text{Base} \log(\text{Base}_i / 4) +
  \alpha_\text{Trt} \text{Trt}_i + \alpha_\text{BT} \text{Trt}_i \log(\text{Base}_i / 4) + \\
  & \quad\quad \alpha_\text{Age} \log(\text{Age}_i) + \alpha_\text{V4} \text{V}_4 + \text{b1}_i +
    \text{b}_{i,j} \\
\text{b1}_i &\sim \text{Normal}(0, \sigma_\text{b1}) \\
\text{b}_{i,j} &\sim \text{Normal}(0, \sigma_\text{b}) \\
\alpha_* &\sim \text{Normal}(0, 100) \\
\sigma^2_\text{b1}, \sigma^2_\text{b} &\sim \text{InverseGamma}(0.001, 0.001),

where y_{ij} are the counts on patient i at visit j, \text{Trt} is a treatment indicator, \text{Base} is baseline seizure counts, \text{Age} is age in years, and \text{V}_4 is an indicator for the fourth visit.

Analysis Program

using Distributed
@everywhere using Mamba

## Data
epil = Dict{Symbol, Any}(
  :y =>
    [ 5  3  3  3
      3  5  3  3
      2  4  0  5
      4  4  1  4
      7 18  9 21
      5  2  8  7
      6  4  0  2
     40 20 21 12
      5  6  6  5
     14 13  6  0
     26 12  6 22
     12  6  8  4
      4  4  6  2
      7  9 12 14
     16 24 10  9
     11  0  0  5
      0  0  3  3
     37 29 28 29
      3  5  2  5
      3  0  6  7
      3  4  3  4
      3  4  3  4
      2  3  3  5
      8 12  2  8
     18 24 76 25
      2  1  2  1
      3  1  4  2
     13 15 13 12
     11 14  9  8
      8  7  9  4
      0  4  3  0
      3  6  1  3
      2  6  7  4
      4  3  1  3
     22 17 19 16
      5  4  7  4
      2  4  0  4
      3  7  7  7
      4 18  2  5
      2  1  1  0
      0  2  4  0
      5  4  0  3
     11 14 25 15
     10  5  3  8
     19  7  6  7
      1  1  2  3
      6 10  8  8
      2  1  0  0
    102 65 72 63
      4  3  2  4
      8  6  5  7
      1  3  1  5
     18 11 28 13
      6  3  4  0
      3  5  4  3
      1 23 19  8
      2  3  0  1
      0  0  0  0
      1  4  3  2],
  :Trt =>
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1],
  :Base =>
    [11, 11, 6, 8, 66, 27, 12, 52, 23, 10, 52, 33, 18, 42, 87, 50, 18, 111, 18,
     20, 12, 9, 17, 28, 55, 9, 10, 47, 76, 38, 19, 10, 19, 24, 31, 14, 11, 67,
     41, 7, 22, 13, 46, 36, 38, 7, 36, 11, 151, 22, 41, 32, 56, 24, 16, 22, 25,
     13, 12],
  :Age =>
    [31, 30, 25, 36, 22, 29, 31, 42, 37, 28, 36, 24, 23, 36, 26, 26, 28, 31, 32,
     21, 29, 21, 32, 25, 30, 40, 19, 22, 18, 32, 20, 30, 18, 24, 30, 35, 27, 20,
     22, 28, 23, 40, 33, 21, 35, 25, 26, 25, 22, 32, 25, 35, 21, 41, 32, 26, 21,
     36, 37],
  :V4 => [0, 0, 0, 1]
)
epil[:N] = size(epil[:y], 1)
epil[:T] = size(epil[:y], 2)

epil[:logBase4] = log.(epil[:Base] / 4)
epil[:BT] = epil[:logBase4] .* epil[:Trt]
epil[:logAge] = log.(epil[:Age])
map(key -> epil[Symbol(string(key, "bar"))] = mean(epil[key]),
    [:logBase4, :Trt, :BT, :logAge, :V4])


## Model Specification
model = Model(

  y = Stochastic(2,
    (a0, alpha_Base, logBase4, logBase4bar, alpha_Trt, Trt, Trtbar,
     alpha_BT, BT, BTbar, alpha_Age, logAge, logAgebar, alpha_V4, V4,
     V4bar, b1, b, N, T) ->
      UnivariateDistribution[(
        mu = exp(a0 + alpha_Base * (logBase4[i] - logBase4bar) +
                 alpha_Trt * (Trt[i] - Trtbar) + alpha_BT * (BT[i] - BTbar) +
                 alpha_Age * (logAge[i] - logAgebar) +
                 alpha_V4 * (V4[j] - V4bar) + b1[i] +
                 b[i, j]);
        Poisson(mu)) for i in 1:N, j in 1:T
      ],
    false
  ),

  b1 = Stochastic(1,
    s2_b1 -> Normal(0, sqrt(s2_b1)),
    false
  ),

  b = Stochastic(2,
    s2_b -> Normal(0, sqrt(s2_b)),
    false
  ),

  a0 = Stochastic(
    () -> Normal(0, 100),
    false
  ),

  alpha_Base = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_Trt = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_BT = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_Age = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_V4 = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha0 = Logical(
    (a0, alpha_Base, logBase4bar, alpha_Trt, Trtbar, alpha_BT, BTbar,
     alpha_Age, logAgebar, alpha_V4, V4bar) ->
      a0 - alpha_Base * logBase4bar - alpha_Trt * Trtbar - alpha_BT * BTbar -
        alpha_Age * logAgebar - alpha_V4 * V4bar
  ),

  s2_b1 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  s2_b = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)


## Initial Values
inits = [
  Dict(:y => epil[:y], :a0 => 0, :alpha_Base => 0, :alpha_Trt => 0,
       :alpha_BT => 0, :alpha_Age => 0, :alpha_V4 => 0, :s2_b1 => 1,
       :s2_b => 1, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T])),
  Dict(:y => epil[:y], :a0 => 1, :alpha_Base => 1, :alpha_Trt => 1,
       :alpha_BT => 1, :alpha_Age => 1, :alpha_V4 => 1, :s2_b1 => 10,
       :s2_b => 10, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T]))
]


## Sampling Scheme
scheme = [AMWG([:a0, :alpha_Base, :alpha_Trt, :alpha_BT, :alpha_Age,
                :alpha_V4], 0.1),
          Slice(:b1, 0.5),
          Slice(:b, 0.5),
          Slice([:s2_b1, :s2_b], 1.0)]
setsamplers!(model, scheme)


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

Results

Iterations = 2502:15000
Thinning interval = 2
Chains = 1,2
Samples per chain = 6250

Empirical Posterior Estimates:
               Mean         SD        Naive SE       MCSE        ESS
      s2_b  0.13523750 0.031819272 0.00028460022 0.0025401820 156.91007
     s2_b1  0.24911885 0.073166731 0.00065442313 0.0044942066 265.04599
  alpha_V4 -0.09287934 0.083666872 0.00074833925 0.0051087325 268.21356
 alpha_Age  0.45830900 0.394536219 0.00352883922 0.0310012419 161.96291
  alpha_BT  0.24217000 0.190566444 0.00170447809 0.0163585849 135.70673
 alpha_Trt -0.75931393 0.397734236 0.00355744316 0.0337796459 138.63592
alpha_Base  0.91104974 0.135354470 0.00121064718 0.0111438503 147.52807
    alpha0 -1.35617079 1.313240197 0.01174597740 0.1021568814 165.25442

Quantiles:
               2.5%        25.0%       50.0%        75.0%        97.5%
      s2_b  0.07155807  0.112591385  0.13626498  0.158032604  0.193715882
     s2_b1  0.13817484  0.197135457  0.23761145  0.289655043  0.422805121
  alpha_V4 -0.25504531 -0.148157017 -0.09313603 -0.036681431  0.072099011
 alpha_Age -0.19666987  0.176356196  0.41608686  0.696647899  1.305075377
  alpha_BT -0.09025014  0.108102968  0.22656174  0.358354280  0.657804969
 alpha_Trt -1.63682212 -1.011390894 -0.75653998 -0.480870874 -0.016113397
alpha_Base  0.66318818  0.817700638  0.90268210  0.997417378  1.200619714
    alpha0 -4.16888778 -2.157932918 -1.26343143 -0.436226488  0.866195785