Bayesian Poisson Regression

Bayesian Poisson Regression

The Bayesian Poisson regression model (also called log-linear model) is
$y_i\mid \mathbf{\beta },{\mathbf{x}}_i\sim \text{Poisson}(\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })),\mathbf{\beta }\sim N({\mathbf{\beta }}_0,{\mathbf{\Sigma }}_0).$
The log-link forces ${\lambda }_i=\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })>0$. The likelihood is
$p(\mathbf{y}\mid \mathbf{\beta })={\prod }_{i=1}^n\frac{e^{-{\lambda }_i}{\lambda }_i^{y_i}}{y_i!},\log p(\mathbf{y}\mid \mathbf{\beta })={\sum }_i[y_i{\mathbf{x}}_i^{\top }\mathbf{\beta }-\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })]+\text{const}.$
Like logistic regression, the Poisson likelihood is not conjugate to any standard prior on $\mathbf{\beta }$; inference uses the Laplace approximation or MCMC.

MAP via Newton (IRLS analog). The log-posterior is concave in $\mathbf{\beta }$. The gradient and Hessian are
$\nabla \log p(\mathbf{\beta }\mid \mathbf{y})={\mathbf{X}}^{\top }(\mathbf{y}-\mathbf{\lambda }(\mathbf{\beta }))-{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0),$
$-{\nabla }^2\log p(\mathbf{\beta }\mid \mathbf{y})={\mathbf{X}}^{\top }\mathrm{diag}(\mathbf{\lambda }(\mathbf{\beta }))\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}.$
Newton’s method converges quadratically, giving the MAP $\hat{\mathbf{\beta }}$ and the Laplace-posterior covariance ${\mathbf{\Sigma }}_{\text{Lap}}=({\mathbf{X}}^{\top }\mathrm{diag}(\hat{\mathbf{\lambda }})\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}$.

Laplace approximation: $p(\mathbf{\beta }\mid \mathbf{y})\approx N(\hat{\mathbf{\beta }},{\mathbf{\Sigma }}_{\text{Lap}})$. This Gaussian serves equally well as an approximate posterior summary or as a proposal distribution for Metropolis-Hastings when Laplace’s accuracy is insufficient — see 3.11.1 Metropolis-Hastings.

Posterior predictive for a new input ${\mathbf{x}}_{\star }$: conditional on $\mathbf{\beta }$, ${\tilde{y}}_{\star }\sim \text{Poisson}(\exp ({\mathbf{x}}_{\star }^{\top }\mathbf{\beta }))$; marginally ${\tilde{y}}_{\star }\mid \mathbf{y}$ is a Poisson-lognormal mixture. A useful approximation uses a Gamma-Poisson (Gamma-Poisson) matching: moment-match $\exp ({\mathbf{x}}_{\star }^{\top }\mathbf{\beta })$ to a Gamma distribution, yielding a Negative-Binomial predictive with explicit overdispersion.

Overdispersion caveat. Real count data typically have $\mathrm{Var}(y)>\mathbb{E}(y)$, violating the Poisson equidispersion assumption. Bayesian remedies: (a) add an observation-level random effect ${\eta }_i\sim N(0,{\tau }^2)$ so that $y_i\sim \text{Poisson}(\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta }+{\eta }_i))$ — the log-normal Poisson model; (b) replace Poisson with Negative Binomial (arising as a Gamma-Poisson mixture over ${\lambda }_i$) at the likelihood level. Both remove the rigid mean-variance coupling.

Prerequisites

• 3.1 Bayesian Framework — Bayesian mechanics.
• 3.8 Bayesian Regression — parent.
• 3.9 Laplace Approximation — standard non-conjugate tool.
• 3.5.3 Gamma-Poisson Model — univariate conjugate baseline.

Intuition

Poisson regression models the log of the rate as a linear function of covariates: each unit increase in $x_j$ multiplies the rate by $e^{{\beta }_j}$. The log link ensures positivity without constraint, at the cost of breaking conjugacy with Gaussian priors on $\mathbf{\beta }$. The log-posterior is the sum of a Poisson log-likelihood (strictly concave in $\mathbf{\beta }$) and a Gaussian log-prior (quadratic), so the posterior is globally log-concave with a unique maximum — ideal for Laplace approximation.

