Marginalizing Censored Likelihoods Over Random Effects

From Conditional BQL Probabilities to Population Likelihoods

A below-quantification (BQL) record does not specify a measured concentration; it specifies that the unobserved concentration fell below the assay's lower limit of quantification (LLOQ). Conditional on a subject's random effects, the model assigns a probability to that event. Population inference requires more: because the random effects are unobserved, the subject's likelihood contribution averages the conditional likelihood over their population distribution.

The order of these operations matters. For one subject, first combine the conditional contributions from all that subject's records, then integrate that joint contribution over the random effects. In general, multiplying separately marginalized record contributions is not equivalent. This is the central computational point when censored records and interindividual variability appear together.

Definition: Let $\theta$ denote population parameters, $\eta_i$ the random-effects vector for subject $i$, and $g(\eta_i;\theta)$ its population density. Let $L_i(\theta,\eta_i)$ be the likelihood of all subject $i$'s records conditional on $\eta_i$. The subject's marginal likelihood contribution is $$ L_i^{\mathrm{marg}}(\theta)=\int L_i(\theta,\eta_i)g(\eta_i;\theta)\,d\eta_i. $$ It is the conditional likelihood averaged over the random-effects distribution, not the conditional likelihood evaluated at the average random effect.

Conditional Contributions for Quantified and BQL Records

For a concrete model, suppose the latent concentration $Y_{ij}$ for subject $i$ at record $j$, conditional on $\eta_i$, is lognormally distributed: $\log Y_{ij}\mid\eta_i\sim N(\log F_{ij}(\theta,\eta_i),\sigma^2)$. Here $F_{ij}$ is the model-predicted median concentration and $\sigma$ is the residual standard deviation on the log scale. A quantified value $y_{ij}$ contributes its density on the concentration scale. A BQL record with LLOQ $L_{ij}$ contributes the probability $P(Y_{ij}<L_{ij}\mid\eta_i)$.

f(y_{ij}\mid\eta_i)=\frac{1}{y_{ij}\sigma\sqrt{2\pi}} \exp\!\left[-\frac{(\log y_{ij}-\log F_{ij}(\theta,\eta_i))^2}{2\sigma^2}\right], \qquad y_{ij}\ge L_{ij},

P(Y_{ij}<L_{ij}\mid\eta_i) =\Phi\!\left(\frac{\log L_{ij}-\log F_{ij}(\theta,\eta_i)}{\sigma}\right).

The second expression is a conditional probability, not a density evaluated at the LLOQ. For example, it increases when the predicted concentration falls relative to the LLOQ. For a different residual-error distribution, use that model's corresponding conditional cumulative probability. The censored-likelihood principles below do not depend on choosing the lognormal model.

Suppose, for now, that records are conditionally independent given $\eta_i$. If $Q_i$ is the set of quantified records and $B_i$ is the set of BQL records, the conditional likelihood is

L_i(\theta,\eta_i)= \prod_{j\in Q_i} f(y_{ij}\mid\eta_i) \prod_{j\in B_i} P(Y_{ij}<L_{ij}\mid\eta_i).

This factorization assumes that the model accounts for dependence through the random effects and that the residuals are conditionally independent. If the residual model specifies dependence between records, use the appropriate joint conditional likelihood instead of this product. If subjects are independent and have independent random-effects draws, the population likelihood is $\prod_i L_i^{\mathrm{marg}}(\theta)$.

Worked Example: Marginalizing One BQL Record

Hypothetical illustration. At one time point, let $L=1.00$ mg/L, $F(\eta)=1.40\exp(\eta)$ mg/L, $\eta\sim N(0,0.20^2)$, and $\sigma=0.25$. Calculate the population probability of a BQL record using three-point Gauss–Hermite quadrature.

For a standard normal expectation, this three-point rule uses nodes $0,\pm\sqrt{3}$ with weights $2/3,1/6,1/6$, respectively. Because $\eta=0.20Z$ for $Z\sim N(0,1)$, the corresponding random-effect nodes are $0,\pm0.3464$. At each node, calculate the conditional BQL probability $p(\eta)=\Phi((\log 1.00-\log(1.40)-\eta)/0.25)$.

