7
$\begingroup$

I try to model quasi-poisson family in bugs language, to handle overdispersion. According to Introduction to WinBUGS for ecologists, this is done by:

$log(\lambda_i) = f(x_i) + \epsilon_i$

$N_i \sim Poiss(\lambda_i)$

where $\epsilon_i$ does the overdispersion. This is equivalent to scenario 8 in Lindén & Mäntyniemi 2011:

$log(\lambda_i) = f(x_i)$

$N_i \sim Poiss(\lambda_i \theta_i)$

However, how do I get the "dispersion parameter for quasipoisson family", as reported by GLM when using quasipoisson family?

> summary(glm.fit.with.OD)
Call:
glm(formula = C.OD ~ x, family = quasipoisson)
[...]
(Dispersion parameter for quasipoisson family taken to be 1.470301)

Here I include full reproducible code, taken from Introduction to WinBUGS for ecologists, see chapter 14.1 (slightly modified):

set.seed(123)

### 14.1. Overdispersion


### 14.1.2. Data generation
n.site <- 10
x <- gl(n = 2, k = n.site, labels = c("grassland", "arable"))
eps <- rnorm(2*n.site, mean = 0, sd = 0.5)# Normal random effect
lambda.OD <- exp(0.69 +(0.92*(as.numeric(x)-1) + eps) )
lambda.Poisson <- exp(0.69 +(0.92*(as.numeric(x)-1)) ) # For comparison

C.OD <- rpois(n = 2*n.site, lambda = lambda.OD)
C.Poisson <- rpois(n = 2*n.site, lambda = lambda.Poisson)

par(mfrow = c(1,2))
boxplot(C.OD ~ x, col = "grey", xlab = "Land-use", main = "With OD", 
ylab = "Hare count", las = 1, ylim = c(0, max(C.OD)))
boxplot(C.Poisson ~ x, col = "grey", xlab = "Land-use", main = "Without OD", 
ylab = "Hare count", las = 1, ylim = c(0, max(C.OD)) )


### 14.1.3. Analysis using R
glm.fit.no.OD <- glm(C.OD ~ x, family = poisson)
glm.fit.with.OD <- glm(C.OD ~ x, family = quasipoisson)
summary(glm.fit.no.OD)
summary(glm.fit.with.OD)
anova(glm.fit.no.OD, test = "Chisq")
anova(glm.fit.with.OD, test = "F")


### 14.1.4. Analysis using WinBUGS
# Define model
sink("Poisson.OD.t.test.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)   
 tau <- 1 / (sigma * sigma)
 maybe_overdisp <- mean(exp_eps[])

# 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])
 }
}
",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("lambda","alpha", "beta", "sigma", "maybe_overdisp")

# MCMC settings
nc <- 3     # Number of chains
ni <- 30000     # Number of draws from posterior per chain
nb <- 10000     # Number of draws to discard as burn-in
nt <- 5     # Thinning rate

# Start Gibbs sampling
out <- bugs(data=win.data, inits=inits, parameters.to.save=params, 
model.file="Poisson.OD.t.test.txt", n.thin=nt, n.chains=nc, 
n.burnin=nb, n.iter=ni, debug = TRUE)

print(out, dig = 3)

Here is the output of the glm() in R:

> summary(glm.fit.with.OD)

Call:
glm(formula = C.OD ~ x, family = quasipoisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2945  -0.5905  -0.4348   1.0798   1.8159  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.4055     0.3131   1.295  0.21166   
xarable       1.2622     0.3546   3.559  0.00224 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 1.470301)

    Null deviance: 51.032  on 19  degrees of freedom
Residual deviance: 28.524  on 18  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

And here is the output of the jags model (using more iterations than in the above code):

> print(out, dig = 3)
Inference for Bugs model at "Poisson.OD.t.test.txt", fit using jags,
 3 chains, each with 300000 iterations (first 100000 discarded), n.thin = 5
 n.sims = 120000 iterations saved
               mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat  n.eff
