The Approximation Behind Penalized Quasi-Likelihood

Why PQL Uses a Working Model

A generalized linear mixed model (GLMM) connects a response distribution to a linear predictor that includes random effects. Its likelihood for the observed data is generally difficult because the random effects must be integrated out of a non-Gaussian conditional response likelihood. Penalized quasi-likelihood (PQL) replaces that calculation with a local approximation: it linearizes the conditional mean and treats a transformed response as Gaussian. The resulting Gaussian working mixed model retains the random-effects structure, but it is not the original model’s likelihood.

The distinction matters. PQL is not an exact integration method for the original GLMM, nor does it simply set random effects to zero. It substitutes a locally constructed Gaussian model for the original conditional response model, then uses mixed-model calculations on that substitute. The earlier material on Laplace approximation and adaptive quadrature concerns approximations or numerical evaluation of the original random-effects integral; PQL changes the likelihood contributions before fitting the working model.

The usual setup assumes a specified conditional response family and variance function, a differentiable inverse link, and normally distributed random effects. For the local construction below, the derivative \(h'(\eta)\) at the expansion point must be finite and nonzero. The approximation is most credible when the relevant linear predictors stay near their expansion points; it is not a claim that the original response distribution is actually Gaussian.

From the Conditional Mean to a Gaussian Working Response

At a current set of parameter values, denote the expansion point by \(\eta_{ij}^{(0)}\), and let \(\mu_{ij}^{(0)}=h(\eta_{ij}^{(0)})\). A first-order Taylor expansion gives

Since \(E(y_{ij}\mid b_j)=h(\eta_{ij})\), this says that, locally, the conditional mean changes linearly with the linear predictor. Rearranging the approximation suggests the working response \(q_{ij}=\eta_{ij}^{(0)}+(y_{ij}-\mu_{ij}^{(0)})/h'(\eta_{ij}^{(0)})\). Equivalently, because \(g'(\mu)=1/h'(g(\mu))\), this is \(q_{ij}=\eta_{ij}^{(0)}+g'(\mu_{ij}^{(0)})(y_{ij}-\mu_{ij}^{(0)})\).

To see why, take the conditional expectation of this response and use the Taylor approximation:

Its conditional variance, evaluated at the expansion point, is \(\operatorname{Var}(q_{ij}\mid b_j)\approx [g'(\mu_{ij}^{(0)})]^2\phi V(\mu_{ij}^{(0)})\). Thus \(q_{ij}\) is treated as a Gaussian observation centered approximately on the original linear predictor, with a variance determined by the original response variance and the link derivative. The normality is part of the working approximation, not an assumption that changes the original response family.

Which Likelihood and Integration Are Replaced?

Let \(f(y_{ij}\mid b_j;\beta)\) be the original conditional response density or probability mass function. Conditional independence gives the cluster contribution to the exact marginal likelihood as

where \(\varphi(b_j;0,G)\) is the multivariate normal density of the random effects. The full likelihood is the product of these cluster contributions. In a non-Gaussian GLMM, the product inside the integral is generally not a Gaussian function of \(b_j\). PQL does not evaluate this original integral exactly or merely approximate its numerical value. It first replaces the original conditional response contributions with Gaussian working contributions.

For fixed expansion points and working variances, write \(q_j\) for the vector of working responses, \(X_j\) and \(Z_j\) for the fixed- and random-effects design matrices, and \(R_j\) for the diagonal matrix of working variances. The working model is

Here \(o_j\) is the offset vector. Conditional on \(b_j\), the working response is Gaussian; after integrating the Gaussian random effects, its working marginal distribution is \(q_j\sim N(X_j\beta+o_j,\;Z_jGZ_j^T+R_j)\). This Gaussian integration is not the original GLMM integration: it integrates a different, locally constructed Gaussian model. The practical fitting criterion is often expressed as a penalized working criterion, with a quadratic penalty for \(b_j\) inherited from its normal distribution. Equivalently, for fixed working quantities, Gaussian mixed-model calculations can account for the random effects through the working marginal covariance.

When the Approximation Is Credible—and When to Be Cautious

The first-order expansion omits curvature in \(h(\eta)\). Its neglected remainder is small only when changes in \(\eta\) are sufficiently local and the inverse link is smooth over the relevant range. The working Gaussian approximation also summarizes a potentially discrete or skewed response by its mean and variance. These approximations can be poor when the mean function is sharply curved over plausible random-effect values, or when the working response is dominated by a few large residuals.

PQL is useful because a Gaussian mixed model is often easier to fit than the original GLMM’s integrated likelihood. It gives a transparent route to estimating fixed effects and random-effect variation while retaining the clustered covariance structure. But its computational convenience does not make its estimates exact maximum-likelihood estimates for the original GLMM. With a Gaussian response and identity link, the working model is exact; for non-Gaussian responses, it is generally only a local surrogate.

• Check the expansion region. Large changes between the current and fitted linear predictors mean the local tangent may not represent the mean function well across the relevant range.
• Check whether the working fit is stable. Strong sensitivity to starting values or to successive updates of expansion points is a warning that one local Gaussian representation may not describe the data reliably.
• Check the response and cluster structure. Sparse observations within clusters or highly variable random effects can make a local approximation based on limited information especially fragile.
• Compare methods when precision matters. Differences from likelihood-based methods using Laplace approximation or adaptive quadrature may indicate material approximation sensitivity; agreement is reassuring but does not prove correctness.

A common mistake is to say, “PQL approximates the integral in the exact GLMM likelihood, so it is just a faster integration rule.” That is not the defining operation. PQL first linearizes the conditional mean and constructs a Gaussian working response. Gaussian integration may then be used for that working model, but it is an integral for a substituted model, not a numerical evaluation of the original non-Gaussian likelihood.