Why Integrate Over Individual Random Effects?
A population PK model describes both typical structural parameters and variation among individuals. For subject \(i\), the predicted concentration may depend on population parameters \(\theta\) and an individual random-effect vector \(\eta_i\). The observed concentrations are therefore informative about \(\eta_i\), but \(\eta_i\) is not usually a directly observed parameter. The population likelihood must account for its distribution rather than treating an estimate of it as known.
This distinction matters because the same subject’s random effect influences all of that subject’s concentrations. Integrating over \(\eta_i\) allows the data to contribute information about population parameters while accounting for uncertainty about the subject-specific effect. Substituting an estimated effect into the conditional likelihood does not, in general, produce the population likelihood.
From Conditional to Population Likelihood
Let \(y_i\) denote all concentration observations for subject \(i\), including their sampling times and any other relevant design information. Let \(g_i(y_i\mid\eta_i;\theta,\psi)\) be their conditional data density, where \(\theta\) denotes structural population parameters and \(\psi\) denotes residual-model parameters. Suppose the random effects have density \(p(\eta_i;\Omega)\), with mean zero and covariance matrix \(\Omega\). Then the marginal density for subject \(i\) is
The integral is over the full vector of random effects assigned to that subject. If subjects are independent conditional on the population parameters and each subject has an independently drawn random-effect vector, the population likelihood is the product of these subject-level contributions:
The subject-level conditional density must include the complete residual model for that subject. For example, if residuals are conditionally independent given \(\eta_i\), it can be written as a product over that subject’s observations. If the residual model includes within-subject correlation, the appropriate joint conditional density replaces that product. Either way, a single \(\eta_i\) is shared across all observations from subject \(i\), unless the model explicitly gives different effects to different occasions or records.
Why the Integral Is a Marginal Likelihood
The justification is the law of total probability. The data can be observed under many possible values of the unobserved random effect. Each conditional data density is weighted by how probable that effect is under the population distribution. Integrating sums those weighted contributions across the possible values of \(\eta_i\). Independence across subjects then gives the product form above.
In general, the integral does not simplify to evaluating the conditional likelihood at a single “best” effect. A short Gaussian derivation shows the difference. Suppose one log-concentration observation \(z_i\) follows \(z_i=\mu+\eta_i+\epsilon_i\), with \(\eta_i\sim N(0,\omega^2)\), \(\epsilon_i\sim N(0,\sigma^2)\), and independent \(\eta_i\) and \(\epsilon_i\). Then \(z_i-\mu\) is the sum of two independent normal variables, so
This is the integral of the conditional normal density over the normal random-effect distribution. The variance \(\omega^2+\sigma^2\) reflects uncertainty from both sources. A conditional likelihood evaluated at an estimated \(\eta_i\) instead treats that effect as fixed at one value and does not produce this marginal density.
Worked Example: An Exact Marginal Density for Repeated Log Concentrations
Consider one hypothetical subject with two log-concentration observations \(z_{i1}=3.1\) and \(z_{i2}=2.9\). Use the model \(z_{ij}=\mu+\eta_i+\epsilon_{ij}\), where \(\mu=3.0\), \(\eta_i\sim N(0,0.3^2)\), and the two conditionally independent residuals satisfy \(\epsilon_{ij}\sim N(0,0.2^2)\). The random effect is shared by both observations.
Transform the two observations to an orthogonal mean coordinate and contrast coordinate: \(q_1=(z_{i1}+z_{i2})/\sqrt{2}=4.2426\) and \(q_2=(z_{i1}-z_{i2})/\sqrt{2}=0.1414\). The mean coordinate has mean \(\sqrt{2}\mu=4.2426\) and variance \(\sigma^2+2\omega^2=0.04+2(0.09)=0.22\). The contrast has mean zero and variance \(\sigma^2=0.04\), because the shared random effect cancels in the difference.
The marginal density is the product of the two independent transformed-coordinate densities:
Substitute the values. Since \(q_1-\sqrt{2}\mu=0\) and \(q_2^2=0.02\),
For comparison, the conditional likelihood at \(\eta_i=0\) is \((2\pi(0.04))^{-1}\exp\{-0.02/(2(0.04))\}\approx 3.099\). The conditional maximum-likelihood estimate of \(\eta_i\) is also zero here, because the subject’s average log concentration equals \(\mu\). That plug-in value is not the marginal likelihood: the marginal calculation averages over the possible shared effects, rather than conditioning on the one value that maximizes the conditional likelihood.
Numerical Integration for Nonlinear PK Models
For nonlinear concentration predictions, integrating the conditional likelihood over a random-effect distribution often has no closed-form solution. A numerical method approximates the integral with weighted evaluations of the conditional likelihood. In one dimension, a three-point Gauss–Hermite rule for a standard normal random variable uses nodes \(0\) and \(\pm\sqrt{3}\), with respective weights \(2/3\) and \(1/6\). For \(\eta_i\sim N(0,\omega^2)\), the corresponding effect values are \(0\) and \(\pm\sqrt{3}\omega\).
These weights and nodes approximate an integral under the assumed random-effect distribution; they are not three observed types of subject. The same principle extends to multiple dimensions, though the number of evaluations can grow rapidly as the number of random effects increases.
Worked Example: A Three-Point Quadrature Approximation
Suppose a hypothetical concentration observation is \(C_i=24\) mg/L. Conditional on a scalar random effect, let \(C_i\mid\eta_i\sim N(20\exp(\eta_i),4^2)\), with \(\eta_i\sim N(0,0.25^2)\). The conditional density, in \((\mathrm{mg/L})^{-1}\), is \(g(24\mid\eta_i)=\{4\sqrt{2\pi}\}^{-1}\exp\{-[24-20\exp(\eta_i)]^2/32\}\).
The three effect values are \(0\) and \(\pm\sqrt{3}(0.25)=\pm0.4330\). At these values the conditional means are approximately \(20.00\), \(12.97\), and \(30.84\) mg/L. The corresponding conditional densities are approximately \(0.06049\), \(0.00223\), and \(0.02314\) \((\mathrm{mg/L})^{-1}\). Weighting and summing gives
The conditional density at \(\eta_i=0\) is \(0.06049\;(\mathrm{mg/L})^{-1}\), not the marginal density. The quadrature result is an approximation: its accuracy depends on how well the chosen nodes and weights represent the integral, particularly if the conditional likelihood is sharply peaked or strongly asymmetric.
Integrate Within Subjects, Then Multiply
The order of operations encodes the model. Each subject has their own random-effect draw, shared across that subject’s records. Therefore, first integrate the joint conditional likelihood for all of one subject’s observations over that subject’s effect; then multiply the resulting contributions across subjects. Integrating separately for every observation would instead allow a new independent random effect at each record, changing the model and generally overstating within-subject independence.
Worked Example: Combining Two Subject Contributions
Suppose two hypothetical subjects each have one log-concentration observation, \(z_1=3.5\) and \(z_2=2.8\). Use the Gaussian random-intercept model with \(\mu=3.0\), \(\omega=0.3\), and \(\sigma=0.2\). Each marginal observation has variance \(\omega^2+\sigma^2=0.13\). Thus,
For subject 1, the density is \(1.106\exp\{-0.5^2/0.26\} =1.106\exp(-0.9615)\approx0.4230\). For subject 2, it is \(1.106\exp\{-(-0.2)^2/0.26\} =1.106\exp(-0.1538)\approx0.9487\). The population likelihood for these two observations is therefore
This product follows from independent subject-specific effects and independent subjects conditional on the population parameters. If the model instead included a shared group-level effect, the subject contributions would not be independent in this way unless that shared effect were also handled appropriately.
Marginal Likelihood Versus Plug-In Effects
An individual effect estimate, such as a conditional mode or an empirical Bayes estimate, summarizes the distribution of that effect given the subject’s data and the population parameters. It is useful for describing subject-specific predictions, but replacing the integral by the conditional likelihood at that estimate is a different calculation. The estimate is data-dependent, its uncertainty is not represented by a single point, and its value can change as the population parameters change.
| Approach | Calculation | What it represents |
|---|---|---|
| Marginal likelihood | Integrate the conditional density over the random-effect distribution | Probability density of the subject’s data under the population model |
| Plug-in conditional likelihood | Evaluate the conditional density at an estimated individual effect | Conditional fit at one estimated effect; not generally the population likelihood |
| Numerical marginalization | Approximate the integral using quadrature, Laplace methods, or simulation | An approximation to the same marginal likelihood, subject to numerical error |
Laplace approximations, Gaussian quadrature, and Monte Carlo integration are computational strategies for evaluating or approximating the marginal integral. They differ in their approximation error and computational cost; none changes the definition of the target likelihood. A plug-in effect calculation, by contrast, targets a conditional quantity rather than approximating the integral itself.
Assumptions, Diagnostics, and Limitations
The marginal likelihood is only as appropriate as the model being integrated. Its interpretation relies on a specified random-effect distribution, a correctly defined conditional concentration model, and a residual model that represents the conditional data adequately. The usual product over subjects also relies on the stated independence structure. The distribution need not be normal in principle, but the density used in the integral must match the chosen model.
Numerical integration can fail even when the model is correctly specified. Useful checks include increasing quadrature accuracy and confirming that parameter estimates and likelihood values are stable; comparing results from different numerical strategies when feasible; checking that densities and weights are evaluated on the intended scale; and inspecting whether likelihood contributions are dominated by a narrow region not well covered by the integration points. In Monte Carlo integration, repeated calculations with different simulation draws can reveal excessive simulation variability.
- Integrating each record separately: This gives each observation a separate random-effect draw. Use one shared subject-level effect for all records governed by that effect.
- Using an estimated effect as if it were observed: A conditional mode is not a substitute for integration. Retain the integral when constructing the population likelihood.
- Ignoring the residual model: The integrand is the full conditional density for the subject’s data, not just the structural prediction or a residual sum of squares unless that quantity is part of the correct density.
- Trusting a numerical approximation without checking it: A computed value is not evidence of adequate integration. Check convergence or stability as numerical accuracy is increased.
Marginal likelihood also has practical limitations. High-dimensional random-effect integrals can be expensive, and a convenient approximation may be inaccurate for a sharply peaked, skewed, or multimodal integrand. Furthermore, a well-computed likelihood does not establish that the structural model, random-effect distribution, or residual assumptions are scientifically adequate. Numerical accuracy and model adequacy are separate questions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In a model with one random effect shared by four concentrations from a subject, what data enter that subject’s conditional density before integration?
- For one log-concentration observation with \(\omega^2=0.09\) and \(\sigma^2=0.04\), what is its marginal variance under the Gaussian random-intercept model?
- Why is the product of two subject-level marginal likelihoods appropriate when subjects are independent, but not automatically appropriate when they share a group-level random effect?
- What does a three-point quadrature rule approximate in a marginal-likelihood calculation?
- State one diagnostic for detecting inadequate numerical integration and one limitation that remains even if integration is numerically accurate.
Solutions. (1) The joint conditional density of all four observations, with the same subject-level effect in each observation’s model, enters the integral. (2) The marginal variance is \(\omega^2+\sigma^2=0.09+0.04=0.13\). (3) Independent subjects with independent effects have a joint density that factors into subject contributions; a shared group effect induces dependence unless it is accounted for in the joint model. (4) It approximates the integral of the conditional data density with respect to the random-effect distribution, using weighted evaluations at selected effect values. (5) One can increase quadrature accuracy and check whether estimates and likelihood values stabilize. Even with accurate integration, the assumed structural, random-effect, or residual model may still be inadequate.