Canonical link functions

Standard generalized linear models are often written using canonical link functions.
For an exponential family distribution, the canonical link is the function of the mean parameter $\mu$ that appears in the exponent of the density.

Although canonical links are convenient, they are not required. For example, in binomial and cumulative multinomial models, the probit link $g(\mu) = \Phi^{-1}(\mu)$ is often used even though it is not canonical.

Offsets

Sometimes it is useful to include a predictor whose coefficient is fixed and known, rather than estimated. Such a predictor is called an offset.

A common example arises in Poisson models. Suppose the rate of occurrence is $\mu$ per unit time. With exposure $T$, the expected count is $\mu T$. We may want to model
$\log \mu = X \beta$, as in the usual Poisson GLM.
However, the model is parameterized in terms of the mean of $y$, which is $\mu T$.

To fit this with a Poisson GLM, we write:

log⁡(μT)=log⁡μ+log⁡T=Xβ+log⁡T,\log(\mu T) = \log \mu + \log T = X \beta + \log T,

so we augment the design matrix with a column containing $\log T$, and fix its coefficient to 1. This augmented column is the offset. A typical application is Poisson regression for rates, where previous exposure or baseline rate is treated as an offset.

Interpreting model parameters

The choice and parameterization of predictors $x$ follow the same principles as in linear models. However, there is an important additional consideration: the linear predictor relates to the link-scale, not directly to the outcome scale.

The linear predictor is used to approximate $g(\mu)$, where $\mu = E(y)$. Therefore, the effect of changing the $j$-th predictor $x_j$ by a fixed amount depends on the current value of $x$, because the transformation back to the outcome scale is nonlinear.

One way to interpret effects on the scale of $y$ is to define a reference case with a predictor vector $x_0$ and corresponding mean

y0=g−1(x0β).y_0 = g^{-1}(x_0 \beta).

Then, for a change $\Delta x$ in the predictors, the new outcome is

ynew=g−1(g(y0)+Δx β),y_{\text{new}} = g^{-1}(g(y_0) + \Delta x \, \beta),

and the effect can be expressed as $y_{\text{new}} – y_0$. This formulation allows interpretation in terms of contrasts relative to a meaningful baseline case.

Latent continuous data interpretation

Discrete-data regression models can often be interpreted in terms of unobserved (latent) continuous variables, which is useful both conceptually and computationally.

Probit model for binary data

Consider the probit model for binary outcomes:

Pr⁡(yi=1)=Φ(Xiβ).\Pr(y_i = 1) = \Phi(X_i \beta).

This is equivalent to introducing a latent variable $u_i$ with