Since $\log(1.40)=0.3365$, the conditional probabilities are approximately $p(-0.3464)=\Phi(0.0397)=0.5158$, $p(0)=\Phi(-1.346)=0.0891$, and $p(0.3464)=\Phi(-2.732)=0.00315$. Thus the quadrature estimate is

L^{\mathrm{marg}}_{\mathrm{BQL}} \approx \frac{1}{6}(0.5158)+\frac{2}{3}(0.0891)+\frac{1}{6}(0.00315) =0.1459.

The estimated population BQL probability is about $0.146$, or 14.6%. It is an average of conditional probabilities over the random-effects distribution. It is not the BQL probability for a hypothetical subject whose random effect is fixed at zero; that conditional probability is $0.0891$.

Why Marginalize the Joint Subject Contribution?

All records from one subject share the same $\eta_i$. Consequently, conditional record contributions are evaluated at the same random-effect vector before being combined. The marginal likelihood then averages that joint contribution. Multiplying separately averaged contributions would effectively allow each record to have a different, independently drawn random effect, which is a different model.

Worked Example: Two Records Share One Random Effect

Hypothetical illustration. For a subject with one BQL record and one quantified record, suppose a numerical integration rule has three random-effect nodes with weights $0.25,0.50,0.25$. At those nodes, the conditional BQL probabilities are $0.80,0.30,0.05$, and the quantified-record densities are $0.10,0.50,0.20$ per mg/L. Assume conditional independence of the two records given the shared random effect.

At each node, multiply the two conditional contributions first. The products are $0.80(0.10)=0.080$, $0.30(0.50)=0.150$, and $0.05(0.20)=0.010$ per mg/L. The marginal joint contribution is therefore

0.25(0.080)+0.50(0.150)+0.25(0.010) =0.020+0.075+0.0025=0.0975\ \mathrm{(mg/L)^{-1}}.

By contrast, the separately marginalized BQL probability would be $0.25(0.80)+0.50(0.30)+0.25(0.05)=0.3625$. The separately marginalized density would be $0.25(0.10)+0.50(0.50)+0.25(0.20)=0.325$ per mg/L. Their product is $0.3625(0.325)=0.1178$ per mg/L, not $0.0975$. That product does not represent the specified shared-random-effect model.

Numerical Integration: Choices and Diagnostics

Except in special models, the integral over $\eta_i$ has no closed-form solution. It must be approximated repeatedly as population parameters change during likelihood optimization. The approximation must represent the integrand $L_i(\theta,\eta_i)g(\eta_i;\theta)$, which can vary substantially across subjects and parameter values.

Justification: The random effects are latent, so the joint contribution of subject $i$'s records and $\eta_i$ is $L_i(\theta,\eta_i)g(\eta_i;\theta)$. Integrating this joint contribution over all possible $\eta_i$ gives the likelihood of the observed records alone: $$ P(\text{records}_i\mid\theta) =\int P(\text{records}_i\mid\eta_i,\theta)g(\eta_i;\theta)\,d\eta_i =\int L_i(\theta,\eta_i)g(\eta_i;\theta)\,d\eta_i. $$ Thus the integration is over the population distribution of random effects, while the censored probabilities remain conditional contributions inside the integral.

Common approaches make different trade-offs. Gaussian quadrature evaluates the integrand at weighted nodes. Adaptive quadrature centers and scales the nodes to better cover the subject-specific region contributing most to the integral. These methods can be accurate in low dimensions, but the number of nodes grows rapidly as the random-effects dimension increases. A Laplace approximation uses a local quadratic approximation around the conditional mode of the integrand; it can be efficient, but its accuracy depends on how well a local approximation represents the full integral. Monte Carlo integration averages integrand values at random draws, while quasi-Monte Carlo uses structured low-discrepancy points. Simulation-based methods can handle more dimensions, but their approximation error can make the likelihood noisy, particularly during optimization.

Approach	Strength	Diagnostic concern
Gaussian or adaptive quadrature	Deterministic and often effective in low dimensions	Check stability as node count or adaptation settings change
Laplace approximation	Can reduce computation when the integrand is well approximated near its mode	Check for sensitivity to poor local approximation or a poorly located mode
Monte Carlo or quasi-Monte Carlo	Can be useful as random-effects dimension grows	Check simulation error and stability across draws, seeds, or point counts

