Separating Approximation Bias From Sampling Variability

Why Separate Bias From Variability?

A PQL estimate from one data set can differ from the target parameter for at least two reasons. The estimate may be systematically displaced because PQL approximates the likelihood-related calculation, or it may differ simply because another sample from the same population would produce a different estimate. Those explanations are not interchangeable: more data may reduce sampling variability without removing an approximation bias.

The distinction is most directly studied by repeatedly generating data from a specified GLMM and applying the same PQL fitting procedure to every generated data set. Since the generating parameters are known in this exercise, each estimate can be compared with its target. The repeated estimates reveal the behavior of the procedure under that particular model and design; they do not establish that the same numerical bias will occur in every application.

These quantities are defined relative to the full procedure being studied: the data-generating model, cluster design, PQL fitting method, and rule for handling failed fits. For example, a binary random-intercept model could generate independent clusters, with

Here \(\beta_1\) is a fixed-effect coefficient and \(\tau^2\) is a random-intercept variance. In a simulation, the generating values of \(\beta_1\) and \(\tau^2\), as well as the covariates and cluster sizes, must be specified before the estimates are examined. The PQL fit is then applied consistently to data generated under those choices. This setup assumes that the simulated observations follow the stated model, the clusters are generated according to the stated independence structure, and the fitting and recording rules are applied consistently across replicates.

Repeated-Sampling Summaries and the MSE Decomposition

Suppose \(R\) data sets are generated independently under the same model and design. Let \(\hat{\theta}_r\) be the PQL estimate from replicate \(r\). The simulation estimates the bias with the average estimation error and describes sampling variability with the spread of the replicate estimates:

The sign of the estimated bias matters. A negative value means the average PQL estimate is below the generating value; a positive value means it is above. The standard deviation \(s_{\hat{\theta}}\) is not another estimate of bias. It summarizes how much the estimates vary from one generated data set to another.

A central identity relates both sources of error to mean squared error (MSE). For any estimator with finite second moment,

To justify it, add and subtract \(E(\hat{\theta}_{\text{PQL}})\) inside the squared error. The cross-product has expectation zero because the centered estimator has expectation zero:

The first term is variability; the second is squared systematic displacement. For a finite simulation, the average squared error around the known target estimates MSE. It is also possible to decompose that finite-simulation quantity exactly using the simulation mean and variance calculated with divisor \(R\). The usual reported sample standard deviation uses divisor \(R-1\), so it does not plug into that finite-simulation identity without a small adjustment.

Variance Components Need Their Own Summaries

A random-effect variance is also a parameter, so the same definitions apply with \(\theta=\tau^2\). However, a variance component is constrained to be nonnegative. Its sampling distribution can be asymmetric, may pile up near zero, and need not be well described by a normal approximation. Consequently, report the replicate distribution or useful quantiles as well as an average and standard deviation when the shape is important. Do not infer that a symmetric interval is appropriate merely because the mean and standard deviation have been calculated.

Design Comparisons and Simulation Precision

Bias and variability can respond differently when the design changes. For instance, a larger number of clusters may reduce the dispersion of estimates, but that alone does not show that a PQL approximation bias has disappeared. The following hypothetical simulation summaries illustrate how to compare both quantities on the same parameter scale.

Simulation summaries are themselves estimates based on a finite number \(R\) of replicates. If the replicate estimates have sample standard deviation \(s_{\hat{\theta}}\), the Monte Carlo standard error of their average is approximately \(s_{\hat{\theta}}/\sqrt{R}\), under independent replicate generation. For example, with \(R=1{,}000\) and \(s_{\hat{\theta}}=0.10\), this standard error is \(0.10/\sqrt{1{,}000}=0.00316\). It describes simulation noise in the estimated mean, not the sampling standard deviation of a single PQL estimate. Increasing \(R\) makes the simulation summary more precise; it does not change the estimator’s underlying bias or variability for the chosen design.

Assumptions, Diagnostics, and Limits

A useful repeated-sampling assessment requires a clearly specified target. The generating model should reflect the features being investigated, including the outcome distribution, link, random-effects distribution, covariates, cluster sizes, and number of clusters. Replicates must be generated independently under that same setup. If the simulation changes several design features at once, observed differences cannot be attributed to one feature alone.

The fitting rule also matters. Record the number of attempted fits, successful fits, convergence warnings, and boundary estimates. If failed fits are excluded, the reported summaries describe successful fits rather than the full procedure; if failures are assigned values or handled by another rule, that rule must be stated. A high failure rate can make a simple average of the successful estimates misleading. It can also indicate that the estimator’s behavior is not adequately summarized by a mean and standard deviation.

Useful checks include plotting or tabulating the replicate estimates, checking for skewness and boundary accumulation, and comparing the average with the generating value. For fixed effects, report bias in coefficient units; for a logit slope, those are log-odds units per one-unit increase in the covariate. For a variance component, report squared units and consider relative bias only when the generating variance is positive. If estimates are skewed, quantiles may communicate their distribution more faithfully than a standard deviation alone.

Repeated-sampling simulation is not the only way to describe uncertainty. A standard error from an individual fitted model targets sampling variability under that model and its approximations; it does not, by itself, measure the difference between the PQL estimator’s repeated-sampling mean and the true parameter. A simulation can estimate both quantities when the generating truth is known. In an applied data analysis, that truth is not known, so a simulation study is informative only to the extent that its assumed model and design are relevant. Neither a small simulation standard error nor a narrow sampling distribution establishes that approximation bias is negligible.

The earlier material on PQL bias explains why systematic displacement can arise, including for fixed effects and random-effect variances. The practical lesson here is narrower: do not use the word “bias” for ordinary scatter, and do not use a standard error as though it measured approximation bias. Keep the target, average error, and replicate spread distinct.