Chapter 10: Models for the Poisson family

In glmbayes: Bayesian Generalized Linear Models (IID Samples)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(glmbayes)

1. Introductory Discussion

Poisson regression is the canonical GLM for modeling count data [@McCullagh1989].
The Poisson likelihood is log‑concave under the log link, making it an ideal setting for envelope‑based Bayesian sampling [@Nygren2006].
This chapter mirrors the structure used for the Binomial vignette: we begin with the classical glm() fit, then introduce the Bayesian glmb() version, discuss prior specification, and compare posterior summaries to their classical counterparts.

We use the well‑known randomized controlled trial dataset of [@Dobson1990], which appears in many GLM textbooks and is included in the base R documentation for glm().
The dataset contains counts of subjects classified by treatment group and outcome category.

2. Poisson Likelihood and Model Structure

Poisson regression models count data of the form

[ Y_i \in {0,1,2,\ldots}, \qquad Y_i \sim \text{Poisson}(\mu_i), ]

where (\mu_i > 0) is the expected count.
Observation weights (w_i) (e.g., exposure times) enter the likelihood through

[ Y_i \sim \text{Poisson}(w_i \mu_i). ]

The mean is linked to a linear predictor via

[ \eta_i = x_i^\top \beta, \qquad \mu_i = g^{-1}(\eta_i). ]

2.1 Weighted Poisson Log‑Likelihood

Up to constants, the weighted log‑likelihood is

[ \ell(\beta) = \sum_{i=1}^n \left[ w_i y_i \log(\mu_i) - w_i \mu_i \right]. ]

Under the log link, the mean is

[ \mu_i = e^{\eta_i}, ]

and the log‑likelihood simplifies to

[ \ell(\beta) = \sum_{i=1}^n \left[ w_i y_i \eta_i - w_i e^{\eta_i} \right]. ]

This is the form used by both glm() and the Bayesian functions glmb() and rglmb().

2.2 Exponential‑Family Representation

The Poisson distribution belongs to the exponential family with canonical parameter

[ \theta_i = \log(\mu_i), ]

and cumulant function

[ b(\theta_i) = e^{\theta_i}. ]

The variance is

[ \mathrm{Var}(Y_i) = \mu_i. ]

Under the log link, the linear predictor equals the canonical parameter:

[ \eta_i = \theta_i. ]

This identity ensures that the Poisson log‑likelihood is globally log‑concave in the linear predictor, which is ideal for envelope construction and accept–reject sampling in glmbayes.

2.3 Log Link and Its Properties

The canonical link for the Poisson family is the log link:

[ g(\mu_i) = \log(\mu_i), \qquad \mu_i = e^{\eta_i}. ]

Key properties:

• automatically enforces (\mu_i > 0)
• yields a concave log‑likelihood in (\eta_i)
• ensures stable IRLS updates in classical GLMs
• supports efficient envelope‑based iid sampling in glmbayes

Noncanonical links (identity, square‑root) do not preserve log‑concavity and are not supported in glmb().

2.4 Likelihood, Prior, and Posterior (Normal Prior on beta)

With the log link, the weighted log‑likelihood is

[ \ell(\beta) = \sum_{i=1}^n \left[ w_i y_i \eta_i - w_i e^{\eta_i} \right]. ]

A multivariate normal prior is specified as

[ \log p(\beta) = -\tfrac{1}{2}(\beta - \mu_0)^\top \Sigma_0^{-1} (\beta - \mu_0) + \text{const}. ]

Combining likelihood and prior yields the log‑posterior:

[ \log p(\beta \mid y) = \sum_{i=1}^n \left[ w_i y_i \eta_i - w_i e^{\eta_i} \right] - \tfrac{1}{2}(\beta - \mu_0)^\top \Sigma_0^{-1} (\beta - \mu_0) + \text{const}. ]

Because both terms are concave in (\beta), the posterior is log‑concave, enabling efficient iid sampling via the envelope‑based accept–reject sampler implemented in glmbayes.

3. Data Setup

We reproduce Dobson’s example exactly.
The data consist of 9 observations: three outcome levels crossed with three treatment groups.

    ## Dobson (1990) Page 93: Randomized Controlled Trial :
    set.seed(333)
    counts <- c(18,17,15,20,10,20,25,13,12)
    outcome <- gl(3,1,9)
    treatment <- gl(3,3)
    print(d.AD <- data.frame(treatment, outcome, counts))

The table shows the observed counts for each treatment–outcome combination.

4. Classical Poisson Regression

The classical glm() call uses the Poisson family with the canonical log link:

    glm.D93 <- glm(counts ~ outcome + treatment,
                   family = poisson(link = "log"))
    summary(glm.D93)

4.1 Interpretation of the Classical Output

The glm summary reports:

• Coefficient estimates on the log‑scale
• Standard errors derived from the Fisher information
• Residual deviance and null deviance
• AIC, which penalizes model complexity

The coefficients represent multiplicative effects on the expected count.
For example, (\exp(\beta)) gives the rate ratio associated with a one-unit change in a predictor.

The residual deviance (approximately 5.13) indicates a good fit relative to the null deviance (approximately 10.58), consistent with Dobson's original analysis.

5. Bayesian Poisson Regression with glmb()

