Poisson Distribution for Counts and Rates

In plain language

The Poisson distribution is a mathematical formula for "how many times will this rare event happen?" when events occur independently of one another — for example, how many new heart-failure hospitalizations will occur across a group of patients over a year of follow-up. Instead of a head count of patients who had the event, it divides the number of events by the total time all patients were actively being watched (summing every patient's individual follow-up length), producing a rate like "20 hospitalizations per 1000 patient-years of observation." That rate is the distribution's single number λ (lambda), which also happens to equal its spread — a convenient property called equidispersion that almost always breaks down in real healthcare data where a small group of high-utilizers inflates the variance, signalling the analyst to upgrade to the negative binomial model instead.

The Poisson distribution, parameterized by a single rate λ (lambda), assigns to each non-negative integer k the probability P(X = k) = e^(−λ)·λ^k / k!. When the observation unit is one patient-year of follow-up, λ is the incidence rate — the most fundamental summary in RWE. Understanding the Poisson distribution as a probability model in its own right, not only as a regression family, clarifies why it fits, why it fails, and which fixes apply to which violations: negative-binomial overdispersion correction, robust sandwich standard errors, or exact confidence intervals for rare counts.

The Poisson process and its four assumptions

The Poisson distribution is the stationary distribution of a Poisson process, a continuous-time counting process characterized by four assumptions. First, independent increments: events in non-overlapping time intervals are statistically independent, so one burst of events does not predict the next. Second, stationary increments: the event rate λ is constant over time within the observation period — the probability of an event in any small window dt is λ·dt regardless of calendar time or how long the process has been running. Third, orderliness: at most one event can occur in an infinitesimally small window; simultaneous events have probability zero. Fourth, non-negative integer support: counts are 0, 1, 2, ... with no upper bound.

In claims and EHR data, all four assumptions are routinely violated. Independent increments fails when a single severe illness triggers a cascade of correlated events — hospitalization begets readmission begets ED visit, a pattern called event clustering that is pervasive in complex-care populations. Stationarity fails because seasonal patterns (winter respiratory illness spikes, COVID-driven disruption of care-seeking in 2020) make λ vary with calendar time. Orderliness fails when same-day claims from multiple providers are grouped under a single service date. Despite these violations, the Poisson model with correctly specified robust standard errors (sandwich variance) remains a useful workhorse for point estimation of rates and rate ratios; it is the standard errors — not the point estimates — that suffer most when the assumptions fail.

Equidispersion: the defining constraint and when it fails

The Poisson distribution has one free parameter: λ governs both the mean and the variance, so E(X) = Var(X) = λ. This identity — equidispersion — is not an assumption layered on top; it is the mathematical consequence of the Poisson process and is therefore always testable from the data. A Pearson dispersion statistic (Pearson chi-squared divided by residual degrees of freedom) near 1.0 is consistent with equidispersion; values substantially greater than 1 signal overdispersion (variance exceeds mean), and values substantially less than 1 signal underdispersion (rare in healthcare data). In real-world healthcare utilization data, overdispersion is almost universal: a minority of high-utilizers — frail elderly patients, individuals with multi-morbidity — generates a long right tail of event counts whose variance far exceeds the mean. Under overdispersion, Poisson regression produces point estimates of the log-rate-ratio that are roughly unbiased, but the standard errors are anticonservatively small, confidence intervals too narrow, and p-values inflated. The canonical fix is negative binomial regression, which adds a dispersion parameter α so that Var(X) = λ + αλ² (NB-2 parameterization); a lighter-weight alternative is Poisson regression with robust sandwich standard errors, which corrects the SEs without committing to the NB likelihood.

Person-time denominators and the offset trick

The central pattern in RWE incidence-rate analysis is Poisson regression with log(person-time) entered as an offset — a covariate whose coefficient is constrained to exactly 1 — rather than as a free predictor. The model becomes log(μ_i) = β₀ + β₁X_i + log(T_i), equivalent to log(μ_i / T_i) = β₀ + β₁X_i, where μ_i/T_i is the event rate for patient i per unit person-time. The exponentiated coefficient exp(β₁) is then an incidence rate ratio (IRR) comparing event rates per unit person-time between exposure groups, conditional on covariates. Omitting the offset — or treating log(person-time) as a regular covariate with a free coefficient — corrupts the estimand: without the offset constraint, the model predicts raw counts rather than rates, and patients followed longer appear to have more events simply because they were observed longer, producing a bias that mirrors immortal-time distortion.

In claims data, person-time is the sum of each patient's continuously observable FFS-enrolled days, censored at disenrollment, death, or the administrative study end. Enrollment gaps, Medicare Advantage-only spans (where inpatient claims are absent and event counts are incomplete), and pre-index washout periods must be excluded from the observable person-time before computing the offset. Constructing this denominator correctly is upstream infrastructure; its quality determines the validity of every downstream Poisson or negative binomial rate model.

Exact Poisson confidence intervals for rare counts

When event counts are small — typically fewer than about 30 — the normal approximation for an incidence rate confidence interval (rate ± 1.96 × SE) is unreliable because the Poisson distribution is right-skewed at small λ. The exact Poisson CI (Garwood interval, or equivalently the gamma/chi-squared inversion) constructs the CI directly from the count D and total person-time T using chi-squared quantiles: lower = χ²(α/2, 2D) / (2T) and upper = χ²(1−α/2, 2D+2) / (2T). This exact CI is the standard for rare safety-signal reporting in pharmacovigilance — spontaneous reporting systems, observed-versus-expected analyses, and sequential signal detection methods such as maxSPRT and TreeScan. In R, poisson.test() supplies the exact CI from aggregate inputs; in Python, scipy.stats.chi2 provides the quantiles; in SAS, PROC STDRATE uses aggregate person-time and event counts to produce exact CIs directly.

The Poisson approximation to the binomial for rare events

When a binary event is rare (event probability p much less than 0.05) and the number of trial opportunities n is large, the binomial distribution B(n, p) converges to a Poisson distribution with λ = np. This approximation underpins the use of Poisson models in cohort studies of rare outcomes: a binary indicator — first stroke yes or no — with a low event rate across large accrued person-time behaves like a Poisson count, and Poisson regression serves as a log-linear approximation to the binomial with excellent accuracy. The approximation is reliable when p is below about 0.05 and np is below about 10. At higher event rates the binomial with logistic regression is the correct model for a binary outcome.

Poisson regression with robust variance: Zou's method for risk ratios on binary outcomes

In clinical trials and cohort studies, the risk ratio — the ratio of proportions — is often more interpretable and policy-relevant than the odds ratio. Logistic regression gives an odds ratio that approximates the risk ratio only when the event is rare; for common outcomes (prevalence above 10 to 20 percent), the odds ratio substantially overstates the risk ratio. A widely adopted workaround described by Zou (2004) is Poisson regression with robust sandwich standard errors applied to a binary outcome: treat the binary event indicator as if it were a count, fit a Poisson GLM with a log link and no offset (because all patients have the same fixed follow-up horizon), obtain the exponentiated coefficient as the risk ratio, and use the sandwich SE to correct for the misspecified variance structure (since a binary outcome is Bernoulli, not Poisson). The RR estimate from this modified-Poisson approach is consistent and the sandwich-corrected CIs have valid asymptotic coverage. This method is now ubiquitous in RWE publications and regulatory submissions for producing adjusted risk ratios from binary endpoints in large cohorts.

Interpreting the output

Consider a Poisson rate regression for a hospitalization outcome comparing an exposed cohort to an unexposed cohort, adjusted for age, sex, and a comorbidity index, with log(person-years) as the offset. The model returns a log-rate-ratio coefficient of 0.693 for the exposure indicator, with a model-based 95% confidence interval on the log scale of 0.41 to 0.98.

Formal interpretation. The exponentiated coefficient exp(0.693) equals approximately 2.0, which is the incidence rate ratio (IRR): the estimated hospitalization rate among the exposed cohort is 2.0 times the rate in the unexposed cohort per unit person-time, holding age, sex, and comorbidity constant. The 95% CI on the IRR scale corresponds to exp(0.41) to exp(0.98), approximately 1.51 to 2.66. This interval does NOT mean there is a 95% probability that the true IRR lies between 1.51 and 2.66; it means that in repeated sampling under the model's assumptions, 95% of such intervals would contain the true parameter. The equidispersion caveat is critical: if the Pearson dispersion statistic substantially exceeds 1.0, the model-based SEs are anticonservatively narrow and the CI is too tight — in that setting, report negative-binomial IRRs or Poisson IRRs with robust sandwich SEs rather than the model-based CI shown here.

Practical interpretation. Hospitalizations occur about twice as often per year of follow-up among the exposed group as among the unexposed group, after accounting for age, sex, and comorbidity burden. This is an association; whether it reflects a causal effect of the exposure depends entirely on the study design and on how well unmeasured confounding has been addressed.

Pros, cons, and trade-offs

Pros of the Poisson distribution as the baseline count model: a single interpretable parameter λ governs both the mean and variance; the log-link offset trick cleanly converts counts to rates; the IRR is directly interpretable and maps naturally to benefit-risk narratives in regulatory and HTA submissions; exact Poisson CIs require no approximation for rare safety events; the Poisson-binomial relationship provides a theoretical bridge to binary-outcome cohort analyses; Zou's robust-variance extension delivers risk ratios from binary outcomes without the rare-event restriction of logistic regression. Cons: equidispersion Var(X) = E(X) fails in virtually all real healthcare utilization data, requiring negative-binomial or quasi-Poisson variance correction; the constant-rate assumption is violated by seasonal clustering and secular trends; the independence assumption is violated by within-person event clustering. Versus negative binomial: the NB adds a dispersion parameter α that absorbs overdispersion at the cost of one additional parameter; reserve Poisson for genuinely equidispersed or rare-event counts and treat NB as the safe default for HCRU. Versus binomial with logistic link: Poisson targets rates with person-time denominators while binomial targets proportions over a fixed horizon; the two are connected through Zou's robust-variance approach when a risk ratio is desired for a common binary outcome.

When to use

Use the Poisson distribution or Poisson regression when computing crude or adjusted incidence rates from person-time denominators in a log-linear offset model; constructing exact Poisson confidence intervals for rare safety-signal counts in pharmacovigilance, observed-versus-expected analyses, or sequential signal detection; applying Zou's modified Poisson regression with robust SEs to estimate risk ratios for binary outcomes in large cohorts where the outcome is too common for the rare-event approximation; diagnosing the equidispersion baseline against which overdispersion is tested before choosing negative-binomial or quasi-Poisson variance correction; or modeling genuinely equidispersed count data such as sentinel event rates in tightly controlled procedural registers.

When NOT to use

Do not use naive Poisson regression with model-based SEs when the Pearson dispersion statistic substantially exceeds 1 and overdispersion is confirmed — use negative-binomial regression or quasi-Poisson SEs instead, because plain Poisson SEs will be anticonservatively narrow and will produce false positives at the conventional significance threshold. Do not model durations, time-to-first-event, or censored survival outcomes with Poisson regression — these belong to the Cox proportional-hazards model or parametric survival families that correctly handle the at-risk process and censoring structure. Do not use a Poisson model for recurrent events with strong within-patient correlation without adding cluster-robust (sandwich) variance or a frailty term; ignoring within-person clustering understates uncertainty, particularly when high-utilizer patients dominate the count distribution. Do not apply a single-rate Poisson model when the hazard is strongly time-varying — a peri-procedural spike followed by a flat tail, for example — without either splitting person-time into clinically meaningful intervals or including follow-up-time splines. Do not interpret a cause-specific Poisson rate as a patient-facing probability when competing events such as death are common in the cohort.

Worked example

Scenario

A pharmacoepidemiologist is measuring the incidence rate of a serious adverse event (SAE) in two cohorts followed from their first prescription. Cohort A (new drug) accumulated 24 SAE events over 1200 person-years of follow-up. Cohort B (standard care) accumulated 36 SAE events over 900 person-years. She wants the crude rate in each cohort per 1000 person-years and the rate ratio comparing the two.

Dataset

Aggregate event counts and person-time for two treatment cohorts. Each patient's individual follow-up days were summed to produce person_years; rate_per_py divides events by person_years.

Steps

• Cohort A crude rate: 24 / 1200 = 0.02 events per person-year.

• Rescale to per 1000 person-years: 0.02 * 1000 = 20 events per 1000 person-years.

• Cohort B crude rate: 36 / 900 = 0.04 events per person-year.

• Rescale: 0.04 * 1000 = 40 events per 1000 person-years.

• Rate ratio (B vs A) = 0.04 / 0.02 = 2.0. The SAE occurs twice as often per person-year in Cohort B as in Cohort A.

• For an exact Poisson 95% CI on Cohort A's rate: use chi-squared inversion on the raw count of 24 events over 1200 person-years. In R, poisson.test(24, T=1200) produces the exact interval; in Python, scipy.stats.chi2.ppf() supplies the same quantiles. The exact CI matters here because the normal approximation is unreliable for small counts.

Result

Rate A = 24 / 1200 = 0.02 per person-year; 0.02 1000 = 20 per 1000 person-years. Rate B = 36 / 900 = 0.04 per person-year; 0.04 1000 = 40 per 1000 person-years. Rate ratio = 0.04 / 0.02 = 2.0: the SAE occurs twice as often per person-year in Cohort B.

Runnable example

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import poisson, chi2
import statsmodels.api as sm
import statsmodels.formula.api as smf

# ── 1. Poisson PMF and CDF ──────────────────────────────────────────────────
lam = 2.5   # e.g., average 2.5 events per patient-year
k_vals = np.arange(0, 10)
pmf_vals = poisson.pmf(k_vals, mu=lam)
print("Poisson PMF (lambda=2.5):")
for k, p in zip(k_vals, pmf_vals):
    print(f"  P(X={k}) = {p:.4f}")
print(f"Mean = Variance = {lam} (equidispersion)")

# ── 2. Exact Poisson CI for an aggregate rate (Garwood / chi-squared inversion) ──
D  = 24          # observed events
T  = 1200.0      # person-years
alpha = 0.05

lo = chi2.ppf(alpha / 2,    2 * D    ) / (2 * T) if D > 0 else 0.0
hi = chi2.ppf(1 - alpha / 2, 2 * (D + 1)) / (2 * T)
print(f"\nExact Poisson 95% CI for {D} events / {T} PY:")
print(f"  Rate = {D/T:.4f} per PY = {D/T*1000:.1f} per 1000 PY")
print(f"  95% CI: ({lo*1000:.1f}, {hi*1000:.1f}) per 1000 PY")
# Equivalent using scipy.stats.poisson.interval (two-sided):
exact = stats.poisson.interval(0.95, D)          # count-space CI
lo2 = exact[0] / T; hi2 = exact[1] / T
print(f"  Via poisson.interval: ({lo2*1000:.1f}, {hi2*1000:.1f}) per 1000 PY")

# ── 3. Poisson GLM with person-time offset: incidence rate ratio ─────────────
# Synthetic analytic table; replace with real data.
np.random.seed(42)
n = 400
analytic = pd.DataFrame({
    "person_id":       np.arange(n),
    "exposed":         np.repeat([1, 0], n // 2),
    "event_count":     np.random.poisson(lam=[0.06 if e else 0.03 for e in np.repeat([1, 0], n // 2)]),
    "person_years":    np.random.uniform(0.5, 2.0, n),
    "age":             np.random.normal(60, 10, n),
    "sex":             np.random.binomial(1, 0.5, n),
    "baseline_comorb": np.random.poisson(2, n),
})
analytic = analytic[analytic["person_years"] > 0].copy()
offset_vec = np.log(analytic["person_years"].to_numpy())

model = smf.glm(
    "event_count ~ exposed + age + sex + baseline_comorb",
    data=analytic,
    family=sm.families.Poisson(),
    offset=offset_vec,
).fit(cov_type="cluster", cov_kwds={"groups": analytic["person_id"]})

print("\nPoisson rate model (cluster-robust SEs):")
irr    = np.exp(model.params["exposed"])
ci_lo  = np.exp(model.conf_int().loc["exposed", 0])
ci_hi  = np.exp(model.conf_int().loc["exposed", 1])
print(f"  IRR (exposed vs unexposed) = {irr:.3f}  95% CI: ({ci_lo:.3f}, {ci_hi:.3f})")
print(f"  Pearson dispersion = {model.pearson_chi2 / model.df_resid:.3f}  (should be ~1 for equidispersion)")

# ── 4. Zou's modified Poisson with robust SEs: risk ratio on a binary outcome ──
analytic_binary = pd.DataFrame({
    "person_id": np.arange(n),
    "exposed":   np.repeat([1, 0], n // 2),
    "outcome":   np.random.binomial(1, [0.25 if e else 0.15 for e in np.repeat([1, 0], n // 2)]),
    "age":       analytic["age"].values,
    "sex":       analytic["sex"].values,
})

zou = smf.glm(
    "outcome ~ exposed + age + sex",   # no offset: fixed follow-up window
    data=analytic_binary,
    family=sm.families.Poisson(),      # Poisson family, log link
).fit(cov_type="HC1")                  # robust (sandwich) SEs correct Bernoulli variance

rr    = np.exp(zou.params["exposed"])
rr_lo = np.exp(zou.conf_int().loc["exposed", 0])
rr_hi = np.exp(zou.conf_int().loc["exposed", 1])
print(f"\nZou modified Poisson (robust SEs) — risk ratio:")
print(f"  RR = {rr:.3f}  95% robust CI: ({rr_lo:.3f}, {rr_hi:.3f})")
print("  (Use this instead of logistic OR when outcome prevalence > 10%)")

library(sandwich)
library(lmtest)

# ── 1. Poisson PMF and CDF ──────────────────────────────────────────────────
lam <- 2.5
k   <- 0:9
cat("Poisson PMF (lambda = 2.5):\n")
cat(sprintf("  P(X=%d) = %.4f\n", k, dpois(k, lambda = lam)))
cat(sprintf("Mean = Variance = %.1f (equidispersion)\n\n", lam))

# ── 2. Exact Poisson CI via poisson.test (Garwood / mid-p exact) ─────────────
D <- 24L;   T_py <- 1200.0
pt <- poisson.test(D, T = T_py, conf.level = 0.95)
cat(sprintf("Exact Poisson 95%% CI for %d events / %.0f PY:\n", D, T_py))
cat(sprintf("  Rate = %.4f per PY = %.1f per 1000 PY\n", D / T_py, D / T_py * 1000))
cat(sprintf("  95%% CI: (%.1f, %.1f) per 1000 PY\n\n",
            pt$conf.int[1] * 1000, pt$conf.int[2] * 1000))

# ── 3. Poisson GLM with person-time offset: incidence rate ratio ─────────────
set.seed(42)
n <- 400L
analytic <- data.frame(
  person_id       = seq_len(n),
  exposed         = rep(c(1L, 0L), each = n %/% 2L),
  person_years    = runif(n, 0.5, 2.0),
  age             = rnorm(n, 60, 10),
  sex             = rbinom(n, 1L, 0.5),
  baseline_comorb = rpois(n, 2)
)
analytic$rate     <- ifelse(analytic$exposed == 1L, 0.06, 0.03)
analytic$event_count <- rpois(n, analytic$rate * analytic$person_years)
analytic <- subset(analytic, person_years > 0)

pois_fit <- glm(
  event_count ~ exposed + age + sex + baseline_comorb,
  family = poisson(link = "log"),
  data   = analytic,
  offset = log(person_years)
)

# Cluster-robust variance by person_id
vc_cl  <- vcovCL(pois_fit, cluster = ~ person_id)
ct_cl  <- coeftest(pois_fit, vcov. = vc_cl)
irr    <- exp(ct_cl["exposed", "Estimate"])
irr_ci <- exp(ct_cl["exposed", "Estimate"] + c(-1, 1) * 1.96 * ct_cl["exposed", "Std. Error"])
cat(sprintf("Poisson rate model (cluster-robust SEs):\n"))
cat(sprintf("  IRR = %.3f  95%% CI: (%.3f, %.3f)\n", irr, irr_ci[1], irr_ci[2]))
disp <- sum(residuals(pois_fit, type = "pearson")^2) / df.residual(pois_fit)
cat(sprintf("  Pearson dispersion = %.3f  (equidispersion -> ~1.0)\n\n", disp))

# ── 4. Zou's modified Poisson: risk ratio on a binary outcome ─────────────────
analytic_bin <- data.frame(
  person_id = seq_len(n),
  exposed   = rep(c(1L, 0L), each = n %/% 2L),
  outcome   = rbinom(n, 1L, prob = rep(c(0.25, 0.15), each = n %/% 2L)),
  age       = analytic$age,
  sex       = analytic$sex
)

zou_fit <- glm(
  outcome ~ exposed + age + sex,
  family = poisson(link = "log"),   # Poisson family, log link, NO offset
  data   = analytic_bin
)
vc_hc1 <- vcovHC(zou_fit, type = "HC1")   # sandwich SEs correct Bernoulli variance
ct_hc1 <- coeftest(zou_fit, vcov. = vc_hc1)
rr     <- exp(ct_hc1["exposed", "Estimate"])
rr_ci  <- exp(ct_hc1["exposed", "Estimate"] + c(-1, 1) * 1.96 * ct_hc1["exposed", "Std. Error"])
cat("Zou modified Poisson (robust SEs) -- risk ratio on binary outcome:\n")
cat(sprintf("  RR = %.3f  95%% robust CI: (%.3f, %.3f)\n", rr, rr_ci[1], rr_ci[2]))
cat("  (Use instead of logistic OR when outcome prevalence > ~10%)\n")