The concavity makes Newton’s method converge quickly from any reasonable starting point. Each Newton step reweights observations by ${\lambda }_i=e^{{\mathbf{x}}_i^{\top }{\mathbf{\beta }}_k}$: observations with high predicted rate carry proportionally more weight, mirroring the Fisher information of a Poisson variate (whose variance equals its mean). The Bayesian addition is the regularization term ${\mathbf{\Sigma }}_0^{-1}$, which guarantees invertibility of the Newton step even when the data are sparse or the design matrix is rank-deficient.

The overdispersion problem is the elephant in the Poisson room. Empirical count data almost always have variance-to-mean ratios well above 1, and fitting Poisson regression to them produces posterior intervals that are too narrow. The cleanest fix is a log-normal random effect per observation — which, integrated out, gives a log-normal-Poisson mixture with adjustable overdispersion. Bayesian software handles this naturally because the random effect fits into the hierarchical-model framework (3.6 Hierarchical Models) and the Laplace approximation extends straightforwardly.

Properties

(Log-Link Forces Positive Rates) ${\lambda }_i=\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })>0$ for every $\mathbf{\beta }$.

This sidesteps any positivity constraint on the linear predictor but means that coefficients are interpreted multiplicatively: $e^{{\beta }_j}$ is the rate ratio for a unit increase in $x_j$. A posterior credible interval for ${\beta }_j$ translates to a credible interval for the rate ratio via exponentiation.

(Log-Posterior is Strictly Concave) The MAP is unique and globally optimal.

The Hessian $-{\nabla }^2\log p(\mathbf{\beta }\mid \mathbf{y})={\mathbf{X}}^{\top }\mathrm{diag}(\mathbf{\lambda })\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}$ is positive definite (positive semi-definite first term, positive definite second term), so the log-posterior is strictly concave. Newton’s method converges quadratically to the MAP from any reasonable starting point. The prior-precision term is what rescues the non-informative case: the raw Poisson likelihood alone can have a non-strictly-concave posterior if the design has linear dependence.

(Laplace as MH Proposal) The Laplace Gaussian is a standard Metropolis-Hastings proposal distribution.

When full MH is required (e.g., with non-Gaussian priors or skewed posteriors), the proposal $N(\hat{\mathbf{\beta }},c^2{\mathbf{\Sigma }}_{\text{Lap}})$ with $c^2=2.38^2/d$ (Roberts-Gelman-Gilks optimal scaling for random-walk MH targeting a Gaussian) gives high acceptance rates in practice. This combines Laplace’s locality with MH’s exactness: if the true posterior is nearly Gaussian, the proposal is near-perfect; if it deviates, MH corrects. See 3.11.1 Metropolis-Hastings.

(Offset Terms Handle Exposure) Known exposures $t_i$ enter as offsets in the linear predictor.

If observation $i$ is a count over exposure $t_i$ (time, area, population at risk), the model becomes $y_i\sim \text{Poisson}(t_i\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta }))$. Equivalently, a fixed term $\log t_i$ is added to the linear predictor:
$\log {\lambda }_i=\log t_i+{\mathbf{x}}_i^{\top }\mathbf{\beta }.$
This offset carries no coefficient — its effect is fixed to 1. Bayesian and frequentist software support this as standard; the posterior structure on $\mathbf{\beta }$ is unchanged.

(Overdispersion Diagnostic) A log-normal-Poisson or Negative-Binomial extension fixes equidispersion failures.

If Pearson residuals from the fitted Poisson model have variance $\gg 1$, the data are overdispersed. Two Bayesian extensions:

• Log-normal Poisson: ${\eta }_i\sim N(0,{\tau }^2)$ with $y_i\sim \text{Poisson}(\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta }+{\eta }_i))$. Integrating out ${\eta }_i$ gives a log-normal-Poisson marginal with unbounded dispersion; fit by Laplace or MCMC.
• Negative-Binomial: ${\lambda }_i\sim \text{Gamma}(r,r/{\mu }_i)$ with $y_i\sim \text{Poisson}({\lambda }_i)$, integrating to $y_i\sim \text{NegBin}(r,{\mu }_i)$. Fit by MCMC; $r$ is a learned overdispersion parameter.
  The log-normal version fits naturally into the hierarchical-models framework and is the default in Bayesian GLMM software.