alpha            0.275   0.324 -0.415  0.073  0.291  0.498  0.861 1.001  23000
beta             1.289   0.385  0.547  1.035  1.281  1.534  2.071 1.001  41000
lambda[1]        1.163   0.586  0.263  0.740  1.092  1.496  2.512 1.001  24000
lambda[2]        1.368   0.654  0.398  0.914  1.273  1.700  2.934 1.001  27000
lambda[3]        2.173   1.081  0.818  1.429  1.919  2.638  4.976 1.001 120000
lambda[4]        1.367   0.654  0.397  0.912  1.273  1.701  2.927 1.001  26000
lambda[5]        1.370   0.655  0.401  0.914  1.275  1.702  2.929 1.001  20000
lambda[6]        1.602   0.759  0.545  1.086  1.464  1.954  3.503 1.001  24000
lambda[7]        1.365   0.656  0.399  0.909  1.268  1.697  2.941 1.001  34000
lambda[8]        1.366   0.653  0.394  0.912  1.272  1.701  2.919 1.001  43000
lambda[9]        1.368   0.654  0.400  0.916  1.270  1.700  2.934 1.001  16000
lambda[10]       1.872   0.904  0.690  1.260  1.679  2.261  4.198 1.001  84000
lambda[11]       4.608   1.548  1.978  3.525  4.494  5.521  8.032 1.001  54000
lambda[12]       5.058   1.628  2.340  3.937  4.903  5.969  8.818 1.001  52000
lambda[13]       6.544   2.064  3.426  5.091  6.189  7.649 11.521 1.001 120000
lambda[14]       4.170   1.492  1.610  3.095  4.078  5.116  7.355 1.001  67000
lambda[15]       4.613   1.554  1.974  3.524  4.499  5.528  8.073 1.001 120000
lambda[16]       7.066   2.258  3.759  5.439  6.675  8.293 12.477 1.001  49000
lambda[17]       7.615   2.462  4.054  5.802  7.190  9.003 13.472 1.001  42000
lambda[18]       3.385   1.458  0.956  2.281  3.264  4.385  6.403 1.001  61000
lambda[19]       6.538   2.066  3.431  5.089  6.185  7.643 11.505 1.001 120000
lambda[20]       3.385   1.459  0.957  2.279  3.260  4.384  6.399 1.001  95000
maybe_overdisp   1.121   0.179  0.898  1.006  1.076  1.189  1.577 1.001 120000
sigma            0.465   0.236  0.051  0.301  0.449  0.606  0.983 1.001  27000
deviance        74.767   5.985 64.292 70.311 74.343 79.042 86.528 1.001 120000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 17.9 and DIC = 92.7
DIC is an estimate of expected predictive error (lower deviance is better).

Edit 2

Modified model:

model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)   
 tau <- 1 / (sigma * sigma)
 maybe_overdisp <- mean(exp_eps[])

# co mi poradili na stats.SE
sigma2 <- sigma*sigma
kappa <- exp(alpha + beta); 
mean.x <- exp(alpha + 0.5*sigma2); 
vari <- kappa*exp(sigma2/2)+kappa*kappa*exp(2*sigma2) - kappa*kappa*exp(sigma2) 
DI <- vari/mean.x

# 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])
 }
}

output:

> print(out, dig = 6)
Inference for Bugs model at "Poisson.OD.t.test.txt", fit using jags,
 3 chains, each with 30000 iterations (first 10000 discarded), n.thin = 5
 n.sims = 12000 iterations saved
                 mu.vect    sd.vect      2.5%       25%       50%       75%     97.5%     Rhat n.eff
