Link Functions, Variance Functions, and Poisson Regression with Exposure

A generalised linear model separates three decisions that the linear model fuses: the distribution of the response, the function that maps the linear predictor to the mean, and the function that relates the mean to the variance. Choosing a Poisson response with a log link makes all three explicit, and the choice has consequences that a squared-loss model on the same data does not have: the fitted values cannot go negative, the coefficients multiply rather than add, and the coefficients depend on the variance function through the estimation equations.

Three components and the canonical link

The response $y_i$ is assumed to follow a distribution in the exponential family with mean $\mu_i$. The linear predictor is $\eta_i = x_i^\top\beta$, and the link $g$ satisfies $g(\mu_i) = \eta_i$. The variance is a function of the mean, $\operatorname{Var}(y_i) = \phi V(\mu_i)$, with $V$ the variance function fixed by the family and $\phi$ a dispersion parameter.

For the Poisson family, $V(\mu) = \mu$: the variance equals the mean. The canonical link is the one that makes the natural parameter linear, and for the Poisson that is the log, $g(\mu) = \log\mu$. With the canonical link, the sufficient statistic for $\beta$ is $X^\top y$ and the estimation equations reduce to $X^\top(y - \mu) = 0$. The log link also has a range that matches the parameter space: $\mu = e^{\eta}$ is positive for every real $\eta$, so no combination of covariates can produce a negative fitted rate, and no constraint is needed on the linear predictor.

The identity link with a Poisson family is a bad pairing for a different reason: it permits negative fitted means, and the likelihood is not defined there, so the fitted model can be outside the parameter space. The square-root link stabilises the variance, since $\operatorname{Var}(\sqrt{y}) \approx 1/4$ for large $\mu$, but has no direct interpretation for the coefficients and is undefined at zero counts. The log link is the default because it satisfies all three requirements at once: correct range, a sufficient statistic, and interpretable coefficients.

Multiplication is the natural scale for rates

With a log link, a unit increase in $x_j$ multiplies the fitted mean by $e^{\beta_j}$. The coefficient is an additive effect on the log scale and a multiplicative effect on the rate scale, and reporting it as either without stating which is the most common source of misinterpretation. A coefficient of $0.69$ is a rate ratio of about 2, meaning the rate doubles; a coefficient of $0.05$ is a rate ratio of about 1.05, a five percent increase.

The multiplicative reading matches how counts arise. If each of $n$ exposures independently produces an event with a small probability $p$, the expected count is $np$ and the logarithm is $\log n + \log p$, an additive decomposition. A covariate that changes the per-exposure probability by a factor multiplies the expected count, so a model on the log scale has the right structure and a model on the count scale does not.

The variance function compounds this. For an observation with mean $\mu$, a change in the linear predictor has a larger effect on variance than on the mean in absolute terms, because both are equal to $\mu$. Fitting therefore weights observations differently from a least-squares fit: the estimation equations are $X^\top W X \beta = X^\top W z$ with $W = \operatorname{diag}(\mu_i)$, so high-rate observations, which have more variance, receive proportionally more weight. A squared-loss fit gives every observation the same weight regardless of its mean, which is inconsistent with a count model whose noise grows with the level.

import numpy as np

def poisson_irls(X, y, exposure=None, iters=50, tol=1e-9):
    n, p = X.shape
    off = np.zeros(n) if exposure is None else np.log(exposure)
    beta = np.zeros(p)
    for _ in range(iters):
        eta = X @ beta + off
        mu = np.exp(eta)                 # can overflow if eta is large
        W = mu
        z = eta - off + (y - mu) / W     # working response on the linear scale
        sw = np.sqrt(W)
        new = np.linalg.lstsq(X * sw[:, None], z * sw, rcond=None)[0]
        if np.max(np.abs(new - beta)) < tol:
            return new
        beta = new
    return beta

The exponentiation is the failure point. A linear predictor above about 709 in double precision produces an infinite mean, and an infinite mean produces an infinite working response and NaN in the next iteration. Clipping $\eta$ bounds the damage but introduces a bias in exactly the region where the fit is struggling; the honest response is to check for divergence and reduce the step, which is the same damping logic as in Newton and IRLS fitting for the logistic model.

Exposure is an offset, not a predictor

A count observed over a period of observation time, an area, or a population is not directly comparable to a count observed over a different period. A region with a larger population has more cases regardless of its rate. The quantity of interest is the rate $\mu/t$, not the count $\mu$, and the model is

$$\log \mu_i = \log t_i + x_i^\top \beta .$$

The term $\log t_i$ is an offset: a predictor whose coefficient is fixed at 1. It belongs on the right-hand side with no coefficient to estimate, because the relationship between exposure and the expected count is exact by construction, not an empirical finding.