ui∼N(Xiβ,1),yi={1if ui>0,0if ui≤0.u_i \sim N(X_i \beta, 1), \quad y_i = \begin{cases} 1 & \text{if } u_i > 0, \\ 0 & \text{if } u_i \le 0. \end{cases}

In substantive applications, $u_i$ may have an interpretable meaning (for example, an underlying “preference” or “propensity”).

This parameterization also supports an efficient Gibbs sampler:

  • Conditional on $u_i$, the model is just a normal linear regression in $\beta$.
  • Conditional on $\beta$ and $y_i$, each $u_i$ follows a truncated normal:
    $u_i \sim N(X_i \beta, 1)$ truncated above or below $0$ depending on $y_i$.

Logistic model via latent variable

A similar construction can be considered for logistic regression. Instead of $u_i \sim N(X_i \beta, 1)$, one defines

ui∼Logistic(Xiβ,1),u_i \sim \text{Logistic}(X_i \beta, 1),

where the logistic distribution has density proportional to $1 / (e^{x/2} + e^{-x/2})$, the derivative of the inverse-logit function. This model is conceptually helpful but is less convenient computationally than the probit version, because linear regression with logistic errors is harder to handle.

Ordered multinomial models

For ordered categorical outcomes, say $y_i \in {0, 1, 2, 3}$, one can introduce latent variables $u_i$ and cutpoints $c_0, c_1, c_2$ such that:

  • $y_i = 0$ if $u_i < c_0$
  • $y_i = 1$ if $c_0 \le u_i < c_1$
  • $y_i = 2$ if $c_1 \le u_i < c_2$
  • $y_i = 3$ if $u_i \ge c_2$

There are multiple ways to parameterize the cutpoints, and the choice affects identifiability and hierarchical extensions. Because adding a constant to all cutpoints leaves the likelihood unchanged, one must fix one cutpoint, for example setting $c_0 = 0$, so that only the remaining cutpoints need to be estimated. The binary probit model is recovered as a special case with a single fixed cutpoint at 0.

Bayesian generalized linear models

We can consider several prior structures for $\beta$ and, if present, a dispersion parameter $\phi$.

Let $\phi$ denote a dispersion parameter (for example, variance in a normal model). A joint prior can be written as

p(β,ϕ)=p(ϕ) p(β∣ϕ).p(\beta, \phi) = p(\phi) \, p(\beta \mid \phi).

Noninformative priors for $\beta$

A flat or noninformative prior for $\beta$ yields a posterior mode at the maximum likelihood estimate. This mode can be computed by iteratively reweighted least squares (IRLS), as implemented in standard software.

A normal approximation around the mode provides approximate posterior inference.

Conjugate prior using hypothetical data

A conceptually useful method for specifying a prior for $\beta$ is via hypothetical data $(y_0, X_0)$, of size $n_0$, assumed to arise from the same generalized linear model:

  • Combine real and hypothetical data into a single stacked dataset:
    • Response: $(y, y_0)$
    • Design matrix: $(X, X_0)$
  • Use a flat prior on $\beta$.

The resulting posterior is equivalent to using a conjugate prior based on the pseudo-data $(y_0, X_0)$. Computationally, the same IRLS machinery used for noninformative priors can then be applied.

Nonconjugate priors

Often, prior information is more naturally expressed directly as a distribution on $\beta$, for example:

β∼N(β0,Σβ),\beta \sim N(\beta_0, \Sigma_\beta),

with specified $\beta_0$ and $\Sigma_\beta$. Normal priors are particularly convenient because they combine well with normal approximations to the likelihood.

Hierarchical models

Hierarchical priors for $\beta$ allow dependence on hyperparameters and naturally accommodate complex structured predictors. Normal hyperpriors are common, and one can treat generalized linear models similarly to hierarchical normal linear models by using a normal approximation to the generalized linear model likelihood.

This approach helps avoid overfitting while still allowing a rich set of predictors.

Normal approximation to the likelihood

For many generalized linear models, exact posterior computations are infeasible. A useful strategy is to approximate the likelihood in $\beta$ (and possibly $\phi$) by a normal distribution, using a local quadratic approximation.

Pseudodata and pseudovariances

For each observation $y_i$, we construct a pseudodatum $z_i$ and a pseudovariance $\sigma_i^2$ such that:

p(yi∣Xiβ,ϕ)≈N(zi∣Xiβ,σi2).p(y_i \mid X_i \beta, \phi) \approx N(z_i \mid X_i \beta, \sigma_i^2).

Stacking the $z_i$ and $\sigma_i^2$ yields an approximate linear regression problem:

  • Response vector $z = (z_1, \ldots, z_n)$
  • Design matrix $X$
  • Known variance matrix $\operatorname{diag}(\sigma_1^2, \ldots, \sigma_n^2)$

This lets us reuse standard linear regression algorithms.

Center of the approximation

The approximation is centered at some point $(\hat{\beta}, \hat{\phi})$, with corresponding linear predictor $\hat{\eta} = X \hat{\beta}$.

  • In the mode-finding phase, $(\hat{\beta}, \hat{\phi})$ is updated iteratively.
  • Once the mode is found, the normal approximation is evaluated at this fixed point.

Determining $z_i$ and $\sigma_i^2$

Let $L(y_i \mid \eta_i, \phi)$ denote the log-likelihood contribution of $y_i$. We approximate:

L(yi∣ηi,ϕ)≈−12σi2(zi−ηi)2+constant,L(y_i \mid \eta_i, \phi) \approx -\frac{1}{2 \sigma_i^2} (z_i – \eta_i)^2 + \text{constant},

which is a quadratic in $\eta_i$.

Using a Taylor expansion around $\hat{\eta}_i = X_i \hat{\beta}$, and writing derivatives as $L'(y_i \mid \eta_i, \phi)$ and $L”(y_i \mid \eta_i, \phi)$, we obtain:

zi=η^i−L′(yi∣η^i,ϕ^)L′′(yi∣η^i,ϕ^),σi2=−1L′′(yi∣η^i,ϕ^).z_i = \hat{\eta}_i – \frac{L'(y_i \mid \hat{\eta}_i, \hat{\phi})}{L”(y_i \mid \hat{\eta}_i, \hat{\phi})}, \quad \sigma_i^2 = -\frac{1}{L”(y_i \mid \hat{\eta}_i, \hat{\phi})}.

Example: binomial-logistic model

For binomial data with logistic link,

L(yi∣ηi)=yiηi−nilog⁡(1+eηi),L(y_i \mid \eta_i) = y_i \eta_i – n_i \log(1 + e^{\eta_i}),

we have:

dLdηi=yi−nieηi1+eηi,d2Ldηi2=−nieηi(1+eηi)2.\frac{dL}{d\eta_i} = y_i – n_i \frac{e^{\eta_i}}{1 + e^{\eta_i}}, \quad \frac{d^2L}{d\eta_i^2} = -n_i \frac{e^{\eta_i}}{(1 + e^{\eta_i})^2}.

Plugging into the formulas above gives the pseudodata $z_i$ and variances $\sigma_i^2$ as functions of $\hat{\eta}_i$ and $y_i$.

Iterative weighted least squares (IRLS)

The normal approximation suggests an algorithm:

  1. Start with an initial guess $(\beta, \phi)$.
  2. Compute $\hat{\eta}$, then $z_i$ and $\sigma_i^2$ for all $i$.
  3. Fit a weighted linear regression of $z$ on $X$ with weights $1 / \sigma_i^2$ to update $\beta$.
  4. If a dispersion parameter $\phi$ is present, update $\phi$ by maximizing its conditional posterior given the current $\beta$.
  5. Repeat until convergence.

This is equivalent to applying Newton’s method to the system $\partial p(\beta \mid y) / \partial \beta = 0$ and typically converges quickly.

One issue is that estimates of components of $\beta$ can diverge to $\pm \infty$ in situations like complete or quasi-complete separation in logistic regression. Reasonable prior information usually prevents this.

Approximate normal posterior distribution

Once the mode $(\hat{\beta}, \hat{\phi})$ is found, a normal approximation to the conditional posterior of $\beta$ is:

p(β∣ϕ^,y)≈N(β∣β^,Vβ),p(\beta \mid \hat{\phi}, y) \approx N(\beta \mid \hat{\beta}, V_\beta),

where

Vβ=(XT diag⁡{−L′′(yi∣η^i,ϕ^)} X)−1.V_\beta = \left(X^T \, \operatorname{diag}\{-L”(y_i \mid \hat{\eta}_i, \hat{\phi})\} \, X\right)^{-1}.

If there is no dispersion parameter (for example, binomial or Poisson models) and the sample size is large, this approximation is often adequate for posterior summaries.

If a dispersion parameter $\phi$ is present, one can approximate the marginal posterior $p(\phi \mid y)$ using a Laplace-type correction at the conditional mode $\hat{\beta}(\phi)$, with a similar normal approximation for $\beta$ given $\phi$.

Beyond simple normal approximations

The normal approximation is a useful starting point, but more advanced methods such as full Markov chain Monte Carlo, variational Bayes, or expectation propagation can yield more accurate approximations to the posterior distribution. For generalized linear models, these methods often still rely on local normal approximations to the likelihood, but they refine the approximation globally across parameter space.