Transcript Lecture 2 - UoB Interactive Server

Lecture 16
MCMC for Poisson response
models
Lecture Contents
•
•
•
•
•
Multilevel Poisson Model
MCMC Algorithms for Poisson Models
MLwiN for Poisson models
WinBUGS for Poisson models
Comparisons
Multilevel Poisson model
The Poisson distribution is used to model count data. It is
often used in disease data as an approximation to the
Binomial distribution where the exposure (number of
trials) is used as an offset variable.
In the practical we will look in detail at a veterinary
epidemiology dataset on TB cases in animal herds in
Canada. Here the response is the number of cases of
TB in the period 1985 to 1994 in herds of cattle, cervids
and bison. Note that the data contains only herds which
were infected in outbreaks of TB and the data has been
falsified to meet confidentiality regulations.
Multilevel Poisson Regression
A multilevel Poisson regression model can be written as
follows:
yij ~ Poisson(ij )
log(ij )  log(expij )  X ij   u j
u j ~ N (0,  u2 )
You will notice the offset (log(expij)) which in the TB
example is the number of animal days at risk in the
group. Including this offset allows comparison of rates
rather than number of cases which is more sensible as
otherwise the model will simply predict that more
cases are seen in herds with more days at risk!
Bayesian Multilevel Poisson model
To extent the model into a Bayesian framework we
need to include priors for the fixed effects and
between herd variance.
We will use the standard ‘diffuse’ priors as shown
below:
yij ~ Poisson(ij )
log(ij )  log(expij )  X ij   u j
u j ~ N (0,  u2 ),
  1,  u2 ~  1 ( ,  )
MCMC Algorithm
Our MLwiN algorithm has 3 steps:
• 1. Generate βi by Univariate Normal MH
Sampling in MLwiN or Gamerman method
in WinBUGS.
• 2. Generate each uj by Univariate Normal
MH Sampling in MLwiN or AR sampling in
WinBUGS.
• 3. Generate 1/σu2 from its Gamma
conditional distribution.
A final model for TB-Real ?
We will consider here only the
following model for the
dataset which contains all
feasible predictors.
The model to the right has been
fitted using 1st order MQL
estimation.
Note: to construct this model you
will need to make predictors
categorical via the Names
window and construct the
offset via the Calculate
window.
MLwiN MH Estimation
Using MH in MLwiN gives the following estimates after 50,000 iterations:
Trajectories plot
The (thinned) trajectories are as follows:
WinBUGS code
model
{
# Level 1 definition
for(i in 1:N) {
reactors[i] ~ dpois(mu[i])
log(mu[i]) <- offs[i] + beta[1] * cons[i]
+ beta[2] * type_2[i]
+ beta[3] * type_3[i]
+ beta[4] * type_5[i]
+ beta[5] * sex_2[i]
+ beta[6] * age_1[i]
+ beta[7] * age_2[i]
+ u2[farm_id[i]] * cons[i]
}
# Higher level definitions
for (j in 1:n2) {
u2[j] ~ dnorm(0,tau.u2)
}
# Priors for fixed effects
for (k in 1:7) { beta[k] ~ dflat() }
# Priors for random terms
tau.u2 ~ dgamma(0.001000,0.001000)
sigma2.u2 <- 1/tau.u2
}
Here we see the use of dpois
for a Poisson distribution.
Note that when code is
generated via MLwiN the
offset is always named ‘offs’.
Fortunately WinBUGS allows
variable names to include the
underscore character.
WinBUGS timing comparison
• WinBUGS took 2
minutes 4 seconds for
55,000 iterations
compared with 15
seconds in MLwiN.
• The chain for the
intercept and two
worst mixing fixed
effects were as
follows:
beta[1]
-8.0
-10.0
-12.0
-14.0
-16.0
5001
20000
40000
iteration
beta[6]
8.0
6.0
4.0
2.0
0.0
5001
20000
40000
iteration
beta[7]
6.0
4.0
2.0
0.0
5001
20000
40000
iteration
WinBUGS/MLwiN estimates
The following estimates were produced:
Parameter 1st MQL
MLwiN
β0
-10.071
-11.231 (-12.96, -9.78) -11.289 (-13.09, -9.68)
β1
-0.388
-0.405 (-1.082, 0.289)
-0.405 (-1.065, 0.248)
β2
-0.412
-0.241 (-1.205, 0.752)
-0.238 (-1.219, 0.721)
β3
-0.319
-0.250 (-2.019, 1.303)
-0.240 (-1.992, 1.322)
β4
-0.350
-0.352 (-0.760, 0.053)
-0.348 (-0.772, 0.051)
β5
2.715
2.873 (1.614, 4.503)
2.919 (1.569, 4.606)
β6
2.457
2.626 (1.400, 4.235)
2.671 (1.374, 4.331)
σ2u
2.904
2.016 (0.964, 3.883)
2.043 (0.964, 4.029)
WinBUGS
Effective sample size comparison
Note that the Ratio in ESS figures below have to be balanced by
the fact that MLwiN is 124/15 = 8.27 times faster!
Parameter
MLwiN
WinBUGS
Ratio
β0
66
215
3.26
β1
520
2829
5.44
β2
304
1328
4.37
β3
865
4362
5.04
β4
3666
22521
6.14
β5
77
237
3.08
β6
72
224
3.11
σ2u
3244
15498
4.78
Information for the practical
In the final practical of the week you are
asked to look at fitting more Poisson
models to the tb-real dataset.
You can use the DIC for Poisson models in
both MLwiN and WinBUGS to find firstly
the best single level model and then the
best random effects model.