The Role Of The Link Function In PQL Bias

Why the Link Function Matters in PQL

In a binary generalized linear mixed model, the link function maps a conditional event probability to a linear predictor. In penalized quasi-likelihood (PQL), this mapping affects more than how coefficients are described: it determines the local slope of the conditional mean curve, and therefore the working response and working weight used to construct the approximate mixed-model problem. Logit and probit links can consequently produce different working-model updates even when they describe the same conditional probability at a particular point.

This comparison helps explain one source of variation in PQL behavior. It does not imply that changing links eliminates PQL bias. PQL still substitutes a local linear approximation for a nonlinear conditional mean in a model with random effects; the quality of that approximation depends on the predictor values and random-effect variation as well as on the link.

Definition: For a binary observation $Y_{ij}$, let $Y_{ij}\mid b_j\sim\operatorname{Bernoulli}(p_{ij})$, where $p_{ij}=h(\eta_{ij})$ and $\eta_{ij}=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j$. Here $h$ is the inverse link, $x_{ij}$ and $z_{ij}$ are the fixed- and random-effect design vectors, and the random effects are conventionally modeled as $b_j\sim N(0,D)$. Conditional independence of observations given the random effects is assumed. For a current finite predictor $\eta_{ij}$, the PQL working response and weight are defined when $h'(\eta_{ij})\ne 0$: $$ \widetilde{Y}_{ij}=\eta_{ij}+\frac{Y_{ij}-p_{ij}}{h'(\eta_{ij})}, \qquad w_{ij}=\frac{\{h'(\eta_{ij})\}^{2}}{p_{ij}(1-p_{ij})}. $$

The weight is the local Bernoulli information for the linear predictor: it combines the change in mean per unit change in predictor, $h'(\eta_{ij})$, with the conditional Bernoulli variance $p_{ij}(1-p_{ij})$. In the working mixed model, larger weights give observations more influence on the local fit, all else equal. These are observation-specific quantities, recalculated at the current predictor during the PQL iterations.

How the Two Links Enter the Working Model

For the logit link, $h(\eta)=\{1+\exp(-\eta)\}^{-1}$. Its derivative is $h'(\eta)=p(1-p)$, so the working weight simplifies to $w_{\text{logit}}=p(1-p)$. For the probit link, $h(\eta)=\Phi(\eta)$, where $\Phi$ is the standard normal cumulative distribution function, and $h'(\eta)=\phi(\eta)$, the standard normal density. Its working weight is $w_{\text{probit}}=\phi(\eta)^2/\{\Phi(\eta)[1-\Phi(\eta)]\}$. Both derivatives are positive at every finite predictor, so the working responses and weights are defined there.

Equal predictor values do not generally mean equal conditional probabilities across links. At $\eta=0$, however, both links give $p=0.5$, making this a useful direct comparison. Away from zero, it is often more informative to compare the links at the same probability, which typically requires different predictor values. Either comparison must state what is being held fixed.

Worked Example: Comparing Weights at a Conditional Probability of 0.5

Suppose a hypothetical binary observation has current predictor $\eta=0$. Under both links, $p=0.5$, and the Bernoulli variance is $0.5(1-0.5)=0.25$. The conditional model and finite predictor assumptions hold; both inverse links are differentiable with nonzero derivative here.

Logit. Its derivative at zero is $0.5(1-0.5)=0.25$. Thus $$ w_{\text{logit}}=\frac{0.25^2}{0.25}=0.25. $$

Probit. Its derivative at zero is $\phi(0)=1/\sqrt{2\pi}\approx0.3989423$. Therefore $$ w_{\text{probit}}=\frac{\phi(0)^2}{0.25} =\frac{0.1591549}{0.25} \approx0.6366. $$

At this common conditional probability, the probit working weight is about $0.6366/0.25=2.546$ times the logit weight. This is a local working-model comparison, not evidence that one link is universally more precise or that its PQL estimates are less biased.

Worked Example: Working Weights and Responses Away From the Center

Consider a hypothetical observation at current predictor $\eta=1$. This time, the links imply different probabilities, so the results compare each link at the same predictor rather than at the same mean. The Bernoulli variance and nonzero derivative conditions hold for both links at this finite value.

Logit. The probability is $p=1/(1+e^{-1})\approx0.7311$, and $h'(\eta)=p(1-p)\approx0.1966$. The weight is $$ w_{\text{logit}}=\frac{0.1966^2}{0.7311(1-0.7311)} =\frac{0.03866}{0.1966} \approx0.1966. $$ If the observed outcome is $Y=1$, the working response is $$ \widetilde{Y}=1+\frac{1-0.7311}{0.1966} \approx1+1.3679=2.368. $$

Probit. Here $p=\Phi(1)\approx0.8413$, $h'(\eta)=\phi(1)\approx0.2420$, and $p(1-p)\approx0.1335$. Thus $$ w_{\text{probit}}=\frac{0.2420^2}{0.8413(1-0.8413)} \approx\frac{0.05855}{0.1335} \approx0.4386. $$ For $Y=1$, the working response is $$ \widetilde{Y}=1+\frac{1-0.8413}{0.2420} \approx1+0.6557=1.656. $$

At the same predictor, these links give different conditional probabilities, weights, and working responses. Such differences are why the link matters to PQL’s working fit. This calculation alone does not establish which link is appropriate; that choice concerns the conditional mean model and scientific interpretation as well as the approximation.

Linearization Error Depends on the Link

The working response comes from a first-order Taylor approximation. For a predictor change $\delta$, the exact mean is $h(\eta+\delta)$, while its local linear approximation around $\eta$ is $h(\eta)+h'(\eta)\delta$. Provided $h$ is twice differentiable between these predictor values, Taylor’s theorem gives $$ h(\eta+\delta)=h(\eta)+h'(\eta)\delta+\frac{1}{2}h''(\xi)\delta^2 $$ for some $\xi$ between $\eta$ and $\eta+\delta$. Thus the error is governed by local curvature and the size of the predictor change; a small first-order increment does not guarantee a small error when the curve bends appreciably or $\delta$ is large.

For logit, $h''(\eta)=p(1-p)(1-2p)$. For probit, $h''(\eta)=-\eta\phi(\eta)$. Both second derivatives are zero at $\eta=0$, so the quadratic term vanishes there; higher-order terms can still differ. This makes the center a useful illustration of how even curves with zero local curvature can depart differently from their tangent lines over a finite change.

Worked Example: Comparing a One-Unit Linearization

At $\eta=0$, take a hypothetical predictor change $\delta=1$. The tangent approximation for each link is its value at zero plus its derivative at zero times one. The conditional mean curves are differentiable, and their derivatives at the expansion point are nonzero.

Logit. The tangent approximation is $0.5+0.25(1)=0.7500$. The exact value is $h(1)=0.7311$. The signed error, exact minus approximate, is $$ 0.7311-0.7500=-0.01894. $$

Probit. The tangent approximation is $0.5+\phi(0)(1)=0.8989$. The exact value is $\Phi(1)=0.8413$. The signed error is $$ 0.8413-0.8989=-0.05760. $$

For this one-unit change, the probit tangent overstates the exact probability by about 0.05760, compared with about 0.01894 for logit. These are errors for this specified expansion point and change, not general rankings of the links. Different predictor locations or smaller changes can yield different comparisons.

What the Comparison Does—and Does Not—Say About PQL Bias

PQL uses the working quantities to fit a penalized working mixed model, then updates the predictors and repeats. The link affects those steps through $h$, $h'$, and $w$. It also defines the scale on which fixed and random effects shift the conditional mean. Consequently, two fitted models with different links need not have directly comparable coefficient magnitudes: a coefficient is a change in the linear predictor per unit of its covariate, and that predictor has a link-specific scale.

PQL bias, however, is not determined by the working weight at one point. It arises in an approximation involving nonlinear conditional means and random effects, and the approximation is revisited across observations and iterations. The differences shown here help locate how link choice changes the approximation; they do not establish a bias correction. The logit link is not automatically safe because its weights have a simpler form, and the probit link is not automatically preferable because it has a different local slope.

A useful comparison states the target explicitly: same predictor, same conditional probability, or models fitted to the same hypothetical data. The first two answer distinct mathematical questions. A full fit also depends on covariates, cluster structure, random-effect variation, and the realized outcomes, so isolated weight calculations cannot predict the direction or magnitude of all coefficient or variance-component biases.

Diagnostics, Alternatives, and Limits

When comparing links in a PQL analysis, inspect the quantities that generate its working approximation, not just whether the algorithm converged. Useful checks include the range of fitted predictors, the distribution of conditional fitted probabilities, and the working weights over observations and iterations. Plotting observed outcomes or residual summaries against the fitted predictor can help reveal regions where a fitted mean curve is poorly aligned with the data. Large predictor changes associated with appreciable curvature warn that a local tangent may be a weak representation over the relevant range.

These checks are warnings, not tests that prove PQL bias or identify its direction. Very small weights can make some observations contribute little to a particular working update; near extreme probabilities, both links may have small derivatives, so numerical behavior also deserves attention. Convergence only indicates that the iterative algorithm met its stopping rule. It does not certify that the PQL approximation is accurate.

The main alternative is to estimate the model using likelihood-based integration over the random effects, for example with Laplace approximation or adaptive quadrature. These methods address the random-effect integration differently rather than relying on the same PQL working-response construction. They also have their own computational and approximation considerations. Comparing PQL fits across links is a sensitivity analysis, not a substitute for checking the adequacy of the approximation or considering an integration-based analysis where feasible.

The comparison has further limits. The formulas above concern individual Bernoulli observations with a correctly specified conditional Bernoulli mean-variance relationship. Grouped or weighted observations require corresponding adjustments to the working quantities. Results at one predictor value do not summarize a whole cluster, and a chosen link should not be changed solely to obtain a preferred variance estimate. The conditional model assumptions—including the random-effects distribution and conditional independence—remain substantive assumptions, not consequences of selecting logit or probit.

Common Mistake

Pitfall: “The probit weight is larger at $p=0.5$, so probit must reduce PQL bias.” This confuses a local working weight with the total approximation error. A weight describes local information for the linear predictor at a specified point; PQL bias also depends on the nonlinear mean over the random-effect distribution and on the full iterative fit. The correction is to use weight and linearization comparisons to understand sensitivity, then assess the fitted model and approximation without treating either link as a remedy.

Key takeaway: Logit and probit links produce different local working weights and can produce different tangent-line errors. These differences explain how the link enters binary PQL, but neither a favorable weight nor a smaller error in one example guarantees less PQL bias overall.

Check Your Understanding

Answer each question, then compare with the solutions below.

At $\eta=0$, why can logit and probit be compared at the same conditional probability?
For a binary observation with $p=0.4$ and logit link, calculate its working weight.
What feature of the inverse-link curve, besides the current mean, contributes to first-order linearization error?
Why does a larger probit working weight at $p=0.5$ not prove that probit PQL has less bias?

Solutions. (1) Both inverse links give $p=0.5$ at $\eta=0$. (2) Under logit, $h'=p(1-p)=0.4(0.6)=0.24$, so $w=0.24^2/0.24=0.24$. (3) Curvature and the size of the predictor change determine how far the exact mean departs from its tangent approximation. (4) A local weight does not measure the full nonlinear random-effect approximation or establish the bias of the fitted model.