Log-scale computation is important regardless of the integration method. A product of many small densities and BQL probabilities may underflow if formed directly. Compute conditional log-likelihoods as sums, and combine weighted contributions using a stable log-sum-exp calculation when appropriate. The likelihood must still be integrated on its proper scale: taking an average of log-likelihoods is not the same as taking the log of the average likelihood.

Worked Example: Assessing Monte Carlo Integration Error

Hypothetical illustration. Suppose four independent draws from the random-effects distribution give conditional BQL probabilities $0.90,0.40,0.10,0.20$. The Monte Carlo estimate of the marginal probability is their average:

\widehat{p}_{\mathrm{BQL}}=\frac{0.90+0.40+0.10+0.20}{4}=0.400.

The sample variance of the four values is $[ (0.90-0.40)^2+(0.40-0.40)^2+(0.10-0.40)^2+(0.20-0.40)^2 ]/(4-1)=0.1267$. The estimated standard error of the Monte Carlo average is therefore $\sqrt{0.1267/4}=0.178$. With only four draws, this uncertainty is large relative to the estimate. The calculation illustrates how to assess simulation error; it is not evidence that $0.400$ is a precise estimate. Increasing the number of draws and checking stability are necessary for reliable integration.

Failure Modes and Practical Interpretation

A numerical approximation can appear to converge while still misrepresenting the likelihood. Compare estimates under more nodes or draws, and inspect changes in the objective and parameter estimates—not only whether the software reports successful optimization. Large differences across integration settings suggest that the likelihood approximation may be influencing inference. Extreme conditional probabilities also deserve care: values near zero can dominate a log-likelihood, and their logarithms must be evaluated stably.

Heavy censoring can make subject contributions especially sensitive to how random effects are integrated: many records contribute probabilities rather than observed concentrations. That does not change the marginalization formula, but it can concentrate or reshape the integrand in ways that challenge a coarse quadrature rule or local approximation. Numerical integration diagnostics address approximation quality; they do not establish that the model parameters are identifiable or that the assumed residual and random-effects distributions are appropriate.

A useful implementation check is to compare the integrated result with a simpler calculation when one is available. In the lognormal example above, $\log Y=\log(1.40)+\eta+\epsilon$, where $\eta\sim N(0,0.20^2)$ and $\epsilon\sim N(0,0.25^2)$ independently. Hence $\log Y$ has variance $0.20^2+0.25^2=0.1025$, and the exact marginal BQL probability is $\Phi[-0.3365/\sqrt{0.1025}]$, approximately $0.147$. This agrees closely with the three-point quadrature estimate $0.1459$, providing a useful check in this particular example. Many nonlinear PK models do not permit this analytic shortcut.

Common Mistake

A common error is to calculate a BQL probability at $\eta_i=0$ and treat it as the subject's population likelihood contribution. That is a conditional probability for a typical random-effect value, not an average over interindividual variability. Another error is to integrate each record separately and multiply the results, even though the records share the same random effects. The correction is to construct the full conditional subject likelihood using the appropriate quantified densities and BQL probabilities, then integrate that joint contribution against the random-effects density.

Key takeaway: For each subject, BQL probabilities enter inside the conditional likelihood alongside quantified-record densities. Combine all records conditional on the shared random effects, then integrate the joint contribution over their population distribution.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a lognormal residual model, what conditional contribution does a BQL record with LLOQ $L$ make?
Why is the product of separately marginalized record contributions generally not the subject's marginal likelihood?
A one-node approximation evaluates a BQL conditional probability of $0.12$ at $\eta=0$. Does that value necessarily equal the population BQL probability?
Name one diagnostic for quadrature and one for Monte Carlo integration.

Solutions. (1) It contributes $\Phi((\log L-\log F(\theta,\eta_i))/\sigma)$, conditional on the subject's random effects. (2) The records share the same random-effects vector, so their conditional contributions must be multiplied before integration; separate integration treats them as if their random effects were independently drawn. (3) No. It is the conditional probability at $\eta=0$; the population probability averages conditional probabilities over the random-effects distribution. (4) For quadrature, increase the node count and check whether the likelihood and estimates stabilize. For Monte Carlo, increase the draw count and assess simulation error or stability across independent runs.