The dispersion parameter in the quasi-Poisson model
Let us first see how the dispersion parameter is calculated in the model using quasi-Poisson likelihood. One assumption in Poisson regression is that the conditional variance and the conditional mean of the response $Y$ are the same:
$$
V(Y_{i}|\eta_{i})=E(Y_{i}|\eta_{i})=\mu_{i}
$$
where $\eta_{i}$ is the linear predictor. Sometimes we observe that the variance is greater than the mean which is called overdispersion. The quasi-Poisson likelihood model is a simple remedy for overdispersed count data because it introduces a dispersion parameter ($\phi$) into the Poisson model, so that the conditional variance of the response is now a linear function of the mean:
$$
V(Y_{i}|\eta_{i})=\phi\mu_{i}
$$
If $\phi>1$, the conditional variance increases more rapidly than its mean. Now, $\phi$ is estimated as:
$$
\widehat{\phi}=\frac{1}{n-k}\sum\frac{(Y_{i}-\hat{\mu_{i}})^2}{\hat{\mu_{i}}}
$$
where $n$ is the sample size, $k$ is the number of estimated parameters (including the intercept) and $\widehat{\mu_{i}}=g^{-1}(\widehat{\eta_{i}})$ is the fitted expectation of $Y_{i}$ ($g^{-1}$ is the inverse link function). Note that the formulation above is simply the Pearson $\chi^{2}$ divided by the residual degrees of freedom: $\widehat{\phi}=\chi^{2}/df$. Let's estimate the dispersion parameter in your example in R:
n <- 20
k <- 2
1/(n-k)*sum((C.OD-fitted(glm.fit.no.OD))^2/fitted(glm.fit.no.OD))
[1] 1.4703
# Via Pearson residuals
sum(residuals(glm.fit.no.OD, type="pearson")^2)/df.residual(glm.fit.no.OD)
[1] 1.4703
The estimated dispersion parameter is $1.47$. This is the same as given by the model output from the quasi-Poisson model:
summary(glm.fit.with.OD)
[...]
(Dispersion parameter for quasipoisson family taken to be 1.470301)
[...]
The dispersion parameter in the negative binomial model
There are several approaches to model count data with overdispersion. One very popular approach negative binomial regression. The conditional variance-mean relationship in the negative binomial model is:
$$
V(Y_{i}|\eta_{i})=\mu_{i}+\mu_{i}^{2}/\phi=\mu_{i}(1+\mu_{i}/\phi)
$$
where the second term provides the overdispersion where smaller $\phi$s denote stronger overdispersion. As $\phi\rightarrow\infty$, the variance approaches the mean and the distribution approaches the Poisson distribution as the second term gets very small. Sometimes, $1/\phi$ is used instead and if $1/\phi\rightarrow 0$, the distribution approaches the Poisson. In contrast to the quasi-Poisson model, the conditional variance is now a quadratic function of the mean (see this paper for more information). Note that the dispersion parameter is the multiplicative factor $1+\mu_{i}\phi$, which depends on $\mu_{i}$ (in contrast to the quasi-Poisson model).
We can fit a negative binomial regression using glm.nb from the MASS package:
library(MASS)
glm.negbin <- glm.nb(C.OD ~ x)
summary(glm.negbin)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4055 0.2847 1.424 0.154360
xarable 1.2622 0.3381 3.734 0.000189 ***
(Dispersion parameter for Negative Binomial(6.9569) family taken to be 1)
Theta: 6.96
Std. Err.: 6.74
The estimate for $\phi$ is $6.96$ (it's called theta in the output). The reciprocal of $\phi$ is sometimes used and is $1/6.96\approx0.144$ in this case. Because $1/\phi>0$ we have overdispersion.
Other methods to model overdispersed count data
I will outline three approaches here which can be found in the book Bayesian Modeling Using WinBUGS by Ioannis Ntzoufras (section 8.31, pages 283-286; section 9.2.3, pages 315-318).
Estimating overdispersion using a negative binomial model
Our model is:
$$
\begin{align}
Y_{i} &\sim \text{NB}(\pi_{i},r_{i}) \\
\pi_{i} &= \frac{r_{i}}{r_{i}+\lambda_{i}} \\
\log(\lambda_{i}) &= \beta_{0} + \sum_{j=1}^{p}\beta_{j}X_{ij}
\end{align}
$$
We have two groups, so the two lambdas are $\lambda_{1}=e^{\beta_{0}}, \lambda_{2}=e^{\beta_{0}+\beta_{1}}$. The dispersion index is given by $\text{DI}=1+\lambda/r$. If $\text{DI}>0$ we have overdispersion. Here is our WinBUGS/JAGS model (I used JAGS to sample from the posterior):
library(rjags)
library(R2jags)
sink("Negbin_model.txt")
cat("
model {
for (i in 1:n) {
y[i] ~ dnegbin(p.ind[i], r.ind[i])
p.ind[i] <- r.ind[i]/(r.ind[i] + lambda.ind[i])
log(lambda.ind[i]) <- beta[1] + beta[2]*x[i]
r.ind[i] <- r[ x[i] + 1 ]
}
lambda[1] <- exp(beta[1])
lambda[2] <- exp(beta[1] + beta[2])
beta1 <- exp(beta[2])
for(j in 1:2) {
logr.cont[j] ~ dunif(0, 10)
log(r.cont[j]) <- logr.cont[j]
r[j] <- round(r.cont[j])
#r[j] ~dgamma(0.001, 0.001)
beta[j] ~ dnorm(0.0, 0.0001)
DI[j] <- (1 + lambda[j])/r[j]
vari[j] <- lambda[j]*DI[j]
p[j] <- r[j]/(r[j]+lambda[j])
}
}
",fill=TRUE)
sink()
# Bundle data
win.data <- list(y = C.OD, x = as.numeric(x)-1, n = length(x))
# Inits function
inits <- function(){ list(beta=rlnorm(2), logr.cont=runif(2, 0,10))}
params <- c("beta", "lambda", "r", "DI", "p", "vari")
# MCMC settings
nc <- 3 # Number of chains
ni <- 50000 # Number of draws from posterior per chain
nb <- 10000 # Number of draws to discard as burn-in
nt <- 1 # Thinning rate
# Start Gibbs sampling
#out <- bugs(data=win.data, inits=inits, parameters.to.save=params,
# model.file="Negbin_model.txt", n.thin=nt, n.chains=nc,
# n.burnin=nb, n.iter=ni, debug = TRUE, program="OpenBUGS")
#print(out, dig = 3)
out <- jags(
data = win.data,
parameters.to.save = params,
model.file = "Negbin_model.txt",
n.chains = nc,
n.iter = ni,
n.burnin = nb,
n.thin=nt,
inits=inits,
progress.bar="text")
#out <- update(out, n.iter=50000)
out
The output is:
Inference for Bugs model at "Negbin_model.txt", fit using jags,
3 chains, each with 50000 iterations (first 10000 discarded)
n.sims = 120000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
DI[1] 0.141 0.347 0.000 0.001 0.010 0.098 1.166 1.001 59000
DI[2] 1.025 1.387 0.000 0.032 0.516 1.520 4.486 1.001 67000
beta[1] 0.378 0.269 -0.176 0.203 0.386 0.562 0.880 1.001 35000
beta[2] 1.287 0.326 0.664 1.066 1.280 1.502 1.950 1.001 22000
lambda[1] 1.512 0.404 0.839 1.225 1.471 1.755 2.412 1.001 35000
lambda[2] 5.381 1.035 3.657 4.704 5.285 5.933 7.700 1.001 27000
p[1] 0.942 0.111 0.585 0.945 0.994 0.999 1.000 1.001 120000
p[2] 0.682 0.271 0.205 0.441 0.699 0.973 1.000 1.001 33000
r[1] 2452.090 4589.758 2.000 25.000 240.000 2326.000 17602.050 1.001 59000
r[2] 1143.373 3351.097 2.000 4.000 12.000 193.000 13614.000 1.001 80000
vari[1] 0.232 0.696 0.000 0.002 0.015 0.143 1.923 1.001 57000
vari[2] 5.964 10.345 0.002 0.170 2.629 7.883 30.075 1.001 56000
deviance 83.371 2.319 80.238 81.701 82.964 84.510 89.135 1.001 120000
The posterior median dispersion index for the group "arable" (DI[2]) is larger than zero indicating overdispersion. On the other hand, the dispersion index for the group "grassland" (DI[1]) is only slightly larger than zero. Let's look at the posterior density of the dispersion index for the group "arable" and calculate the 95% Highest Posterior Density intervals (HDP) for the dispersion indices:
library(ggplot2)
library(runjags)
jagsfit.matrix <- rbind(as.matrix(as.mcmc(out)[[1]]),
as.matrix(as.mcmc(out)[[2]]),
as.matrix(as.mcmc(out)[[3]]))
name <- "DI[2]"
vect <- jagsfit.matrix[, name]
vect.plot <- vect[vect<=20]
mcmc.combined <- combine.mcmc(as.mcmc(out))
hpd.ints <- HPDinterval(mcmc.combined, prob=0.95)
hpd.ints
#hdr(mcmc.combined[,"DI"], prob=95, h=hdrbw(mcmc.combined[,"DI"], gridsize=1000000, HDRlevel=0.95), nn=5000)
plot.frame <- data.frame(dispersion=vect.plot)
ggplot(plot.frame, aes(x=vect.plot)) +
geom_density(alpha=0.5, fill="#1B4F97", color="#1B4F97") +
geom_vline(xintercept = c(0, 3.561465e+00), alpha=0.6, size=1) +
xlim(c(0,20)) +
ylab("Density") +
xlab("Dispersion index") +
ggtitle("Posterior distribution of the dispersion index for the group \"arable\"") +
theme(axis.title.y =element_text(vjust=0.4, size=20, angle=90)) +
theme(axis.title.x =element_text(vjust=0, size=20, angle=0)) +
theme(axis.text.x =element_text(size=15, colour = "black")) +
theme(axis.text.y =element_text(size=17, colour = "black")) +
theme(panel.background = element_rect(fill = "grey85", colour = NA),
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "grey90", size = 0.25))
The 95% HDP for the dispersion parameter for the group "arable" ranges from $0$ to $3.56$ (marked by the vertical grey lines in the graphic above).
Estimating overdispersion using a Poisson-log-normal model
First, let's define our model:
$$
\begin{align}
Y_{i} &\sim \text{Poisson}(\lambda_{i}) \\
\log(\lambda_{i}) &= \mu_{i}+b_{i} \\
\mu_{i} &= \beta_{0} + \sum_{j=1}^{p}\beta_{j}X_{ij} \\
b_{i} &\sim \mathcal{N}(0, \sigma^{2})
\end{align}
$$
The mean and the variance then are
$$
\begin{align}
E(Y|\lambda,\sigma_{b}^{2}) &= \lambda e^{\sigma_{b}^{2}/2} \\
V(Y|\lambda, \sigma_{b}^{2}) &= \lambda e^{\sigma_{b}^{2}/2} + \lambda^{2}e^{2\sigma_{b}^{2}} - \lambda^{2}e^{\sigma_{b}^{2}}
\end{align}
$$
So we simply add $b_{i}$ to the linear predictor to take the overdispersion into account.
We can build the dispersion index (DI) into the model and estimate it for each group separately:
sink("Poisson.OD.t.test.txt")
cat("
model {
# Priors
alpha ~ dnorm(0,0.001)
beta ~ dnorm(0,0.001)
sigma ~ dunif(0, 10)
sigma2 <- sigma*sigma
tau <- 1 / sigma2
maybe_overdisp <- mean(exp_eps[])
kappa[1] <- exp(alpha)
kappa[2] <- exp(alpha + beta)
mean.x[1] <- exp(alpha + 0.5*sigma2)
mean.x[2] <- exp(alpha + beta + 0.5*sigma2)
vari[1] <- kappa[1]*exp(sigma2/2)+kappa[1]*kappa[1]*exp(2*sigma2) -
kappa[1]*kappa[1]*exp(sigma2)
vari[2] <- kappa[2]*exp(sigma2/2)+kappa[2]*kappa[2]*exp(2*sigma2) -
kappa[2]*kappa[2]*exp(sigma2)
DI[1] <- vari[1]/mean.x[1]
DI[2] <- vari[2]/mean.x[2]
# Likelihood
for (i in 1:n) {
C.OD[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha + beta*x[i] + eps[i]
eps[i] ~ dnorm(0, tau)
exp_eps[i] <- exp(eps[i])
#di.index.ind[i] <- 1 + exp(eps[i])*(es-1)*sqrt(es)
}
}
",fill=TRUE)
sink()
# Bundle data
win.data <- list(C.OD = C.OD, x = as.numeric(x)-1, n = length(x))
# Inits function
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1), sigma=rlnorm(1))}
# Parameters to estimate
params <- c("alpha", "beta", "sigma", "sigma2",
"maybe_overdisp", "DI", "mean.x", "vari", "kappa")
# MCMC settings
nc <- 3 # Number of chains
ni <- 50000 # Number of draws from posterior per chain
nb <- 10000 # Number of draws to discard as burn-in
nt <- 1 # Thinning rate
# Start Gibbs sampling
out2 <- jags(
data = win.data,
parameters.to.save = params,
model.file = "Poisson.OD.t.test.txt",
n.chains = nc,
n.iter = ni,
n.burnin = nb,
n.thin=nt,
inits=inits,
progress.bar="text")
The output is
Inference for Bugs model at "Poisson.OD.t.test.txt", fit using jags,
3 chains, each with 50000 iterations (first 10000 discarded)
n.sims = 120000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
DI[1] 1.746 16.499 1.004 1.132 1.312 1.666 4.054 1.001 100000
DI[2] 3.573 22.604 1.013 1.486 2.133 3.404 12.319 1.001 8500
alpha 0.273 0.324 -0.416 0.069 0.293 0.496 0.857 1.001 4800
beta 1.294 0.384 0.559 1.042 1.286 1.538 2.076 1.001 7400
kappa[1] 1.382 0.436 0.660 1.072 1.340 1.642 2.355 1.001 4800
kappa[2] 4.911 1.076 2.927 4.199 4.868 5.560 7.167 1.001 23000
maybe_overdisp 1.119 0.176 0.900 1.006 1.075 1.187 1.570 1.001 49000
mean.x[1] 1.581 0.565 0.794 1.220 1.506 1.847 2.802 1.001 7000
mean.x[2] 5.632 1.505 3.591 4.728 5.413 6.225 9.022 1.001 31000
sigma 0.461 0.232 0.050 0.300 0.442 0.601 0.970 1.009 1300
sigma2 0.266 0.260 0.002 0.090 0.195 0.361 0.940 1.009 1300
vari[1] 6.195 846.592 0.925 1.501 1.986 2.844 9.543 1.001 45000
vari[2] 32.314 1013.170 4.830 7.456 11.123 19.352 98.875 1.001 11000
deviance 74.771 5.949 64.349 70.338 74.364 78.967 86.505 1.001 5100
First note the posterior medians of the groups means.x are close to the observed means of the data ($1.51$ vs. $1.5$ and $5.41$ vs. $5.3$). The posterior variances vari for the groups are also close to the observed variances ($1.99$ vs. $1.39$ and $11.12$ vs. $10.68$). Importantly, the posterior median of the estimate $\hat{\beta_{1}}$ (beta) ($1.286$) is very close to the estimate calculated by glm using family="quasipoisson" which was $1.262$. The posterior mean and median of the dispersion index (DI) is $1.75$ and $1.31$ for the group "grassland" and $3.57$ and $2.13$ for the group "arable". It seems that the data for the group "arable" is more overdispersed than the group "grassland". The dispersion parameter estimated by glm with quasi-Poisson likelihood was around $1.47$ which is in between the posterior medians of the two dispersion indices, so our estimations look reasonable. Let's look at the posterior density of the dispersion index for the group "arable" and calculate the 95% Highest Posterior Density intervals (HDP) for the dispersion indices:
library(ggplot2)
library(runjags)
jagsfit.matrix <- rbind(as.matrix(as.mcmc(out2)[[1]]),
as.matrix(as.mcmc(out2)[[2]]),
as.matrix(as.mcmc(out2)[[3]]))
name <- "DI[2]"
vect <- jagsfit.matrix[, name]
vect.plot <- vect[vect<=20]
mcmc.combined <- combine.mcmc(as.mcmc(out2))
hpd.ints <- HPDinterval(mcmc.combined, prob=0.95)
hpd.ints
lower upper
DI[1] 1.000000e+00 3.0691670
DI[2] 1.000001e+00 8.5232449
alpha -3.729040e-01 0.8911411
beta 5.598711e-01 2.0766436
deviance 6.397199e+01 86.0493853
kappa[1] 5.908833e-01 2.2448459
kappa[2] 2.821622e+00 7.0389753
maybe_overdisp 8.565536e-01 1.4860521
mean.x[1] 6.856225e-01 2.5876715
mean.x[2] 3.300125e+00 8.3619700
sigma 4.044946e-04 0.8686939
sigma2 9.929906e-08 0.7545791
vari[1] 6.201473e-01 6.6375401
vari[2] 3.416068e+00 61.4973745
attr(,"Probability")
[1] 0.95
plot.frame <- data.frame(dispersion=vect.plot)
ggplot(plot.frame, aes(x=vect.plot)) +
geom_density(alpha=0.5, fill="#1B4F97", color="#1B4F97") +
geom_vline(xintercept = c(1, 8.5232449), alpha=0.6, size=1) +
xlim(c(0,20)) +
ylab("Density") +
xlab("Dispersion index") +
ggtitle("Posterior distribution of the dispersion index for the group \"arable\"") +
theme(axis.title.y =element_text(vjust=0.4, size=20, angle=90)) +
theme(axis.title.x =element_text(vjust=0, size=20, angle=0)) +
theme(axis.text.x =element_text(size=15, colour = "black")) +
theme(axis.text.y =element_text(size=17, colour = "black")) +
theme(panel.background = element_rect(fill = "grey85", colour = NA),
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "grey90", size = 0.25))
The 95% HDP for the dispersion parameter for the group "arable" ranges from $1$ to $8.52$ (marked by the vertical grey lines in the graphic above). It is important to note that we used a uniform prior for the standard deviation of $b_{i}$. There are other possibilities and the posterior distribution of the dispersion index can vary depending on the prior. Other priors include but are not limited to: uniform on the variance ($\sigma^{2}$), half-normal prior on $\sigma$ or $\sigma^{2}$, half-Cauchy on $\sigma$ and others.
Estimating overdispersion using a Poisson-gamma model
We can also model the data using a Poisson-gamma model:
$$
\begin{align}
Y_{i} &\sim \text{Poisson}(\lambda_{i}u_{i}) \\
u_{i} &\sim \text{Gamma}(r_{i}, r_{i})
\end{align}
$$
The WinBUGS model (or OpenBUGS, JAGS) is as follows:
sink("gamma_mix.txt")
cat("
model{
for(i in 1:n){
y[i] ~ dpois(mu.ind[i])
mu.ind[i] <- mu[i]*u[i]
log(mu[i]) <- beta[1]+beta[2]*x[i]
u[i] ~ dgamma(r[x[i]+1], r[x[i]+1])
}
mean.u <- mean(u[])
lambda[1] <- exp(beta[1])
lambda[2] <- exp(beta[1] + beta[2])
for (j in 1:2){
r[j] ~ dgamma(0.001, 0.001)
beta[j] ~ dnorm(0.0, 0.0001)
DI[j] <- (1+lambda[j]/r[j])
vari[j] <- lambda[j]*DI[j]
p[j] <- r[j]/(r[j]+lambda[j])
}
}
",fill=TRUE)
sink()
# Bundle data
win.data <- list(y = C.OD, x = as.numeric(x)-1, n = length(x))
# Inits function
inits <- function(){ list(beta=rlnorm(2), r=rlnorm(2))}
# Parameters to estimate
params <- c("beta", "lambda", "r", "DI", "mean.x", "vari", "tau", "s", "s2", "mean.u")
# MCMC settings
nc <- 3 # Number of chains
ni <- 50000 # Number of draws from posterior per chain
nb <- 10000 # Number of draws to discard as burn-in
nt <- 1 # Thinning rate
out <- jags(
data = win.data,
parameters.to.save = params,
model.file = "gamma_mix.txt",
n.chains = nc,
n.iter = ni,
n.burnin = nb,
n.thin=nt,
inits=inits,
progress.bar="text")
The output is
Inference for Bugs model at "gamma_mix.txt", fit using jags,
3 chains, each with 50000 iterations (first 10000 discarded)
n.sims = 120000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
DI[1] 1.180 0.492 1.001 1.008 1.038 1.163 2.196 1.002 5600
DI[2] 2.280 1.734 1.017 1.248 1.774 2.665 6.580 1.001 8500
beta[1] 0.375 0.286 -0.206 0.189 0.383 0.567 0.909 1.002 3100
beta[2] 1.292 0.352 0.612 1.058 1.287 1.523 1.995 1.002 3500
lambda[1] 1.515 0.445 0.814 1.209 1.467 1.763 2.481 1.002 3100
lambda[2] 5.414 1.180 3.524 4.649 5.287 6.009 8.087 1.001 12000
mean.u 1.001 0.093 0.820 0.952 0.998 1.045 1.207 1.001 34000
r[1] 184.429 373.385 1.358 8.872 37.503 172.345 1378.710 1.002 2100
r[2] 37.280 103.862 1.076 3.162 6.733 20.971 303.187 1.001 28000
vari[1] 1.861 2.106 0.863 1.294 1.600 1.998 4.259 1.001 27000
vari[2] 13.170 17.311 4.585 6.442 9.057 14.297 45.626 1.001 6500
deviance 76.089 5.763 65.979 71.679 75.802 80.428 87.121 1.001 16000
Again, the posterior means (lambda) and variances (vari) are very close to the observed ones. The estimate for the coefficient beta[2] is again practically identical to the estimates we've got using the Poisson-log-normal approach (i.e. $\approx 1.29$). The dispersion indices for the two groups are about $1.04$ for "grassland" and $1.77$ for "arable". These are very close to the observed overdispersion which are $0.926$ for "grassland" and $2.01$ for "arable", respectively. The HDPs and the posterior density of the dispersion index for the group "arable" is
jagsfit.matrix <- rbind(as.matrix(as.mcmc(out)[[1]]),
as.matrix(as.mcmc(out)[[2]]),
as.matrix(as.mcmc(out)[[3]]))
name <- "DI[2]"
vect <- jagsfit.matrix[, name]
vect.plot <- vect[vect<=20]
mcmc.combined <- combine.mcmc(as.mcmc(out))
hpd.ints <- HPDinterval(mcmc.combined, prob=0.95)
hpd.ints
lower upper
DI[1] 1.0002195 1.7847379
DI[2] 1.0022835 5.1658954
beta[1] -0.1984061 0.9146267
beta[2] 0.5995819 1.9807803
deviance 65.7749045 86.8262578
lambda[1] 0.7592323 2.3828823
lambda[2] 3.2764333 7.6725830
mean.u 0.8060325 1.1907674
r[1] 0.1508372 900.3718394
r[2] 0.2372253 193.6838708
vari[1] 0.6524941 3.2896264
vari[2] 3.5796354 33.4126294
attr(,"Probability")
[1] 0.95
plot.frame <- data.frame(dispersion=vect.plot)
ggplot(plot.frame, aes(x=vect.plot)) +
geom_density(alpha=0.5, fill="#1B4F97", color="#1B4F97") +
geom_vline(xintercept = c(1, 5.1658954), alpha=0.6, size=1) +
xlim(c(0,20)) +
ylab("Density") +
xlab("Dispersion index") +
ggtitle("Posterior distribution of the dispersion index for the group \"arable\"") +
theme(axis.title.y =element_text(vjust=0.4, size=20, angle=90)) +
theme(axis.title.x =element_text(vjust=0, size=20, angle=0)) +
theme(axis.text.x =element_text(size=15, colour = "black")) +
theme(axis.text.y =element_text(size=17, colour = "black")) +
theme(panel.background = element_rect(fill = "grey85", colour = NA),
panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "grey90", size = 0.25))
The 95% HDP for the dispersion parameter for the group "arable" ranges from $1$ to $5.17$ (marked by the vertical grey lines in the graphic above). The interval is smaller than the interval obtained by the Poisson-log-normal approach, which was ranging from $1$ to $8.52$.
Calculate dispersion parameter as in the quasi-Poisson model
The dispersion parameter in the quasi-Poisson GLM is estiamted as follows:
$$
\widehat{\mathrm{DI}}=\frac{1}{n-k}\sum_{i}^{n}r_{P,i}^{2}
$$
where $n$ is the sample size, $k$ the number of estimated parameters and $r_{P,i}$ are the Pearson residuals:
$$
r_{P}=\frac{y-\mu}{\sqrt{\mu}}
$$
The dispersion index can be estimated using normal Poisson regression in WinBUGS:
sink("poisson_dispersion.txt")
cat("
model{
for(i in 1:n){
y[i] ~ dpois(mu[i])
log(mu[i]) <- beta[1] + beta[2]*x[i]
fitted.y[i] <- exp(beta[1]+beta[2]*x[i])
}
DI.index <- 1/(n-2)*sum(pow((y[]-fitted.y[]),2)/fitted.y[])
for (j in 1:2){
beta[j] ~ dnorm(0.0, 0.0001)
}
}
",fill=TRUE)
sink()
# Bundle data
win.data <- list(y = C.OD, x = as.numeric(x)-1, n = length(x))
# Inits function
inits <- function(){ list(beta=rlnorm(2))}
# Parameters to estimate
params <- c("beta", "DI.index")
# MCMC settings
nc <- 3 # Number of chains
ni <- 50000 # Number of draws from posterior per chain
nb <- 10000 # Number of draws to discard as burn-in
nt <- 1 # Thinning rate
out <- jags(
data = win.data,
parameters.to.save = params,
model.file = "poisson_lognormal.txt",
n.chains = nc,
n.iter = ni,
n.burnin = nb,
n.thin=nt,
inits=inits,
progress.bar="text")
out
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
DI.index 1.637 0.272 1.354 1.441 1.560 1.750 2.359 1.001 4500
beta[1] 0.375 0.261 -0.170 0.208 0.388 0.556 0.855 1.002 2700
beta[2] 1.282 0.294 0.729 1.080 1.273 1.474 1.886 1.001 4300
deviance 84.491 2.068 82.536 83.054 83.869 85.265 89.928 1.001 44000
The posterior median of the dispersion index is 1.56 and the 95%-HDI is ranging from $1.345$ to $2.172$ and the value of $1.47$ as estimated by glm is well within the 95%-HDI. Heres the density plot of the posterior distribution of the dispersion index:
