Why Separate Variance Bias From Fixed-Effect Bias?
FO, FOCE, and Laplace approximate the likelihood in different ways, and their differences can affect different parameters by different amounts. A method might estimate a typical structural parameter accurately while estimating the spread of individual random effects too low. Conversely, the estimated variance might be close to its generating value while a fixed effect is biased. A comparison that reports only one parameter, or only an objective value, can miss this distinction.
The target here is the random-effect variance component. The central question is not whether one approximation must always overestimate or underestimate it: there is no universal direction. The direction and size of bias depend on the model, the sampling design, the amount of information per subject, and the sample size. Repeated simulation is useful because it separates a method's average tendency from the random variation in any one dataset.
Model, Assumptions, and What Each Approximation Changes
For a concrete scalar example, suppose subject \(i\) has observations \(y_{ij}\) at scheduled times \(t_{ij}\), with model
Here \(\theta\) is a fixed-effect parameter, \(\omega^2\) is the between-subject random-effect variance, and \(\sigma^2\) is the residual variance. For the discussion and examples, assume subjects are independent, residual errors are independent conditional on the random effects, the model is correctly specified, and the data-generating values are known for the simulation. The fitting procedure must also converge to the intended optimum; a failed or local optimization is not evidence of approximation bias.
FO uses a first-order Taylor approximation to the individual prediction, commonly around \(\eta_i=0\). FOCE linearizes around a subject-specific conditional estimate of the random effect. Laplace approximates the random-effect integral locally around its conditional mode, using the curvature there. Thus, these methods differ in how they represent the contribution of random effects to the marginal likelihood. FOCE and Laplace are not interchangeable: FOCE uses conditional linearization, while Laplace includes a local curvature contribution to its integral approximation.
Each method maximizes its own approximate likelihood, so it produces an estimator whose repeated-sampling behavior can differ from the estimator obtained by maximizing the exact marginal likelihood. In a simulation comparison, fitting all methods to the same generated datasets makes the contrast more informative. Still, a difference between methods is not automatically approximation bias alone: Monte Carlo error, optimization problems, and model misspecification must also be considered.
Why a First-Order Approximation Can Affect a Variance
A short expansion makes the mechanism visible. At a fixed observation time, write \(g\) and \(h\) for the first and second derivatives of \(f\) with respect to \(\eta\), evaluated at zero. A second-order local expansion is
For \(\eta\sim N(0,\omega^2)\), the centered normal moments are \(E(\eta)=0\), \(E(\eta^2)=\omega^2\), \(E(\eta^3)=0\), and \(E(\eta^4)=3\omega^4\). Therefore the expansion gives
For the variance result, the covariance between the linear and quadratic terms is zero because \(E(\eta^3)=0\). The variance of \(\eta^2\) is \(3\omega^4-\omega^4=2\omega^4\), so the quadratic term contributes \(\frac{1}{4}h^2(2\omega^4)=\frac{1}{2}h^2\omega^4\). A first-order approximation retains the \(g^2\omega^2\) contribution but omits this second-order term and the corresponding shift in the mean. Those omissions change the likelihood information about \(\omega^2\); they do not imply that FO must bias the variance downward in every model.
Worked Example: Curvature in a Log-Scale Random-Effect Model
Suppose the hypothetical structural response is \(f(\eta)=10\exp(\eta)\), with \(\eta\sim N(0,0.36)\). The fixed effect \(\theta=10\) is the response at \(\eta=0\); it is not the average response over the random-effect distribution.
For a lognormal random effect, \(E\{\exp(\eta)\}=\exp(\omega^2/2)\), so the marginal mean response is
The exact variance of \(\exp(\eta)\) is \(\exp(\omega^2)\{\exp(\omega^2)-1\}\). Substituting \(\omega^2=0.36\) gives
Consequently, the variance of the structural response is approximately \(10^2(0.6211)=62.11\) response-units squared. By contrast, the first-order expansion at zero has \(g=10\) and predicts random-effect variance \(g^2\omega^2=10^2(0.36)=36.00\) response-units squared. This calculation illustrates how curvature can make a first-order representation of the random-effect variation incomplete. It is not itself a fitted variance estimate, and it does not establish the direction of bias for every dataset or design.
Quantifying Bias in a Simulation
A simulation comparison should specify the generating model and design, generate repeated datasets, fit each approximation to each dataset, and summarize each parameter separately. The average estimate estimates the estimator's expectation under that particular setup. The simulation standard deviation describes the spread of estimates across simulated datasets; it is not the same quantity as bias.
Worked Example: Comparing Variance and Fixed-Effect Bias
Consider a hypothetical simulation using \(y_{ij}=\theta\exp(\eta_i)+\epsilon_{ij}\), with \(\theta=10\) response units, \(\omega^2=0.360\), and \(\sigma^2=0.250\). Suppose 1,000 datasets are generated, each with 100 subjects and two observations per subject, and the same datasets are fitted by FO, FOCE, and Laplace. The following is a hypothetical summary for illustrating the calculations, not a universal performance result.
| Method | Mean \(\hat{\omega}^2\) | SD of \(\hat{\omega}^2\) | Mean \(\hat{\theta}\) |
|---|---|---|---|
| FO | 0.245 | 0.110 | 10.06 |
| FOCE | 0.331 | 0.090 | 10.02 |
| Laplace | 0.354 | 0.085 | 10.01 |
For FO, the estimated variance bias is \(0.245-0.360=-0.115\) variance-units. Its relative bias is \((-0.115/0.360)\times100\%\approx-31.9\%\). For FOCE, the corresponding bias is \(0.331-0.360=-0.029\), or about \(-8.06\%\). For Laplace, it is \(0.354-0.360=-0.006\), or about \(-1.67\%\). In this particular hypothetical summary, FO has the largest negative variance bias and Laplace the smallest; the result should not be generalized beyond the stated model and design.
The estimated biases for \(\theta\), calculated separately, are \(10.06-10=0.06\) response units for FO, \(10.02-10=0.02\) for FOCE, and \(10.01-10=0.01\) for Laplace. These small fixed-effect biases do not erase the larger variance-component contrast. Nor does a modest fixed-effect bias certify that the variance estimate is adequate.
Simulation estimates are themselves uncertain. For FOCE, the Monte Carlo standard error of the mean variance estimate is approximately \(0.090/\sqrt{1000}=0.00285\). This describes simulation error in the estimated mean, not uncertainty for a variance estimate from one actual study. It helps assess whether small differences in simulated average estimates are distinguishable from Monte Carlo noise.
When Approximation Differences Can Disappear
Nonlinearity is central to many differences between approximations, but it is important not to assume the methods must disagree. In a linear-Gaussian model, the random-effect integral can be evaluated exactly, and the relevant local approximations can coincide with the exact likelihood. The next example demonstrates this boundary case.
Worked Example: A Linear-Gaussian Case With No Approximation Gap
Suppose each of five hypothetical subjects has one observation from \(y_i=\theta+\eta_i+\epsilon_i\), where \(\eta_i\sim N(0,\omega^2)\), \(\epsilon_i\sim N(0,0.040)\), and the residual variance is known. The observations are 9.2, 9.8, 10.1, 10.7, and 10.2 response units. Because the model is linear in \(\eta_i\) and both error terms are normal, the marginal distribution is exactly \(y_i\sim N(\theta,\omega^2+0.040)\).
The sample mean is \(\bar{y}=10.0\), and the sum of squared deviations is
For the normal likelihood with unknown mean and variance, the maximum-likelihood estimate of the marginal variance uses divisor \(n=5\), giving \(1.22/5=0.244\). Subtract the known residual variance to obtain the random-effect variance estimate: \(\hat{\omega}^2=0.244-0.040=0.204\) response-units squared. This is positive, so the nonnegative variance constraint does not alter it.
Here the linearization is exact, and the log integrand for the random effect is quadratic. FO, FOCE, and Laplace therefore recover the same marginal likelihood and the same estimates, including \(\hat{\theta}=10.0\) and \(\hat{\omega}^2=0.204\). This one dataset cannot establish bias: bias is an average over repeated datasets relative to a specified generating value. It does show why approximation differences should not be presumed when the model is linear-Gaussian.
Design Dependence, Comparisons, and Diagnostics
Sparse sampling can leave less within-subject information to distinguish random-effect variability from residual variability. In nonlinear models, that may make the variance estimate more sensitive to how each method approximates the random-effect contribution. More observations per subject can improve information, but it does not guarantee that all approximations will agree: the result still depends on the model's curvature, timing, error structure, and the range of observed responses.
Worked Example: A Hypothetical Sparse-versus-Rich Design Comparison
Continue with the same generating model and true variance \(\omega^2=0.360\). Suppose a second hypothetical simulation uses the same number of subjects and datasets but changes the observations per subject. Its mean variance estimates are summarized below.
| Observations per subject | FO | FOCE | Laplace |
|---|---|---|---|
| 1 (sparse) | 0.220 | 0.315 | 0.351 |
| 8 (richer) | 0.318 | 0.350 | 0.359 |
With one observation per subject, FO's estimated bias is \(0.220-0.360=-0.140\); with eight observations, it is \(0.318-0.360=-0.042\). The absolute FO bias is smaller by \(0.140-0.042=0.098\) variance-units in this hypothetical comparison. For Laplace, the corresponding biases are \(-0.009\) and \(-0.001\). The example illustrates that approximation performance can depend on design. It does not prove that adding observations will always reduce bias, nor that Laplace is always closest to the generating value.
A practical comparison should keep the following distinctions in view:
| Approach | What it approximates | What to examine for variance bias |
|---|---|---|
| FO | Prediction linearized around the population random-effect value | Whether omitted curvature materially changes the marginal variation represented by the fit |
| FOCE | Prediction linearized around subject-specific conditional random-effect estimates | Whether the conditional linearization improves the repeated-simulation variance estimates for this design |
| Laplace | Random-effect integral approximated using local curvature around a conditional mode | Whether the local curvature provides an adequate approximation for the fitted model and data |
| Exact integration or a sufficiently accurate numerical benchmark | The marginal random-effect integral without these local approximations | Whether an apparent approximation contrast persists relative to a more accurate likelihood calculation |
Useful warning signs include repeated estimates at the zero-variance boundary, large method-to-method differences across the same simulated datasets, unusually variable estimates, or convergence and Hessian warnings. These signals do not by themselves identify approximation bias. Check that all fits use the same model, parameter constraints, data, and convergence criteria. In simulation, report the mean estimate and bias along with the spread of estimates and, when useful, the Monte Carlo standard error. A numerical integration benchmark can help distinguish approximation effects from ordinary sampling variation, but it must itself be accurate enough for the comparison.
Do not compare raw objective values from different approximations as if they were directly interchangeable measures of variance bias. Each objective may approximate the likelihood differently. The relevant comparison for this topic is the behavior of the parameter estimates across repeated datasets, ideally also checked against a credible likelihood benchmark.
Common Mistakes and Limitations
- Inferring a universal direction. Seeing FO underestimate a variance in one simulation does not establish that FO always underestimates variance. The local expansion explains a possible mismatch, not a universal sign.
- Calling a single estimate “bias.” The difference between one fitted value and the generating value is a realization error. Bias is the expected estimation error, assessed through repeated sampling or a theoretical calculation.
- Confusing variance bias with fixed-effect bias. A fixed effect and a variance component have different targets and units. Report their biases separately; closeness of one does not imply closeness of the other.
- Mixing variance and standard-deviation scales. If \(\omega^2=0.360\), the random-effect standard deviation is \(\omega=0.600\). Bias in \(\omega\) is not numerically the same as bias in \(\omega^2\).
- Treating all differences as approximation effects. Monte Carlo error, boundary estimates, optimization failures, misspecification, and design can all affect observed results.
Simulation conclusions remain conditional on the chosen generating model and design. A study with different curvature, residual error, sample size, or observation schedule can show a different ranking. Furthermore, estimates near the nonnegative boundary can have irregular sampling behavior, so a mean bias alone may not describe the full distribution of estimates.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A simulation has true \(\omega^2=0.250\) and mean estimated variance \(0.190\). What are the estimated bias and relative bias?
- Why can FO, FOCE, and Laplace give identical estimates in a linear-Gaussian model?
- A method estimates \(\omega=0.500\), while the generating variance is \(\omega^2=0.360\). What is the estimated variance bias? Explain why comparing \(0.500\) directly with \(0.360\) is invalid.
- Across \(R=900\) simulated datasets, the standard deviation of estimated variances is \(0.120\). What is the Monte Carlo standard error of their mean?
Solutions. (1) Bias is \(0.190-0.250=-0.060\); relative bias is \((-0.060/0.250)\times100\%=-24.0\%\). (2) Linearization is exact for a linear predictor, and the random-effect integral has a quadratic Gaussian log integrand, which Laplace integrates exactly. (3) The estimate \(0.500\) is a standard deviation, so its implied variance estimate is \(0.500^2=0.250\). The variance bias is \(0.250-0.360=-0.110\) variance-units; the two original numbers are on different scales. (4) The Monte Carlo standard error is \(0.120/\sqrt{900}=0.00400\), in variance units.