Approximation Bias With Skewed Conditional Effects

Why Conditional Skewness Matters

A marginal likelihood integrates over the possible random effects for each individual. Laplace approximation replaces the integrand near its conditional mode with a Gaussian-shaped curve, using the mode’s height and local curvature. This can work well when the important part of the integrand is close to symmetric and quadratic. It can be less accurate when the conditional random-effect distribution is strongly skewed: a local description centered on the mode may not represent the longer tail or the changing shape farther from the mode.

The skewness at issue is conditional on the individual’s observations and model parameters. Even if the population random effects have a normal distribution, the conditional distribution can be asymmetric when the individual likelihood depends nonlinearly on those effects. This tutorial focuses on how that asymmetry can affect the integrated likelihood, not on a general comparison of all approximation methods.

Definition: For individual $i$, let $\eta_i$ be the random-effects vector, $y_i$ the observations, and $\theta$ the model parameters. Define the individual log-integrand $$ h_i(\eta_i;\theta)=\log p(y_i\mid\eta_i;\theta)+\log p(\eta_i\mid\theta). $$ The individual marginal likelihood is $L_i(\theta)=\int \exp\{h_i(\eta_i;\theta)\}\,d\eta_i$. The conditional random-effect distribution is proportional to this integrand as a function of $\eta_i$. It is skewed when its shape is not symmetric about its center; the mode, median, and mean need not coincide.

In the common normal random-effects model, $p(\eta_i\mid\theta)$ is normal, but the product of that density and a nonlinear individual likelihood need not be. A nonlinear concentration-time model, for example, can make a modest change in a clearance random effect have a very different impact on the likelihood on either side of its best-fitting value. The resulting conditional distribution may have a long tail even though the population distribution is symmetric.

The Local Gaussian Approximation

Suppose $h_i$ has a unique interior mode $\hat\eta_i$, is sufficiently smooth nearby, and has a negative-definite Hessian there. Write $J_i=-\nabla^2 h_i(\hat\eta_i;\theta)$, which is positive definite under this condition. A second-order Taylor approximation around the mode gives a quadratic log-integrand, and integrating its exponential gives the Laplace approximation.

h_i(\eta_i;\theta)\ \approx\ h_i(\hat\eta_i;\theta) -\frac{1}{2}(\eta_i-\hat\eta_i)^\mathsf{T}J_i(\eta_i-\hat\eta_i)

L_i(\theta)\ \approx\ L_{i,\mathrm{Laplace}}(\theta) =\exp\{h_i(\hat\eta_i;\theta)\}(2\pi)^{q/2}|J_i|^{-1/2},

where $q$ is the dimension of $\eta_i$. This result is most persuasive when the integral is dominated by a neighborhood of the mode, the mode is well located, and the log-integrand is close to quadratic throughout the neighborhood that contributes materially to the integral. A strongly skewed conditional distribution warns that this last condition may fail: the quadratic approximation is symmetric about the mode, while the actual integrand may put appreciably different mass on either side.

The local approximation uses the mode and second derivatives, not the entire shape of the integrand. A short justification comes from writing $\eta_i=\hat\eta_i+u$ and expanding further:

h_i(\hat\eta_i+u;\theta) =h_i(\hat\eta_i;\theta)-\frac{1}{2}u^\mathsf{T}J_i u +\text{terms of order three and higher in }u.

The quadratic term produces the Gaussian integral in the Laplace formula. The omitted terms describe departures from that local shape. In a symmetric Gaussian neighborhood, a cubic term alone has zero first-order integral because it is odd in $u$; that does not mean skewness has no effect. The full integrand is exponential, and higher-order terms and interactions among terms can change its total mass. Thus the local formula has no general guarantee of overestimating or underestimating the integral when the conditional distribution is skewed.

Worked Calculations With a Skewed Integrand

The following gamma-shaped integrands provide transparent numerical illustrations of local approximation error. They are mathematical examples of skewed integration kernels, not claims that a gamma distribution is the random-effects distribution in a Gaussian NLME model. In an NLME model with Gaussian population random effects, nonlinear individual likelihoods can nonetheless produce skewed conditional integrands.

Worked Example: One Skewed Integral

Consider the normalized integrand $g(\eta)=\eta^3e^{-\eta}/\Gamma(4)$ for $\eta>0$, and $g(\eta)=0$ otherwise. Its exact integral is 1, since it is a gamma density with shape 4 and rate 1. Find its one-dimensional Laplace approximation.

The log-integrand is $h(\eta)=3\log\eta-\eta-\log\Gamma(4)$. Differentiating gives $h'(\eta)=3/\eta-1$, so the mode is $\hat\eta=3$. The second derivative is $h''(\eta)=-3/\eta^2$, hence $J=-h''(3)=1/3$. Since $\Gamma(4)=3!=6$, $$ \exp\{h(3)\}=\frac{3^3e^{-3}}{6}=\frac{4.5}{e^3}\approx 0.2240. $$ The Laplace approximation is therefore $$ L_{\mathrm{Laplace}}=\exp\{h(3)\}\sqrt{\frac{2\pi}{J}} =\frac{4.5}{e^3}\sqrt{6\pi} \approx 0.9727. $$

The exact integral is 1, so this approximation is about 2.73% below the exact value. The density has skewness $2/\sqrt{4}=1$, and its right tail is longer than its left side. This example shows that using the mode and curvature can miss some total mass; it does not establish that skewness always causes underestimation.

