PQL And Likelihood-Based Integration For Binary Models

Two Different Routes From a Binary Model to Estimates

For a binary outcome, the conditional mean is nonlinear in the random effects. That creates a practical estimation problem: the response is Bernoulli, but the random effects must be accounted for when estimating the fixed effects and their variability. Penalized quasi-likelihood (PQL) addresses this by repeatedly replacing the binary observations with a Gaussian working response. Likelihood-based integration instead averages the conditional binary likelihood over the random-effect distribution.

These are not merely two computational routes to the same objective. PQL fits a sequence of local, linearized working models. The integrated likelihood evaluates the probability of the observed binary data after averaging over possible random effects. Since linearization and integration generally do not commute, their estimates can differ, particularly when the binary mean changes sharply over the range of plausible random effects.

Model: For cluster \(j\) and observation \(i\), let \(Y_{ij}\) be binary, with \(Y_{ij}\mid b_j\) conditionally independent Bernoulli variables and \[ \operatorname{logit}(\mu_{ij})=\eta_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_j. \] Here \(\mu_{ij}=P(Y_{ij}=1\mid b_j)\), \(x_{ij}\) and \(z_{ij}\) are known covariate vectors, and \(b_j\sim N(0,D)\), independently across clusters. The conditional model and the random-effect distribution together specify the GLMM. The marginal probability of an outcome is obtained by averaging its conditional probability over the distribution of \(b_j\).

The assumptions in this formulation matter. Conditional on the random effects, the binary outcomes are independent and have the stated logistic means; the random effects follow the specified mean-zero multivariate normal distribution; and clusters are independent. The integrated likelihood below relies on these assumptions. PQL also relies on the conditional mean and variance to construct its working response, but then makes an additional Gaussian working-model approximation.

What PQL Linearizes

At a current estimate of the linear predictor \(\eta_{ij}\), write \(\mu_{ij}=\operatorname{logit}^{-1}(\eta_{ij})\). The derivative of the inverse-logit mean with respect to \(\eta_{ij}\) is \(\mu_{ij}(1-\mu_{ij})\). A first-order expansion of the mean around the current predictor gives \[ Y_{ij}-\mu_{ij}\ \approx\ \mu_{ij}(1-\mu_{ij})(\eta^{\mathrm{new}}_{ij}-\eta_{ij}). \] Rearranging defines the working response \(Y^*_{ij}\) and its weight \(w_{ij}\):

Y^*_{ij}=\eta_{ij}+\frac{Y_{ij}-\mu_{ij}}{\mu_{ij}(1-\mu_{ij})}, \qquad w_{ij}=\mu_{ij}(1-\mu_{ij}).

The working response is treated as approximately Gaussian with mean \(x_{ij}^{T}\beta+z_{ij}^{T}b_j\) and variance \(1/w_{ij}\). Consequently, one PQL iteration fits a weighted linear mixed model to \(Y^*\), using weights calculated from the current fitted means. Since \(Y^*\) and the weights change as the fitted means change, this is an iterative procedure, not a single fixed linear model.

With a Gaussian random-effect penalty, the corresponding working objective, up to constants, is

Q(\beta,b)=\frac{1}{2}\sum_j \left(Y^*_j-X_j\beta-Z_jb_j\right)^T W_j \left(Y^*_j-X_j\beta-Z_jb_j\right) +\frac{1}{2}\sum_j b_j^T D^{-1}b_j.

Here \(W_j\) is diagonal with entries \(w_{ij}\). The first term measures weighted disagreement with the working responses; the second penalizes random effects according to their assumed covariance. In practice, PQL alternates between updating working quantities and fitting the working mixed model, including its variance components. The objective is useful for understanding the update, but it is not the exact Bernoulli log-likelihood integrated over the random effects.

Worked Example: Minimizing One Working-Model Objective

Consider a single random-effect coordinate \(b\) in a hypothetical PQL iteration. Suppose the aggregate working-data contribution has weight 1 and working target 2, while the Gaussian penalty has precision 8. The part of the objective depending on \(b\) is therefore \[ Q(b)=\frac{1}{2}(2-b)^2+4b^2. \] These are working-model quantities for this calculation, not observed binary outcomes.

