Reporting Binary PQL Estimates With Caveats

Why Reporting Needs More Than a Coefficient

A binary PQL fit can produce familiar-looking coefficients and standard errors, but they come from an approximation to the generalized linear mixed model, not an exact evaluation of its likelihood. A report that gives only an odds ratio or a p-value leaves readers unable to tell what outcome was modeled, how observations were grouped, or how much confidence to place in the approximation. The aim is to report the model and its estimates together, with qualifications that reflect the data and fitting diagnostics.

The relevant model is a binary random-intercept model. Let \(Y_{ij}\) be the outcome for individual \(i\) in cluster \(j\), with \(Y_{ij}=1\) for the event and \(Y_{ij}=0\) otherwise. Let \(x_{ij}\) denote the covariate vector, and let \(b_j\) be a cluster-specific random intercept.

These assumptions describe the statistical model; PQL adds an estimation approximation. In particular, the random intercepts are modeled as normally distributed and independent across clusters, and observations within a cluster are conditionally independent after accounting for the covariates and that intercept. The model does not say that observed binary outcomes are independent within a cluster marginally: shared \(b_j\) induces dependence. A fixed-effect coefficient \(\beta_k\) describes a change in conditional log odds for a one-unit increase in its covariate, holding the other covariates and cluster effect fixed.

What PQL Estimates and What to Report

PQL begins from the conditional linear predictor \(\eta_{ij}=x_{ij}^{\mathsf T}\beta+b_j\). At a current estimate, it uses a first-order expansion of the link relationship to form a working response and weights, then fits a weighted linear mixed model. For a binary logit model, the working response and weight at the current fitted mean \(\mu_{ij}\) are

The estimates are updated iteratively. This construction is useful computationally, but the linearization is least reassuring when fitted probabilities are near zero or one, or when clusters have few events. Therefore, the output should be identified as PQL, and conclusions should be framed as approximate model-based results. The earlier material on the role of the link function in PQL bias and PQL sensitivity analysis explains why link choice and alternative reasonable fits can matter; reporting should make the chosen specification and any material instability visible.

For a fixed-effect estimate \(\hat\beta_k\), the conditional odds ratio for a one-unit increase in \(x_k\) is \(\exp(\hat\beta_k)\), with other terms held fixed. If a model-based standard error \(SE(\hat\beta_k)\) is reported, a conventional approximate 95% Wald interval on the coefficient scale is \(\hat\beta_k\pm1.96SE(\hat\beta_k)\). Exponentiating both endpoints gives an approximate interval for the conditional odds ratio. Because the standard error and estimate come from a PQL fit, the interval is approximate in more than the usual sampling sense: it also inherits the fitting approximation.

For the random intercept, report \(\hat\tau^2\), or \(\hat\tau\), with its scale and any boundary or convergence issue. This is heterogeneity in the cluster-specific log odds, not a direct percentage of outcome variation. A sometimes-reported latent-scale summary is \(\hat\tau^2/(\hat\tau^2+\pi^2/3)\). It uses the variance \(\pi^2/3\) of a standard logistic latent error; it is a model-based latent-scale summary, not the ordinary correlation of the observed binary outcomes.

Describe Cluster Heterogeneity on Its Own Scale

A random-intercept variance is not a second fixed-effect coefficient. It summarizes the modeled spread of cluster-specific intercepts around their mean, on the log-odds scale. An estimate \(\hat\tau\) can also be described through \(\exp(\hat\tau)\): a one-standard-deviation increase in the random intercept multiplies the conditional odds by that factor, relative to a cluster at the mean random intercept, with covariates held fixed. This is a scale description, not a claim that every pair of clusters differs by exactly that amount.

The latent-scale summary can help readers compare model-based heterogeneity, but should be labeled carefully. It relies on the logistic latent-variable representation and does not replace reporting \(\hat\tau^2\) or \(\hat\tau\). In particular, do not call it the observed within-cluster correlation unless that correlation has been separately defined and estimated for the intended outcome and design.

Make the Approximation and Diagnostics Visible

A useful report separates the numerical result from the evidence about whether that result is stable enough to interpret. State which PQL implementation was used, including whether it was first-order or a specified higher-order version. Report the number of clusters and a useful description of cluster sizes, especially when sizes vary substantially. For a binary outcome, give the event count and, where informative, how events and non-events are distributed across clusters. These facts help readers judge how much information supports the fixed effects and the heterogeneity estimate.

Also report convergence status and material warnings. A claimed convergence is not sufficient if the algorithm stopped at an iteration limit, the variance estimate is at or near zero, or the fixed-effect standard errors are extremely large. Sparse events, quasi-separation, fitted probabilities near 0 or 1, and a small number of clusters can make the working approximation and variance estimation fragile. These are reasons to temper conclusions, not automatic proof that a particular estimate is wrong.

When the analysis has compared reasonable specifications or estimation approaches, summarize any consequential differences rather than selectively reporting the most favorable result. The earlier material on sensitivity analysis provides relevant context for such comparisons. If conclusions change materially across plausible fits, say so explicitly. Do not use an ordinary likelihood-ratio test or likelihood-based information criterion as though a PQL fit supplied an exact likelihood: PQL is based on an approximation and does not generally provide that ordinary likelihood basis.

Alternatives, Limitations, and Common Mistakes

The main alternative emphasized in the earlier material on likelihood-based integration for binary models is to estimate the same GLMM using an approach that numerically integrates over the random effects. Such methods can avoid PQL's particular linearization, though their computation and numerical behavior also require assessment. Comparing estimates across fitting approaches can reveal approximation sensitivity. A single result from either approach is not automatically trustworthy merely because it is available.

An ordinary GLM is another possible analysis, but it does not account for the dependence induced by a shared cluster effect unless dependence is handled in some other way. It therefore does not answer the same conditional mixed-model question. The choice between models must reflect the design and estimand, rather than treating a GLM estimate as a direct correction for PQL.

PQL's central limitation for binary data is that its working linear mixed model approximates a nonlinear conditional mean. The approximation can affect fixed-effect estimates, standard errors, and especially random-effect variance estimates. The direction and size of the difference are not guaranteed for every data set. A report should neither dismiss every PQL result nor present PQL output with the same unqualified certainty as an exact likelihood calculation.

• Mistake: reporting only an odds ratio. Without the event coding, link, cluster structure, and conditional interpretation, readers cannot identify what the odds ratio describes. Give the model context and label the odds ratio conditional on the random effects and other covariates.
• Mistake: treating the PQL interval as exact. A Wald interval is already an approximation, and the estimate and standard error also reflect PQL's approximation. Call it an approximate model-based interval and note material diagnostics or sensitivity.
• Mistake: interpreting a latent-scale summary as observed correlation. The formula uses a logistic latent variance; it is not the correlation between observed binary responses. Report the random-effect scale as well.
• Mistake: turning a non-significant test into proof of no effect. Failure to reject \(H_0\) does not establish that the coefficient is zero. Report the estimate and interval, then describe the uncertainty in context.