DI             20.691012 574.969428  2.748532  4.908448  7.479553 13.014255 64.410436 1.001834  2100
alpha           0.276339   0.319067 -0.410506  0.081158  0.293487  0.491261  0.859875 1.001180  6500
beta            1.285987   0.378753  0.556607  1.042379  1.276188  1.523705  2.054021 1.001084  9500
lambda[1]       1.173864   0.584403  0.268459  0.749624  1.106648  1.508932  2.499185 1.001344  4300
lambda[2]       1.364399   0.653852  0.395182  0.913774  1.272952  1.693450  2.944849 1.001782  2200
lambda[3]       2.148461   1.056564  0.814731  1.422458  1.892689  2.620223  4.912488 1.001393  3900
lambda[4]       1.373973   0.651289  0.404761  0.922161  1.283126  1.707821  2.930640 1.001074 10000
lambda[5]       1.364438   0.646131  0.397013  0.915541  1.271812  1.700321  2.908050 1.000922 12000
lambda[6]       1.596684   0.737480  0.562846  1.098588  1.463223  1.931152  3.472518 1.001139  7500
lambda[7]       1.372055   0.667968  0.394750  0.912321  1.270109  1.702314  3.030247 1.001004 12000
lambda[8]       1.368805   0.641905  0.393292  0.926187  1.281149  1.694738  2.933341 1.001159  7000
lambda[9]       1.367516   0.644055  0.410916  0.914202  1.279055  1.693786  2.892084 1.001237  5500
lambda[10]      1.874077   0.903156  0.706316  1.261625  1.683053  2.256230  4.224203 1.001380  4000
lambda[11]      4.604620   1.551257  1.983965  3.515142  4.484644  5.530132  8.040010 1.000890 12000
lambda[12]      5.052838   1.631220  2.325490  3.949769  4.904707  5.942355  8.784829 1.000883 12000
lambda[13]      6.506786   2.050675  3.424231  5.043272  6.159473  7.621714 11.341993 1.000896 12000
lambda[14]      4.168782   1.468719  1.617240  3.104178  4.099795  5.114893  7.324785 1.002270  1400
lambda[15]      4.588128   1.529048  1.986595  3.537436  4.479550  5.476045  7.997842 1.000895 12000
lambda[16]      7.035670   2.235785  3.776619  5.428774  6.617842  8.235651 12.389102 1.001082  9600
lambda[17]      7.586794   2.444835  4.078003  5.780102  7.145878  8.986080 13.364148 1.000973 12000
lambda[18]      3.413945   1.459583  0.955955  2.306502  3.321237  4.407808  6.427175 1.002377  1300
lambda[19]      6.559620   2.077231  3.418885  5.098825  6.196538  7.664223 11.574183 1.001324  4400
lambda[20]      3.391437   1.459840  0.949251  2.262276  3.281995  4.407545  6.371617 1.001482  3300
maybe_overdisp  1.119677   0.178304  0.900895  1.006214  1.073253  1.185895  1.565368 1.002841  1000
sigma           0.458185   0.235224  0.048157  0.294464  0.441186  0.597933  0.972298 1.009139   960
deviance       74.818941   5.985537 64.100886 70.435067 74.417487 79.115472 86.443865 1.001312  4600

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 17.9 and DIC = 92.7
DIC is an estimate of expected predictive error (lower deviance is better).
$\endgroup$
18
$\begingroup$

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))

Posterior density for the dispersion index 2, NB model

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))

Posterior density for the dispersion index 2

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))

Posterior density of the DI for "arable"

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:

Dispersion index quasi-Poisson

$\endgroup$
12
  • 1
    $\begingroup$ Thank you for your very thorough answer!! But I am sorry I cannot make a big use of it - it is like complete manual for professional statisticians, beyond my capabilities. I expected something simple - like you will have a look at my bugs model and show me how to compute the dispersion parameter from my parameters shown in print(out). Or that you will say: your model has to be extended by adding this and this to be able to find out the dispersion parameter. Sorry if your answer already contains this, I am not able to get the information out of it - too much information and new terms to me... $\endgroup$ – Tomas Jul 23 '13 at 11:42
  • 1
    $\begingroup$ I have updated my answer by adding output both from glm() and bugs. $\endgroup$ – Tomas Jul 23 '13 at 11:42
  • $\begingroup$ Try to add the following code to your WinBUGS code: kappa <- exp(alpha + beta); mean.x <- exp(alpha + 0.5*sigma2); vari <- kappa*exp(sigma2/2)+kappa*kappa*exp(2*sigma2) - kappa*kappa*exp(sigma2); DI <- vari/mean.x Then set a node for DI and paste the output into your question again. $\endgroup$ – COOLSerdash Jul 23 '13 at 11:51
  • 2
    $\begingroup$ @Tomas: Please read my comments. The WinBUGS model is a Poisson-Log-Normal model and not a quasi-Poisson model. The dispersion parameters cannot be compared, they will never be the same. So don't try to get 1.47 using WinBUGS, you won't. The overdispersion parameter in the quasi-Poisson glm is calculated as follows: 1/(n-k)*sum((C.OD-fitted(glm.fit.no.OD))^2/fitted(glm.fit.no.OD)) where n is the sample size and k the number of parameters (i.e. 2). $\endgroup$ – COOLSerdash Jul 23 '13 at 12:44
  • 1
    $\begingroup$ @Tomas The Poisson log-normal method (or negative binomial etc.) is perfectly fine for modelling overdispersion, why do you want a quasi-Poisson model in WinBUGS? Quasi-Likelihood is a method to model the relationship between the mean and the variance. Of course the model has a mathematical formulation (just google it) but I doubt that this will help you in any way. My suggestion is to use the perfectly fine Bayesian methodology to fit your model. $\endgroup$ – COOLSerdash Jul 23 '13 at 13:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.