To fit the Bayesian analogue, we must specify a prior distribution for the regression coefficients.
The recommended workflow uses Prior_Setup(), which constructs a Zellner‑type g‑prior aligned with the model matrix.

5.1 Setting the Prior

    ps <- Prior_Setup(counts ~ outcome + treatment,
                      family = poisson())
    mu <- ps$mu
    V  <- ps$Sigma

The printed Prior_Setup() output shows:

• prior means (intercept from null model, effects centered at 0)
• prior standard deviations (scaled by ( g = (1 - pwt) / pwt ))
• the implied pseudo-sample size ( n_{\text{prior}} )

With the default pwt = 0.01, the prior is weakly informative.

5.2 Calling glmb()

    glmb.D93 <- glmb(counts ~ outcome + treatment,
                     family = poisson(),
                     pfamily = dNormal(mu = mu, Sigma = V))

This call mirrors glm() but adds the pfamily argument.

6. Printing and Summarizing the Bayesian Output

    print(glmb.D93)
    summary(glmb.D93)

The summary includes:

• Posterior means and posterior modes
• Posterior standard deviations
• Monte Carlo error estimates
• Tail probabilities relative to the prior mean
• Effective number of parameters (pD)
• Expected residual deviance ((\bar{D}))
• DIC, the Bayesian analogue of AIC

6.1 Posterior Mean Coefficients

Posterior means are close to the classical MLEs because:

• the prior is weak (pwt = 0.01)
• the Poisson likelihood is strongly log‑concave
• the sample size is modest but informative

6.2 Effective Number of Parameters

The pD value is typically slightly below the nominal number of coefficients, reflecting mild shrinkage from the prior.

6.3 DIC

DIC is approximately 56.4, very close to the classical AIC value of 56.76.
This agreement is expected when the prior is weak and the posterior distribution is close to Gaussian.

7. Comparison of Classical and Bayesian Estimates

| Quantity | Classical glm | Bayesian glmb | |----------------|-------------------|------------------------| | Point estimate | MLE | Posterior mean | | Uncertainty | SE | Posterior SD | | Model fit | Deviance, AIC | (\bar{D}), (p_D), DIC | | Interpretation | log-rate ratios | same, with shrinkage |

The Bayesian model provides richer diagnostics (tail probabilities, posterior quantiles) while preserving the familiar GLM structure.

8. Credible Intervals and Posterior Interpretation

Posterior percentiles (2.5%, 50%, 97.5%) provide Bayesian credible intervals for each coefficient.
These intervals closely track the classical Wald intervals because the posterior is nearly normal.

9. Concluding Remarks

This vignette demonstrates that:

• Poisson GLMs with log links are ideal for envelope‑based Bayesian sampling
• glmb() mirrors glm() while providing full posterior inference
• Prior_Setup() offers a principled default prior
• Posterior summaries align closely with classical results under weak priors

In later chapters, we extend these ideas to:

• overdispersed quasi‑Poisson models
• informative priors
• prediction and posterior predictive checks
• GPU‑accelerated Poisson models for large datasets

---

Appendix A. Bayes Rules! companion — equality_index Poisson regression

Book: [@JohnsonLauEtAl2022], Chapter 12 — Poisson regression
Data: bayesrules::equality_index (California removed, as in the book)
Model: laws ~ percent_urban + historical (log link)
Priors in this appendix: Prior_Setup() defaults only. The book uses an informative intercept prior plus weak / autoscaling priors on slopes (partially informative specification); we do not mix in the book intercept prior here.

Requires bayesrules. The table below compares book posterior summaries (tidy(equality_model), Ch. 12) to glmb() under glmbayes default priors.

library(bayesrules)
equality <- bayesrules::equality_index
equality <- equality[equality$laws < max(equality$laws), ]

ps_eq <- Prior_Setup(laws ~ percent_urban + historical,
                     family = poisson(), data = equality)

## Bayes Rules! Ch. 12 tidy(equality_model) posterior estimates (log link)
book_br10 <- data.frame(
  parameter = c("(Intercept)", "percent_urban", "historicalgop", "historicalswing"),
  book_mean = c(1.71, 0.0164, -1.52, -0.610),
  book_sd   = c(0.303, 0.00353, 0.134, 0.103),
  check.names = FALSE
)

set.seed(2026)
glmb_eq <- glmb(
  laws ~ percent_urban + historical,
  family  = poisson(),
  pfamily = dNormal(mu = ps_eq$mu, Sigma = ps_eq$Sigma),
  data    = equality,
  n       = 2000
)
print(glmb_eq)

br10_compare <- data.frame(
  parameter = book_br10$parameter,
  `Book mean` = book_br10$book_mean,
  `Book SD`   = book_br10$book_sd,
  `glmb Post.Mean` = as.numeric(glmb_eq$coef.means[book_br10$parameter]),
  `glmb Post.Sd`   = sapply(book_br10$parameter, function(p)
    sd(glmb_eq$coefficients[, p, drop = TRUE])),
  check.names = FALSE
)
knitr::kable(br10_compare, digits = 4,
  caption = "Bayes Rules! Ch. 12 (informative + weak book priors) vs. glmb() with Prior_Setup() defaults")

Posterior means will differ from the book because the book's fit combines an informative intercept with weakly informative slope priors, while this fit uses a single weak Prior_Setup() specification throughout.

glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.