A Simulation Demonstration Of Binary PQL Bias

Why Use a Simulation?

First-order penalized quasi-likelihood (PQL) replaces the nonlinear binary mixed-model calculation with a sequence of linear mixed-model approximations. This is useful computationally, but the approximation can affect estimates even when the data follow the model being fitted. A simulation gives a controlled way to examine that effect: generate data from a known model, fit PQL repeatedly, and compare the resulting estimates with the known parameter values.

The aim is not to show that one particular numerical pattern must occur in every setting. Rather, it is to make the simulation design and the interpretation of its output explicit. In particular, a mean estimate can be systematically displaced from its target (bias), while individual estimates also vary from one generated data set to another (sampling variability).

Model and Simulation Design

Model: For clusters $i=1,\ldots,K$ and observations $j=1,\ldots,m$, let $Y_{ij}$ be a binary outcome. The random-intercept logistic model is $$ Y_{ij}\mid u_i \sim \operatorname{Bernoulli}(p_{ij}),\qquad \operatorname{logit}(p_{ij})=\beta_0+u_i,\qquad u_i\sim N(0,\tau^2). $$ Conditional on the random effects, outcomes are independent; the random effects are independent across clusters and independent of any covariates. The parameters $\beta_0$ and $\tau^2$ are the conditional log-odds intercept and the random-intercept variance. In the simulation, data are generated from this model and fitted using first-order PQL.

A useful controlled experiment varies one or both design quantities that are expected to affect the approximation: cluster size $m$ and true random-intercept variance $\tau^2$. Hold the number of clusters $K$, fixed-effect value, data-generating distribution, fitting method, and other choices constant across conditions. Then changes in the results can be compared across the specified grid. In the illustration below, $K=200$, $\beta_0=-1$, $m$ is 5, 20, or 50, and $\tau^2$ is 0.25, 1, or 2.25. Each condition is repeated 1,000 times.

For each repetition, draw a new random effect for every cluster, then draw the binary outcomes conditional on those effects. Fit the same random-intercept PQL model and save its estimate $\hat{\tau}^2_r$, where $r$ indexes the repetition. Repeated data sets must be generated independently under the stated model. The mean of these fitted estimates describes the simulation’s average behavior under that particular design; it is not an estimate of a universal PQL bias.

The design holds $K$ fixed, so increasing $m$ also increases the total number of observations $Km$. The results therefore describe the combined change in cluster size and total observations under a fixed number of clusters. A simulation that holds the total number of observations fixed answers a different design question because the number of clusters must then change as $m$ changes.

Definition: For an estimator $\hat{\theta}$ of a fixed simulation target $\theta$, its bias is $E(\hat{\theta})-\theta$, where the expectation is over repeated data generation and fitting under the specified design. The Monte Carlo estimate of bias from $R$ repetitions is $\overline{\hat{\theta}}-\theta$, where $\overline{\hat{\theta}}=R^{-1}\sum_{r=1}^R\hat{\theta}_r$. The sample standard deviation of the $R$ estimates describes their simulated sampling variability.

What Bias and Variability Measure

Bias and sampling variability are different features of an estimator. An estimate can be highly variable but centered near its target, or it can be relatively stable while being systematically displaced. For a finite-variance estimator, mean squared error (MSE) combines both:

\begin{aligned} \operatorname{MSE}(\hat{\theta}) &= E\left[(\hat{\theta}-\theta)^2\right]\\ &= E\left[\{(\hat{\theta}-E(\hat{\theta}))+(E(\hat{\theta})-\theta)\}^2\right]\\ &= \operatorname{Var}(\hat{\theta})+\{E(\hat{\theta})-\theta\}^2. \end{aligned}

The final equality follows because $E[\hat{\theta}-E(\hat{\theta})]=0$, so the cross-product has expectation zero. Thus a small standard deviation across repeated fits does not, by itself, establish good estimation: the squared bias may still be large. This decomposition motivates reporting both the mean’s displacement from the known target and the spread of the repeated estimates.