Differentiate and set the derivative to zero: \[ \frac{dQ}{db}=-(2-b)+8b=-2+9b=0, \qquad \hat b=\frac{2}{9}\approx 0.2222. \] The objective at \(b=0\) is \(Q(0)=\tfrac12(2)^2=2\). At the updated value, \[ Q\left(\frac{2}{9}\right) =\frac{1}{2}\left(\frac{16}{9}\right)^2+4\left(\frac{2}{9}\right)^2 =\frac{128}{81}+\frac{16}{81} =\frac{144}{81} =\frac{16}{9}\approx 1.778. \] Thus the update reduces this working objective. It does not establish that the integrated binary likelihood has increased: the two objectives are different.

Likelihood-Based Integration

For a given cluster, the conditional probability of its observed binary data at random effect \(b_j\) is \[ L_j(\beta\mid b_j)=\prod_i \mu_{ij}^{Y_{ij}}(1-\mu_{ij})^{1-Y_{ij}}. \] The likelihood contribution that does not condition on the unobserved random effect averages this quantity over its assumed distribution:

L_j(\beta,D)=\int L_j(\beta\mid b)\,\phi(b;0,D)\,db, \qquad L(\beta,D)=\prod_j L_j(\beta,D).

In this expression, \(\phi(b;0,D)\) is the multivariate normal density with mean zero and covariance \(D\). The integration is within each cluster because each cluster has its own random effect, shared by its observations. In general, the integral has no simple closed form for a logistic model, so numerical integration is used. Gaussian quadrature is one standard option. Crucially, it approximates an integral of the conditional likelihood; it does not replace the binary likelihood by a Gaussian one.

Worked Example: Integrating a Cluster’s Binary Likelihood

Suppose a hypothetical cluster has two successes, and consider the model with intercept \(\beta=0\), random intercept \(B\sim N(0,1)\), and no other covariates. Conditional on \(B=b\), each success has probability \(\operatorname{logit}^{-1}(b)\), so the probability of two successes is \(\operatorname{logit}^{-1}(b)^2\). The cluster’s marginal probability is \[ P(Y_1=1,Y_2=1)=E\!\left[\operatorname{logit}^{-1}(B)^2\right]. \]

For illustration, approximate this expectation with three-point Gauss-Hermite quadrature for a standard normal: use nodes \(0,\sqrt{3},-\sqrt{3}\) and weights \(2/3,1/6,1/6\). The resulting calculation is \[ \frac{2}{3}(0.5)^2+\frac{1}{6}\{\operatorname{logit}^{-1}(\sqrt{3})^2+ \operatorname{logit}^{-1}(-\sqrt{3})^2\}. \] Since \(\operatorname{logit}^{-1}(\sqrt{3})\approx0.8497\) and \(\operatorname{logit}^{-1}(-\sqrt{3})\approx0.1503\), the approximation is \[ \frac{2}{3}(0.25)+\frac{1}{6}(0.8497^2+0.1503^2) \approx 0.1667+0.1239=0.2906. \] At \(b=0\), the conditional probability of two successes is \(0.25\), not \(0.2906\). The marginal likelihood contribution averages over the full random-intercept distribution; evaluating the conditional likelihood only at its mean does not perform that integration. The quadrature value is an approximation, with its accuracy depending on the number and placement of quadrature points.

Why the Estimates Can Differ

PQL first transforms each response using the current fitted mean, then solves a weighted Gaussian mixed-model problem. Integration instead averages a nonlinear likelihood over the random-effect distribution. The discrepancy can be seen even at the level of a marginal mean: in general, \[ E\{\operatorname{logit}^{-1}(x^T\beta+Z^Tb)\} \ne \operatorname{logit}^{-1}(x^T\beta). \] The left side is an average of conditional probabilities; the right side evaluates the conditional mean at the mean random effect. Nonlinearity makes these distinct quantities.

Worked Example: Conditional Probability Versus an Integrated Mean