(Evidence for Model Comparison) Laplace approximation gives closed-form log-marginal-likelihood.

$\log p(\mathbf{y})\approx \log p(\mathbf{y}\mid \hat{\mathbf{\beta }})+\log p(\hat{\mathbf{\beta }})+\frac{d}{2}\log (2\pi )-\frac{1}{2}\log |{\mathbf{X}}^{\top }\mathrm{diag}(\hat{\mathbf{\lambda }})\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}|.$
This enables Bayes-factor-based model selection and evidence-tuned prior precisions; see 3.7 Marginal Likelihood. The asymptotic BIC surrogate drops the prior-density term and uses $\log n$ in place of the Hessian log-determinant.

(Sparse Counts Require Informative Priors) When most $y_i=0$, flat priors lead to unstable estimates.

In sparse count data (e.g., rare disease incidence), the MLE can be extreme or undefined for covariate levels that appear only with $y_i=0$. A weakly informative prior (e.g., $\mathbf{\beta }\sim N(0,2.5^2\mathbf{I})$ on standardized covariates) stabilizes the MAP and preserves sensible credible intervals. This is the count-data analog of logistic regression’s separability problem.

Examples

Example 1: Single-Covariate Poisson Regression

Data $\{(x_i,y_i){\}}_{i=1}^{50}$ with $x_i\sim N(0,1)$ and $y_i\sim \text{Poisson}(\exp (1+0.3x_i))$ (baseline rate $e\approx 2.72$, rate ratio $e^{0.3}\approx 1.35$). Prior $\mathbf{\beta }=({\beta }_0,{\beta }_1)\sim N(0,10\mathbf{I})$.

Newton iteration: converges in 4 iterations to $\hat{\mathbf{\beta }}\approx (1.02,0.31)$.
Laplace posterior standard deviations: from ${\mathbf{\Sigma }}_{\text{Lap}}$ diagonals, ${\sigma }_{{\beta }_0}\approx 0.09$, ${\sigma }_{{\beta }_1}\approx 0.07$.
95% CI for rate ratio $e^{{\beta }_1}$: $e^{0.31\pm 1.96\cdot 0.07}=[1.18,1.55]$, bracketing the true 1.35.
Posterior predictive at $x_{\star }=1$: posterior mean rate ${\hat{\lambda }}_{\star }=\exp (1.02+0.31)=3.78$; predictive mean matches, but predictive distribution is Poisson-lognormal with inflated variance relative to $\text{Poisson}(3.78)$.

Example 2: Exposure-Adjusted Disease Incidence

A public-health study tallies disease cases in $J=50$ counties over one year, with county populations $t_j$ (offset) and a single risk covariate $x_j$ (e.g., average pollutant concentration).

Model: $y_j\sim \text{Poisson}(t_j\exp ({\beta }_0+{\beta }_1{x}_j))$. The offset $\log t_j$ is fixed; coefficient ${\beta }_1$ measures the log-rate increase per unit pollution.
Prior: ${\beta }_0\sim N(-10,5)$ (placing the baseline rate near $e^{-10}\approx 5\times 10^{-5}$ per person-year), ${\beta }_1\sim N(0,1)$ (weakly informative).
Posterior: Laplace approximation with Newton’s method.
Interpretation: if posterior mean ${\hat{\beta }}_1=0.15$ with SD 0.05, the rate ratio for one unit of pollution is $e^{0.15}=1.16$ with 95% CI $[1.05,1.28]$, a 16% per-unit increase significantly different from 1.
Model check: compare fitted variances to observed variance across counties. If Pearson residuals show large dispersion, add a county-level random effect ${\eta }_j\sim N(0,{\tau }^2)$ to upgrade to a log-normal Poisson hierarchical model.

Example 3: Laplace Proposal for Metropolis-Hastings

Suppose the analyst prefers full MCMC over Laplace point estimates (perhaps for better tail summaries or for a non-Gaussian prior). Use the Laplace-approximate posterior as the MH proposal:

1. Preprocess: Newton-iterate to MAP $\hat{\mathbf{\beta }}$; compute Laplace covariance ${\mathbf{\Sigma }}_{\text{Lap}}$.
2. MH step $t$: propose ${\mathbf{\beta }}^{\ast }\sim N({\mathbf{\beta }}^{(t-1)},c^2{\mathbf{\Sigma }}_{\text{Lap}})$ with $c=2.38/\sqrt{d}$ (Roberts-Gelman-Gilks optimal scaling).
3. Accept with probability $\alpha =\min (1,\exp [\log p({\mathbf{\beta }}^{\ast })-\log p({\mathbf{\beta }}^{(t-1)})])$ — proposal is symmetric, so no proposal-density ratio is needed.

