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
where are the counts on patient at visit , is a treatment indicator, is baseline seizure counts, is age in years, and is an indicator for the fourth visit.
Analysis Program¶
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[
begin
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)
end
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