In a simulation, $\overline{\hat{\theta}}-\theta$ is only an approximation to the estimator’s theoretical bias because it is computed from a finite number of repetitions. Its Monte Carlo standard error is the standard deviation of the simulated estimates divided by $\sqrt{R}$. This quantifies numerical uncertainty in the simulated mean; it is not the standard error of an estimate from one study.

Worked Example: Calculating Bias From Repeated Fits

Suppose a hypothetical simulation uses $K=200$, $m=20$, and true $\tau^2=1$, with $R=1{,}000$ independent repetitions. The fitted PQL variance-component estimates have mean 0.72 and sample standard deviation 0.22.

The estimated bias is $$ \overline{\hat{\tau}^2}-\tau^2=0.72-1.00=-0.28. $$ The negative sign means the mean fitted variance is below the generating variance. Relative bias, defined here as estimated bias divided by the true value, is $-0.28/1.00=-0.28$, or $-28\%$.

The standard deviation 0.22 describes how much the estimates vary across generated data sets. It is not the bias and is not the Monte Carlo standard error of the mean. The latter is $$ SE_{\mathrm{MC}}=\frac{0.22}{\sqrt{1{,}000}}=0.00696. $$ Thus the simulation’s mean is estimated with much less Monte Carlo uncertainty than the 0.28 displacement from the target. The calculation is evidence of bias in this hypothetical design, not a claim about every PQL analysis.

Reading a Design Grid

The following table gives hypothetical simulation summaries for the design just specified. The values are deliberately illustrative: they demonstrate how to calculate and compare summaries, and are not results from a reported empirical study. Each row summarizes 1,000 converged PQL fits. The standard deviation column is the sample standard deviation of the fitted variance estimates across repetitions.

True $\tau^2$	Cluster size $m$	Mean $\hat{\tau}^2$	Estimated bias	SD of $\hat{\tau}^2$
0.25	5	0.18	-0.07	0.15
0.25	20	0.22	-0.03	0.09
0.25	50	0.23	-0.02	0.06
1.00	5	0.55	-0.45	0.31
1.00	20	0.72	-0.28	0.22
1.00	50	0.80	-0.20	0.17
2.25	5	0.90	-1.35	0.55
2.25	20	1.35	-0.90	0.48
2.25	50	1.60	-0.65	0.40

Worked Example: Comparing Bias Across Conditions

For true $\tau^2=2.25$ and $m=5$, the table gives mean $\hat{\tau}^2=0.90$. The estimated bias is $$ 0.90-2.25=-1.35. $$ Dividing by the target gives relative bias $-1.35/2.25=-0.60$, or $-60\%$. Under this hypothetical condition, the mean fitted variance is 60% below the generating variance.

For the same true variance and $m=50$, the mean is 1.60, so the estimated bias is $1.60-2.25=-0.65$, or $-0.65/2.25=-0.289$, approximately $-28.9\%$ relative bias. The larger cluster-size condition is less biased in this illustration, but remains downward biased. This comparison supports a conclusion about these simulation conditions only; it does not guarantee that increasing cluster size will remove PQL bias in a different model or design.

Across the displayed conditions, the mean variance-component estimate is below its target. The discrepancy is larger for the larger true random-intercept variances, and the illustrated bias decreases as cluster size increases. These patterns are consistent with the concern that first-order linearization may not adequately represent the nonlinear binary response, especially when heterogeneity is substantial. They should be treated as patterns to investigate, not as mathematical laws inferred from a single grid.

Worked Example: Distinguishing Bias From Monte Carlo Noise

Return to the hypothetical condition with $m=20$ and true $\tau^2=1$. The simulated mean is 0.72, the across-repetition standard deviation is 0.22, and $R=1{,}000$, so $SE_{\mathrm{MC}}=0.00696$. An approximate 95% Monte Carlo interval for the simulation’s mean estimate is $$ 0.72\pm 1.96(0.00696)=[0.706,\;0.734]. $$