Entering $\log t_i$ as an ordinary predictor is a specific and common error. The model then estimates a coefficient for it, and the fitted coefficient is generally not 1. It is whatever value best fits the sample, and in small samples or with correlation between exposure and the covariates it can be far from 1. With the coefficient free, predictions for a new exposure are scaled by the fitted value rather than the true proportionality, and the other coefficients absorb part of the exposure effect. The result is a model that fits the training data and predicts rates incorrectly for a differently-sized unit.

import numpy as np

def rate_ratio(beta, X0, X1):
    """Predicted rate ratio between two covariate vectors, exposure-independent."""
    return np.exp((X1 - X0) @ beta)

# The offset cancels in a ratio only when the two units share the same exposure;
# otherwise the comparison must be of mu/t, not of mu.

The offset also clarifies what “rate” means for interpretation. Coefficients describe the effect of a covariate on the rate per unit exposure, and a coefficient from a model without an offset is an effect on the count, which depends on the exposure distribution in the sample. Two studies of the same phenomenon with different exposure ranges and no offset will produce different coefficients for the same covariate, with no modelling error in either.

Overdispersion and what to do about it

Poisson regression assumes $V(\mu) = \mu$ exactly. When the observed variance exceeds the mean, which is the common case for real count data because unobserved heterogeneity, clustering, or dependence between events inflates the variance, the model is misspecified in the variance function while remaining consistent for the mean. The coefficient estimates still converge to the true values under correct specification of the mean, but the standard errors computed from the inverse information matrix are too small by a factor of roughly the square root of the dispersion, and confidence intervals are correspondingly too narrow.

Estimating the dispersion from the Pearson residuals, $\hat\phi = \frac{1}{n-p}\sum (y_i - \hat\mu_i)^2/\hat\mu_i$, and multiplying the covariance matrix by it is the quasi-Poisson correction. It leaves the coefficients unchanged and widens the intervals. A likelihood-based comparison of models is no longer available because the quasi-likelihood is not a likelihood, so nested-model tests have to use an $F$ or Wald approximation rather than a deviance difference.

The alternative is to change the family. The negative binomial adds a dispersion parameter through a gamma-distributed mixing effect, giving $V(\mu) = \mu + \mu^2/\theta$; as $\theta \to \infty$ it approaches the Poisson. Fitting it requires estimating $\theta$ as well as $\beta$, which is typically done by alternating between the mean fit and a moment or likelihood update for the dispersion parameter. The variance is now quadratic in the mean, so the weighting of observations changes and the fitted coefficients differ from the Poisson fit, unlike the quasi-Poisson correction which leaves them alone.

The choice between the two is a choice about what to assume. Quasi-Poisson assumes the mean model is right and only the variance is wrong, and it takes no position on the form of the excess variance. Negative binomial assumes the excess variance has a specific quadratic form, and in exchange it provides a genuine likelihood for model comparison and for prediction intervals. Neither repairs a misspecified mean, and a mean model missing a covariate that modifies the rate will show up as overdispersion regardless of which variance function is used.

Zero inflation is a separate pathology that looks like overdispersion when the counts are small. Excess zeros can arise from a process that is structurally incapable of producing a positive count, which is a mixture and not a variance-function problem; the mixture models it with a second component and estimates the mixing probability alongside the rate. Fitting a negative binomial to such data absorbs the excess zeros into a large dispersion parameter and produces rate estimates that are pulled toward zero, since the mean model is forced to explain the zeros as small positive rates.

When the mean model is the weak link

Poisson regression is often described as requiring the mean to be small, but no such requirement is in the model; the constraint is that the variance equals the mean, and that fails at large means for reasons unrelated to the mean’s magnitude. What does break at large means is the numerical conditioning of $X^\top W X$, since the weights are the fitted means and a handful of large counts can dominate the information matrix, effectively reducing the sample size to the number of high-count observations. This is the same phenomenon as the vanishing weights at saturation in the logistic model, inverted: there the informative observations were the uncertain ones, and here they are the large ones.

A rate model is also not a model of proportions. When the count has an upper bound, as in the number of successes in a fixed number of trials, the variance is $n\pi(1-\pi)$, which is not the mean, and the Poisson fit on the raw counts will be wrong in the direction of overstating the variance for high-probability trials. The binomial family with a logit link is the correct one, and the link question becomes the same one as in logistic regression as conditional likelihood. The general lesson is that the family, not the link, determines the variance, and choosing the link for its range properties while ignoring the variance function produces a model that fits the mean and gets its uncertainty wrong.

The separation of mean and variance also decides how the fit should be computed. Both the logistic model and the Poisson model are canonical-link GLMs fitted by the same weighted least-squares iteration, with the weights supplied by the variance function, so the difference between them is entirely in $V(\mu)$ and hence in which observations drive the estimate. That shared structure is the reason Newton and IRLS fitting applies unchanged when the family changes, and why a solver written for one canonical-link model transfers to the other with a single substitution.