Why Sensitivity Analysis Matters for Binary PQL
Penalized quasi-likelihood (PQL) is an approximation for fitting generalized linear mixed models. With binary outcomes, the approximation can be sensitive to sparse events, cluster structure, and the chosen random-effects specification. A single converged fit therefore does not show that its fixed effects or variance components are stable. A sensitivity analysis asks whether the substantive conclusion persists across a small set of scientifically plausible choices.
The goal is not to search for the specification that produces the preferred result. Specify the reference model and a limited set of alternatives before interpreting their estimates. Keep the data, outcome definition, and target comparison fixed whenever possible. When a refit omits or reweights clusters, label it as an influence analysis because it no longer represents the same complete dataset.
Set Up a Comparable Set of Fits
For a binary response \(Y_{ij}\) from observation \(i\) in cluster \(j\), a common starting point is a conditional logistic random-intercept model:
This model assumes the specified conditional mean and logit link are appropriate, observations are conditionally independent given the random effects and covariates, and the cluster effects follow the stated distribution with variance \(\tau^2\). The PQL fit additionally relies on an approximation to the binary mixed-model estimation problem. These are modeling assumptions, not facts guaranteed by obtaining output.
Choose a reference model based on the scientific question and design. Candidate alternatives might add a prespecified confounder, replace a common slope with a scientifically justified random slope, or separate within-cluster and between-cluster exposure information. Do not vary several decisions at once if the purpose is to learn which choice matters. Use the same observations, outcome coding, reference category, and predictor scaling across fits unless a specific check requires changing them. Record convergence and warnings for every fit.
Cluster summaries help plan and interpret the comparisons. At minimum, examine each cluster's size \(n_j\), event count \(E_j=\sum_iY_{ij}\), and event proportion \(E_j/n_j\). Note all-zero and all-one clusters, very small clusters, and clusters with unusual covariate patterns. These summaries do not justify replacing the individual binary responses with proportions in an ordinary Gaussian analysis. Their role here is to reveal where the data may challenge the approximation or influence the fitted model.
Compare Estimates on Meaningful Scales
For a coefficient \(\hat{\beta}_1\) on a logit scale, \(\exp(\hat{\beta}_1)\) is the conditional odds ratio for a one-unit increase in \(x\), holding the other modeled quantities fixed. Compare the coefficient itself and, when useful, its odds ratio. A change in sign or a substantial change in magnitude may matter even if both fits yield intervals that include zero. Do not treat overlap or non-overlap of separate intervals as a formal test of the difference between fits.
A simple descriptive comparison with reference estimate \(\hat{\beta}_{1,R}\) is
where \(A\) denotes an alternative fit and \(SE_R\) is the reference standard error. \(D_\beta\) expresses the shift in units of the reference standard error; it is a scale for describing sensitivity, not a test statistic with a universal cutoff or a p-value. For a positive reference variance estimate \(\hat{\tau}^2_R\), the relative change is
This ratio is undefined when the reference estimate is zero and can exaggerate changes when it is very small. Report the variance estimates themselves as well. If an alternative adds a random slope, its variance components do not have a one-to-one comparison with a random-intercept-only model: the model has changed how heterogeneity is represented.
Worked Example: Comparing Three Plausible Specifications
Suppose a hypothetical binary-outcome data set has 24 clinics, 15 observations per clinic, and an exposure measured in one-unit increments. The outcome totals 96 events among 360 observations. Cluster summaries show two clinics with no events, one with all 15 outcomes, and the remaining clinics with a range of event counts. A logistic random-intercept PQL model is the reference. Two prespecified alternatives add an important measured confounder and then add a random exposure slope. Assume all three fits converge with no reported numerical warnings.
| Specification | Exposure estimate | SE | Random-intercept variance | Random-slope variance |
|---|---|---|---|---|
| Reference: random intercept | -0.55 | 0.24 | 0.72 | Not fitted |
| Add measured confounder | -0.43 | 0.25 | 0.51 | Not fitted |
| Random intercept and slope | -0.31 | 0.29 | 0.44 | 0.18 |
In the reference model, the exposure odds ratio is \(\exp(-0.55)=0.577\). After adding the confounder it is \(\exp(-0.43)=0.651\), approximately. The coefficient shift is \(-0.43-(-0.55)=0.12\), so \(D_\beta=0.12/0.24=0.50\). The direction remains negative, but the estimated association moves toward no conditional association.
The random-intercept variance changes from 0.72 to 0.51 in the confounder model. Its relative change is \((0.51-0.72)/0.72=-0.292\), a decrease of about 29.2%. The random-slope model gives an exposure estimate 0.24 higher than the reference, or \(0.24/0.24=1.00\) reference standard error. Its exposure odds ratio is \(\exp(-0.31)=0.733\). This is a noticeable specification-dependent shift, not evidence that one fit is correct. The random-intercept variance of 0.44 in this model should not be treated as directly equivalent to the reference variance because the added slope reallocates heterogeneity.
Use Cluster Summaries to Guide Influence Checks
A cluster-summary table can indicate whether a small number of clusters drive an apparent result. A targeted leave-one-cluster-out analysis refits the same model after removing one cluster at a time. This answers a limited question: how much does the fitted estimate change when that cluster is absent? It does not estimate the effect of collecting more clusters, and it does not license deleting a cluster from the primary analysis merely because its omission is convenient.
A useful justification for this check is direct: if a cluster carries unusual outcome or covariate information, removing it changes the empirical information contributing to the fit. Comparing the estimates across these refits reveals whether the full-data result is concentrated in that information. The check is descriptive; refits are dependent because they share most of their observations.
Worked Example: Checking a Large, High-Event Cluster
Consider a hypothetical study with eight clinics. Their sizes are 4, 5, 6, 8, 10, 12, 15, and 20, and their event counts are 0, 1, 1, 2, 4, 5, 9, and 16, respectively. Thus \(n=80\) and the total event count is \(0+1+1+2+4+5+9+16=38\), giving an overall event proportion \(38/80=0.475\). The largest clinic has an event proportion \(16/20=0.80\), while the smallest has \(0/4=0\).
Suppose the reference random-intercept PQL fit gives an exposure coefficient of \(-0.48\), \(SE=0.22\), and random-intercept variance \(0.64\). Omitting the largest clinic gives \(-0.29\), \(SE=0.24\), and variance \(0.39\). The coefficient shift is \(-0.29-(-0.48)=0.19\), which is \(0.19/0.22=0.864\) reference standard errors. The corresponding odds ratios are \(\exp(-0.48)=0.619\) and \(\exp(-0.29)=0.748\).
The variance change is \((0.39-0.64)/0.64=-0.391\), or a decrease of about 39.1%. Together, these shifts suggest that the full-data estimates are sensitive to the largest clinic. The appropriate conclusion is to investigate the clinic's covariate and design context and report the sensitivity; it is not to discard the clinic without a defensible data-quality or design reason.
When the Cluster Mean Changes the Question
Cluster summaries can also reveal a modeling distinction, rather than merely identify influential clusters. If exposure varies both within and between clusters, one common slope can combine two different associations. For a cluster-varying predictor, define the cluster mean \(\bar{x}_j\), the overall mean \(\bar{x}\), the within-cluster deviation \(x_{ij}-\bar{x}_j\), and the centered cluster mean \(\bar{x}_j-\bar{x}\). Including both terms allows their coefficients to differ. This is a change in the estimand and model specification, so it should be interpreted as such.
Worked Example: Separating Within- and Between-Clinic Exposure
Suppose a hypothetical data set has four clinics with five patients each. The numbers exposed are 0, 1, 4, and 5. The cluster exposure means are therefore 0, 0.2, 0.8, and 1.0; across the balanced 20-patient data set the overall exposure mean is \(10/20=0.50\). A reference model uses a single exposure coefficient and estimates \(\hat{\beta}=0.60\), with random-intercept variance 0.80. A prespecified within-between model uses \(x_{ij}-\bar{x}_j\) and \(\bar{x}_j-0.50\), estimating coefficients 0.35 and 1.10, respectively, with random-intercept variance 0.63.
The one-slope odds ratio is \(\exp(0.60)=1.82\) per one-unit increase in exposure under that model. In the decomposed model, the within-clinic odds ratio is \(\exp(0.35)=1.42\) per one-unit within-clinic increase, while the between-clinic odds ratio is \(\exp(1.10)=3.00\) per one-unit increase in the clinic mean. These coefficients answer distinct conditional questions. The variance changes by \((0.63-0.80)/0.80=-0.2125\), a decrease of 21.3%. That change may reflect both the different exposure representation and the PQL fit; it does not by itself show which model is preferable.
Diagnostics, Alternatives, and Limits
Sensitivity comparisons are most informative when every fit is technically credible. Record convergence status, warnings, boundary or near-zero variance estimates, very large standard errors, and any evidence of separation or unstable coefficients. A model that fails to converge is not a useful numerical competitor. Refit it only after checking data coding, scaling, and whether the candidate structure is supported by the number of clusters and observations. If a result changes materially across reasonable fits, describe the conclusion as specification-sensitive rather than selecting the most favorable estimate.
A compact sensitivity table should include the candidate specification, the exposure estimate and standard error, the corresponding odds ratio when helpful, relevant variance components, convergence notes, and the cluster-summary check performed. Keep interpretation tied to the scientific scale: distinguish a direction change from a magnitude change, and distinguish uncertainty from model-to-model variation. The standard alternative of likelihood-based integration can provide a useful comparison in appropriate settings, but it is a different fitting approach, not a guarantee that every discrepancy is approximation bias.
The main limitation is that a set of plausible fits cannot reveal the true parameter or quantify all uncertainty. Similar estimates across candidates are reassuring about those particular choices, not proof of correctness. Conversely, large differences flag dependence on specification or influential data structure; they do not establish that PQL alone caused the difference. Cluster summaries are diagnostic descriptions, not substitutes for the individual-level response model.
State the binary outcome, conditional model, link, random effects, and target coefficient.
Tabulate sizes, event counts, event proportions, and notable covariate patterns before selecting targeted checks.
Change one scientifically defensible specification at a time and record why it is plausible.
Report estimate shifts, variance changes, convergence, and influence findings without treating sensitivity measures as hypothesis tests.
Check Your Understanding
Answer each question, then check the solutions below.
- A reference coefficient is \(-0.40\) with \(SE=0.20\), and an alternative estimate is \(-0.10\). Find the shift in reference standard-error units and describe what it says.
- A random-intercept variance changes from 0.50 to 0.35. Calculate the relative change and interpret it.
- Why should a leave-one-cluster-out result not automatically replace the full-data estimate?
- In a within-between model, what does the coefficient of the cluster mean represent that a within-cluster deviation coefficient does not?
Solutions. (1) The shift is \(-0.10-(-0.40)=0.30\), and \(0.30/0.20=1.50\) reference standard errors. This is a sizable descriptive shift; it is not a p-value or formal test. (2) \((0.35-0.50)/0.50=-0.30\), a 30% decrease in the estimated variance. (3) The omitted-cluster fit uses a changed dataset and is an influence check; omission needs a substantive justification, not merely a preferred estimate. (4) The cluster-mean coefficient describes the modeled between-cluster association, while the deviation coefficient describes the within-cluster association. They refer to different comparisons.