Take a hypothetical random-intercept model with \(B\sim N(0,1)\) and fixed linear predictor \(\beta=0.485\). The conditional probability at \(B=0\) is \[ \operatorname{logit}^{-1}(0.485)\approx0.6189. \] To approximate the probability averaged over random intercepts, apply the same three-point quadrature: \[ E\{\operatorname{logit}^{-1}(0.485+B)\} \approx\frac{2}{3}\operatorname{logit}^{-1}(0.485) +\frac{1}{6}\left\{\operatorname{logit}^{-1}(0.485+\sqrt{3}) +\operatorname{logit}^{-1}(0.485-\sqrt{3})\right\}. \] The three inverse-logit values are approximately \(0.6189\), \(0.9018\), and \(0.2231\). Thus the approximation is \[ \frac{2}{3}(0.6189)+\frac{1}{6}(0.9018+0.2231) \approx0.4126+0.1875=0.6001. \] The conditional probability at a zero random effect is about 0.619, whereas the random-effect-averaged probability is about 0.600. These are different summaries from the same specified model, not contradictory answers. The quadrature calculation illustrates why a method based on local conditional linearization can target results that differ from those based on the integrated likelihood.

The calculations also explain why a marginal probability alone should not be mistaken for a full integrated-likelihood estimate. A complete likelihood uses the joint pattern of outcomes within each cluster, not only an overall event proportion. In clustered data, within-cluster dependence carries information about the random-effect variation as well as the fixed effects.

Assumptions, Diagnostics, and Alternatives

Both approaches require a defensible conditional model, link, random-effect structure, and distributional assumptions. Integrated likelihood additionally requires accurate evaluation of the cluster integrals. PQL requires that its changing Gaussian working approximation be adequate for the estimation task; it is not exact merely because the working mixed model is fitted carefully.

Practical warning signs include extreme fitted probabilities, very small working weights \(\mu_{ij}(1-\mu_{ij})\), unstable variance-component updates, convergence problems, and marked changes in estimates across iterations or numerical settings. Small weights make the working response especially sensitive to a residual \(Y_{ij}-\mu_{ij}\), because that residual is divided by \(\mu_{ij}(1-\mu_{ij})\). Such warnings do not by themselves prove that an estimate is wrong, but they call for checking the fit and considering whether the approximation is suitable.

Approach	What is evaluated	Main approximation or requirement
PQL	A sequence of weighted Gaussian working-model objectives	First-order linearization and Gaussian working responses
Likelihood-based integration	The conditional binary likelihood averaged over random effects	Numerical integration must adequately approximate each integral

PQL can be computationally convenient, especially when the working mixed model is easy to fit. Likelihood-based integration directly respects the specified binary response model after accounting for random effects, but may require more computation and careful numerical integration. Neither method repairs a misspecified random-effects distribution, omitted clustering structure, or inadequate covariates. Integration also cannot remove the sampling uncertainty caused by limited data.

Common Mistake

Pitfall: “PQL fits a mixed model, so it maximizes the same likelihood as the binary GLMM.” This is incorrect. PQL fits an iteratively updated Gaussian working model based on a local approximation. The integrated binary likelihood instead averages the conditional Bernoulli likelihood over the random-effect distribution. The objectives, and therefore the resulting estimates, need not agree.

When reporting or comparing fitted results, first identify which probability is being described: conditional on a specified random effect, or averaged over its distribution. Then identify how the estimates were obtained. A difference between PQL and integrated-likelihood estimates is not automatically evidence of a coding error; it can arise because the procedures optimize or approximate different quantities.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a logistic binary GLMM, what distribution is integrated in a cluster’s marginal likelihood, and what is being averaged?
If \(\mu=0.2\), calculate the PQL working weight and the multiplier applied to \(Y-\mu\) in the working response.
Why does evaluating the conditional likelihood at \(b=0\) generally not equal the integrated likelihood contribution?
A colleague says that a lower PQL working objective proves a higher integrated binary likelihood. What is the error?

Solutions. (1) Integrate over the assumed distribution of the cluster’s random effect \(b_j\); the conditional likelihood of that cluster’s observed outcomes is averaged. (2) The weight is \(w=\mu(1-\mu)=0.2(0.8)=0.16\). The multiplier on \(Y-\mu\) is \(1/w=6.25\), so \(Y^*=\eta+6.25(Y-\mu)\). (3) The conditional likelihood is a nonlinear function of \(b\), so its average over the random-effect distribution generally differs from its value at the mean \(b=0\). (4) The PQL objective is a working Gaussian objective, whereas the integrated binary likelihood is a different quantity; reducing one does not establish that the other increased.