This interval quantifies Monte Carlo uncertainty about the mean PQL estimate under the specified simulation design. It is not a confidence interval for the true $\tau^2$, and it does not describe the range of estimates expected from a single data set. Since the interval lies well below 1, the observed downward displacement is not plausibly explained by ordinary Monte Carlo fluctuation in the simulated mean alone. The simulation still estimates behavior only for its chosen model, fitting method, and design.

Assumptions, Diagnostics, and Comparisons

The interpretation depends on both the data-generating process and the fitting process. For this demonstration, the assumed process is the stated logistic random-intercept model: binary outcomes, normal independent random intercepts with the specified variance, and conditional independence within clusters given those intercepts. The fitting procedure must apply first-order PQL to the intended random-intercept model in every repetition. If generation or fitting differs from these specifications, the summaries no longer isolate the behavior being studied.

A simulation should also check its own implementation. Useful checks include verifying that generated random effects have approximately the target variance, that outcomes are drawn using the stated conditional probabilities, and that the stored estimates correspond to the intended model and parameter. Record convergence and boundary estimates rather than silently dropping unsuccessful fits. If some repetitions fail, report how many and clarify whether summaries are conditional on successful fits; otherwise, the reported mean may represent a selected subset of the simulated data sets.

Monte Carlo precision can be assessed by increasing the number of repetitions or by calculating the Monte Carlo standard error for each reported mean. A wide Monte Carlo error relative to the estimated bias means the simulation is too noisy to distinguish the pattern reliably. Conversely, a small Monte Carlo error makes the simulated mean precise for its design but does not correct model misspecification or make that design representative of all applications.

The main alternative is to compare PQL with an estimator based on the integrated likelihood, such as likelihood-based numerical integration. The latter targets the likelihood after accounting for the random-effects distribution, while PQL uses a first-order working approximation; computational demands and finite-sample behavior can differ. Such comparisons require fitting methods to be applied to the same generated data sets. The simulation’s purpose is not to declare one method universally preferable, but to reveal how estimates behave under a specified model and design.

The main limitation is scope: a random-intercept model with one intercept, one cluster-size pattern per condition, and a fixed-effect intercept cannot establish what happens with covariates, random slopes, misspecified random-effect distributions, or other link functions. Nor does the simulation identify a threshold at which PQL becomes reliable. Those are separate questions requiring separately designed comparisons.

Common Pitfall

Common mistake: “The standard deviation of the estimates is small, so PQL has little bias.” The standard deviation measures repetition-to-repetition spread, not displacement from the target. A collection of tightly grouped estimates can still be centered far below the true variance. Report the mean-minus-target bias and the across-repetition standard deviation separately.

Check Your Understanding

Answer each question using the definitions and hypothetical summaries above. Brief solutions follow the questions.

In the model, what does $\tau^2$ represent, and what does $m$ represent?
If the true variance is 1.00 and the simulated mean estimate is 0.80, calculate the estimated bias and relative bias.
Why is the across-repetition standard deviation not the Monte Carlo standard error of the simulated mean?
What additional design change occurs when $m$ increases while $K$ stays fixed?
What should be reported if some PQL fits fail to converge?

Solutions. (1) $\tau^2$ is the variance of the cluster-specific random intercepts; $m$ is the number of binary observations per cluster in that condition. (2) Bias is $0.80-1.00=-0.20$; relative bias is $-0.20/1.00=-20\%$. (3) The across-repetition standard deviation describes the spread of individual fitted estimates, whereas the Monte Carlo standard error of their mean is that standard deviation divided by $\sqrt{R}$. (4) The total number of observations $Km$ increases. (5) Report the number or proportion of failed fits and state whether summaries exclude them; clarify that estimates summarized only among successful fits may be selected.