Worked Example: Accumulation Across Individuals

Suppose 20 hypothetical individuals have the same integrand shape as in the previous calculation, but with an additional multiplicative factor $C=0.08$ for each individual. The exact individual marginal likelihood is $0.08$, while the Laplace value is $0.08(0.9727)\approx0.07782$. Assuming conditional independence across individuals, compare the products of the individual contributions.

The ratio of the approximate to exact data-set likelihood is $$ \left(\frac{0.08(0.9727)}{0.08}\right)^{20} =(0.9727)^{20}\approx 0.5749. $$ Thus, in this constructed example, the approximate product is about 57.49% of the exact product. On the $-2\log L$ scale, the difference (approximate minus exact) is $$ -2\log\{(0.9727)^{20}\}=-40\log(0.9727)\approx 1.107. $$ The small per-individual discrepancy has accumulated across the 20 contributions. This calculation concerns likelihood values, not a probability that one model is true.

Worked Example: Shape Can Change a Likelihood Comparison

Consider two hypothetical parameter settings whose normalized conditional kernels are gamma densities with shape $4$ and shape $9$, each with rate 1. Both exact integrals are 1. Compare their Laplace approximations. For shape $k$, the mode is $k-1$, and the Laplace value is $$ r_k=\frac{(k-1)^{k-1}e^{-(k-1)}}{\Gamma(k)} \sqrt{2\pi(k-1)}. $$

For $k=4$, the calculation from the first example gives $r_4\approx0.9727$. For $k=9$, the mode is 8, $\Gamma(9)=8!=40320$, and $$ r_9=\frac{8^8e^{-8}}{40320}\sqrt{16\pi}\approx0.9896. $$ The shape-4 kernel has skewness $2/\sqrt{4}=1$, while the shape-9 kernel has skewness $2/\sqrt{9}\approx0.6667$. The exact marginal likelihoods are tied, but the local approximations differ. For 20 repeated contributions, the difference between their approximate $-2\log L$ values is $$ 20\{-2\log(r_4)+2\log(r_9)\}\approx20(0.05536-0.02082)=0.6908. $$ The approximation therefore makes the shape-9 setting look better on this objective, although their exact integrals are equal by construction. This is an illustration of how shape-dependent approximation error can affect comparisons, not a general prediction about which setting an NLME analysis will favor.

What to Check in an NLME Analysis

In practice, the exact individual integrals are usually unavailable, so a discrepancy cannot be diagnosed from skewness alone. Useful warning signs are indirect: the conditional distribution appears markedly asymmetric in plots or numerical summaries; the conditional mode lies far from a conditional mean estimate; or modest changes in integration settings or approximation method materially change parameter estimates or the objective. These observations are prompts to investigate, not proof that a particular estimate is biased.

The earlier material on FOCE and Laplace approximation distinguishes their constructions: FOCE uses a first-order expansion around conditional estimates, while Laplace uses local curvature at the conditional mode. Both rely on local information, so neither should be assumed immune to an integrand that is strongly non-Gaussian. A more direct numerical integration can serve as a comparison when feasible, though its accuracy must also be checked. The choice among approximations depends on the model, data, computational cost, and how stable the relevant conclusions are across reasonable calculations.

Several limitations matter when interpreting an apparent approximation bias:

Skewness alone does not determine the sign or size of the integral error. Tail weight and the full shape also matter.
A skewed conditional distribution can reflect genuine information in the observations and is not, by itself, evidence that the population model is misspecified.
Approximation error in the likelihood is not the same as sampling bias in an estimator. It can contribute to differences in estimates, but the direction and magnitude depend on the data and model.
Strong skewness is a warning about a local Gaussian representation, not a quantitative error bound.

Common pitfall: “The conditional distribution is skewed, so Laplace must overestimate the marginal likelihood.” This conclusion is not justified. Laplace substitutes a symmetric local Gaussian for the integrand, but the omitted shape can add or remove mass relative to that Gaussian. The direction must be checked for the particular integrand, preferably by comparison with a reliable integration calculation when possible.

Key takeaway: Laplace approximation retains the integrand’s height and local curvature at its mode, but not its full conditional shape. Strong conditional skewness can therefore produce meaningful marginal-likelihood error, especially when modest individual discrepancies accumulate or alter a comparison between parameter settings.

Check Your Understanding

Answer each question, then compare with the solutions below.

Can a normal population distribution for $\eta_i$ produce a skewed conditional distribution given $y_i$? Explain briefly.
A one-dimensional log-integrand has an interior mode and $h''(\hat\eta)=-4$. What is the local Laplace width factor $\sqrt{2\pi/J}$?
If a Laplace approximation is 0.98 times the exact individual likelihood for each of 10 independent individuals, what is the approximate-to-exact likelihood ratio for the full data set?
Does a skewed conditional distribution imply that Laplace overestimates the marginal likelihood? Why or why not?

Solutions. (1) Yes. A nonlinear individual likelihood multiplied by a normal population density need not be symmetric as a function of $\eta_i$. (2) Here $J=-h''(\hat\eta)=4$, so the factor is $\sqrt{2\pi/4}=\sqrt{\pi/2}\approx1.253$. (3) The ratio is $0.98^{10}\approx0.8171$. (4) No. Skewness signals that a symmetric local approximation may miss mass, but it does not determine the direction of the error; that depends on the full integrand.

Approximation Bias With Skewed Conditional Effects

What You'll Learn