For the Poisson regression in Example 1, this typically achieves acceptance rates near the 23% optimum for $d=2$ and produces effective sample sizes comparable to conjugate Gibbs samplers. See 3.11.1 Metropolis-Hastings.

Derivation

Goal: derive the Laplace approximation and IRLS-analog Newton updates for Poisson regression.

The log-posterior (up to $\mathbf{\beta }$-free constants):
$\mathcal{L}(\mathbf{\beta })={\sum }_{i=1}^n[y_i{\mathbf{x}}_i^{\top }\mathbf{\beta }-\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })]-\frac{1}{2}(\mathbf{\beta }-{\mathbf{\beta }}_0{)}^{\top }{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0).$
Differentiating once:
$\nabla \mathcal{L}={\sum }_i{\mathbf{x}}_i(y_i-{\lambda }_i)-{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0)={\mathbf{X}}^{\top }(\mathbf{y}-\mathbf{\lambda })-{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0).$
Differentiating twice:
${\nabla }^2\mathcal{L}=-{\sum }_i{\mathbf{x}}_i{\mathbf{x}}_i^{\top }{\lambda }_i-{\mathbf{\Sigma }}_0^{-1}=-{\mathbf{X}}^{\top }\mathrm{diag}(\mathbf{\lambda })\mathbf{X}-{\mathbf{\Sigma }}_0^{-1},$
which is negative definite. The MAP is the unique solution to $\nabla \mathcal{L}=0$, found by Newton’s method:
${\mathbf{\beta }}_{k+1}={\mathbf{\beta }}_k-[{\nabla }^2\mathcal{L}{]}^{-1}\nabla \mathcal{L}={\mathbf{\beta }}_k+({\mathbf{X}}^{\top }{\mathbf{W}}_k\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}[{\mathbf{X}}^{\top }(\mathbf{y}-{\mathbf{\lambda }}_k)-{\mathbf{\Sigma }}_0^{-1}({\mathbf{\beta }}_k-{\mathbf{\beta }}_0)],$
with ${\mathbf{W}}_k=\mathrm{diag}({\mathbf{\lambda }}_k)$. Defining the working response $z_{k,i}={\mathbf{x}}_i^{\top }{\mathbf{\beta }}_k+(y_i-{\lambda }_{k,i})/{\lambda }_{k,i}$ casts this as a weighted least-squares fit, the Poisson analog of IRLS.

At convergence, ${\hat{\mathbf{\lambda }}}_i=\exp ({\mathbf{x}}_i^{\top }\hat{\mathbf{\beta }})$ and the Laplace covariance is
${\mathbf{\Sigma }}_{\text{Lap}}=({\mathbf{X}}^{\top }\mathrm{diag}(\hat{\mathbf{\lambda }})\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}.$

Log-marginal-likelihood: substituting into the standard Laplace formula (3.9 Laplace Approximation),
$\log p(\mathbf{y})\approx \log p(\mathbf{y}\mid \hat{\mathbf{\beta }})+\log p(\hat{\mathbf{\beta }})+\frac{d}{2}\log (2\pi )+\frac{1}{2}\log |{\mathbf{\Sigma }}_{\text{Lap}}|.$

See Also

• 3.8 Bayesian Regression — parent index.
• 3.9 Laplace Approximation — underlying approximation technique.
• 3.5.3 Gamma-Poisson Model — univariate conjugate baseline; also the basis of Negative-Binomial overdispersion extensions.
• 3.8.3 Bayesian Logistic Regression — non-conjugate binary-response counterpart with the same Laplace recipe.
• 3.11.1 Metropolis-Hastings — the Laplace Gaussian is a standard MH proposal in this setting.
• 3.11.2 Gibbs Sampling — hierarchical extensions fit into Gibbs via the log-normal random-effect formulation.
• 3.6 Hierarchical Models — overdispersion via observation-level random effects.
• 3.7 Marginal Likelihood — Laplace-approximate evidence for model comparison.