M3: Using a Censored Likelihood

Why M3 Uses a Probability for a BQL Record

A quantified concentration gives an observed value, so its likelihood contribution describes how much density the model assigns at that value. A below-quantification-limit (BQL) record is different: the assay establishes only that the measured result fell below the lower limit of quantification (LLOQ). It does not establish the unobserved numerical value. M3 represents that information with the probability of the observed event, rather than by inventing a concentration.

This distinction matters because a BQL record can still inform the model. If the model predicts a high concentration with little residual variability, a result below the LLOQ may be unlikely. If it predicts a low concentration, that same BQL event may be quite plausible. M3 lets the likelihood reflect this evidence, conditional on the pharmacokinetic parameters and individual effects.

Definition: Let \(Y_{ij}\) be the modelled observation for individual \(i\) at record \(j\), with conditional density \(f(y\mid\theta,\eta_i)\). Here \(\theta\) denotes population parameters and \(\eta_i\) the individual random effects. For a quantified observation \(y_{ij}\), the contribution is \(f(y_{ij}\mid\theta,\eta_i)\). For a BQL observation with LLOQ \(L_{ij}\), the contribution is \(P(Y_{ij}<L_{ij}\mid\theta,\eta_i)\), equivalently the conditional cumulative distribution function evaluated at the LLOQ.

For a continuous observation model, the probability of a value being exactly equal to the LLOQ is zero. Thus using \(Y_{ij}<L_{ij}\) or \(Y_{ij}\leq L_{ij}\) gives the same likelihood contribution under that model. The threshold must be expressed on the same scale as \(Y_{ij}\). If an observation model is specified on a transformed scale, the event and threshold must be transformed consistently.

From the Observation Model to the Likelihood

The BQL contribution follows directly from the probability of the reported event. If the conditional observation model has density \(f\), then the probability that a measurement is below the LLOQ is the integral of that density over all values below the threshold:

P(Y_{ij}<L_{ij}\mid\theta,\eta_i) = \int_{-\infty}^{L_{ij}} f(y\mid\theta,\eta_i)\,dy.

The lower integration limit is the lower end of the observation model's support; for a model supported only on positive values, the integral begins at zero instead. A quantified value contributes a density, not a probability of observing that exact point. A BQL record contributes a probability, because the reported datum is an interval event.

Suppose an individual's records are conditionally independent given that individual's random effects and the model parameters. If \(Q_i\) is the set of quantified records and \(B_i\) is the set of BQL records, the conditional likelihood is the product of the corresponding contributions:

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

Conditional independence here concerns the observation errors after conditioning on the individual effects and model predictions. It does not say that all records from an individual are unconditionally independent: they share \(\eta_i\). In a population model, the individual likelihood is therefore obtained by integrating the joint product over the distribution of random effects:

L_i(\theta) = \int L_i(\theta,\eta_i)\,g(\eta_i\mid\Omega)\,d\eta_i, \qquad L(\theta)=\prod_i L_i(\theta).

Here \(g(\eta_i\mid\Omega)\) is the random-effects density, with \(\Omega\) denoting its distribution parameters. This order is essential: form the individual's conditional likelihood from all of that individual's records, then integrate over the shared random effects. The overall likelihood multiplies the marginal individual likelihoods when individuals are independent.

Worked Example: A Normal Residual Model and One BQL Result

Consider a hypothetical measurement with conditional model \(Y\mid\theta,\eta_i\sim N(6.0,1.5^2)\), where concentrations are in mg/L. The predicted mean is 6.0 mg/L, the residual standard deviation is 1.5 mg/L, and the LLOQ is 4.0 mg/L. Calculate the BQL contribution. Then calculate the density contribution if a separate quantified observation under the same conditional model is 7.0 mg/L.

For the BQL event, standardize the LLOQ using the normal model:

z=\frac{4.0-6.0}{1.5}=-1.333\quad\text{(rounded)}.

The probability below this threshold is \(\Phi(-1.333)\approx0.09127\). Thus the BQL record contributes approximately 0.09127, a dimensionless probability. For the quantified value \(y=7.0\) mg/L, standardization gives \(z=(7.0-6.0)/1.5=0.6667\). The normal density on the concentration scale is:

f(7.0) = \frac{1}{1.5\sqrt{2\pi}} \exp\!\left(-\frac{(7.0-6.0)^2}{2(1.5)^2}\right) \approx 0.21297\ {\rm (mg/L)}^{-1}.

If these were two conditionally independent records for the same individual at the stated parameter values, their conditional likelihood contribution would be \(0.21297\times0.09127\approx0.01944\ {\rm (mg/L)}^{-1}\). The density has inverse-concentration units, while the BQL probability is dimensionless. The example assumes a normal additive residual model on the original concentration scale; it is the specified observation model, not a claim that every assay result must follow that model.

Worked Likelihoods With Mixed Records

The numerical contribution for a BQL record is not a fixed fraction of the LLOQ. It depends on the model prediction, residual-error parameters, and threshold. Likewise, the density contribution for a quantified record depends on the full observation model, not merely on whether the value is above the LLOQ.

Worked Example: Combining One Quantified and One BQL Record

For a hypothetical individual at fixed \(\theta\) and \(\eta_i\), suppose a quantified observation has conditional density \(0.12099\ {\rm (mg/L)}^{-1}\), while a second record is BQL and the model assigns it probability 0.6915 of falling below its LLOQ. Assume the two observation errors are conditionally independent given these model values.

The conditional likelihood contribution is the density multiplied by the event probability:

L_i(\theta,\eta_i) = 0.12099\times0.6915 = 0.083664585 \approx0.08366\ {\rm (mg/L)}^{-1}.

The corresponding log-likelihood contribution is \(\log(0.083664585)\approx-2.481\). This is a likelihood contribution, not a probability for the exact quantified value: its units are inverse mg/L because one density and one dimensionless probability were multiplied. In a fitted population model, this conditional value is not generally the final individual contribution; integration over \(\eta_i\) is still required.

The following two-node calculation illustrates why that integration must follow the joint product. It is a deliberately small numerical approximation to the random-effects integral. Each node has weight 0.5, and the given node-specific values are conditional contributions for one quantified and one BQL record.

Worked Example: Integrating the Joint Contribution Over Random Effects

At random-effects node 1, suppose the quantified density is \(0.19947\ {\rm (mg/L)}^{-1}\) and the BQL probability is 0.15866. At node 2, suppose they are \(0.06476\ {\rm (mg/L)}^{-1}\) and 0.69146. Use equal weights of 0.5.

First multiply the two record contributions within each node, then average the node-specific products:

\begin{aligned} L_i(\theta) &\approx 0.5(0.19947\times0.15866) + 0.5(0.06476\times0.69146)\\ &= 0.5(0.0316479102)+0.5(0.0447789496)\\ &=0.0382134299 \approx0.03821\ {\rm (mg/L)}^{-1}. \end{aligned}

The log contribution is approximately \(\log(0.0382134299)=-3.265\). If instead the density and probability were each averaged first, the approximate marginal density would be \(0.5(0.19947+0.06476)=0.132115\ {\rm (mg/L)}^{-1}\), and the approximate marginal probability would be \(0.5(0.15866+0.69146)=0.42506\). Multiplying those averages gives about \(0.05616\ {\rm (mg/L)}^{-1}\), not 0.03821. That calculation loses the fact that the records from this individual share the same random effect. The joint conditional product must be integrated as a whole.

Assumptions, Diagnostics, and Practical Limits

M3 is appropriate when the BQL status means that the underlying assay measurement was below a known threshold and the observation model describes that measurement process. Its contribution is valid only relative to the chosen residual-error distribution and the model prediction. The likelihood also relies on using the correct threshold for each record and on the stated conditional-independence structure. If assay processing or reporting rules create additional censoring or selection, those mechanisms must be represented rather than silently ignored.

Several checks help identify implementation or model problems:

Check the event and scale. Confirm that the CDF is evaluated at the LLOQ on the observation model's scale. A threshold on the original scale cannot be inserted unchanged into a model defined on a transformed scale.
Check contributions by record type. Quantified observations should use densities; BQL observations should use cumulative probabilities. A BQL record should not be assigned a fabricated concentration and then treated as exact.
Check numerical tails. A very small below-LLOQ probability can cause underflow if probabilities are multiplied directly. Implementations commonly work with log probabilities and stable log-CDF calculations. A zero produced by numerical underflow is not the same as a theoretically zero probability.
Check the random-effects calculation. The marginal likelihood should integrate the within-individual product. Verify quadrature settings or other integration methods, especially when the likelihood changes materially with numerical settings.
Check model fit against the BQL pattern. Systematic runs of BQL records or patterns of BQL status that the fitted predictions cannot reproduce may indicate an inadequate structural or residual-error model. A numerical optimizer can converge while the observation model remains poorly specified.

A useful derivation also clarifies why a BQL probability changes with prediction. Under a normal additive-error model with conditional mean \(\mu_{ij}\) and residual standard deviation \(\sigma\), \(Y_{ij}\mid\theta,\eta_i\sim N(\mu_{ij},\sigma^2)\). Standardizing the event gives:

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

Thus, holding the LLOQ and residual standard deviation fixed, a higher predicted mean lowers the standardized threshold and usually lowers the probability of a BQL observation. This conclusion follows from the monotonic increase of the standard normal CDF: as \(\mu_{ij}\) increases, \((L_{ij}-\mu_{ij})/\sigma\) decreases. The formula is specific to this normal model; another residual-error distribution requires its own CDF.

How M3 Differs From Substitution Approaches

M3 uses the information actually reported: a quantified value contributes its density, while a BQL value contributes the probability of the censoring event. In contrast, deletion discards that record's likelihood information, and half-LLOQ substitution treats an assigned numerical value as though it had been measured. Those approaches can change parameter estimates because they do not represent the same observation process as censoring. M3 is not automatically correct simply because it is a likelihood: it can still be biased or misleading if the structural model, residual model, threshold, or censoring assumptions are wrong.

M3 also has practical costs. Evaluating CDFs and integrating over random effects can make estimation more demanding, and strong dependence on the residual model can matter when many observations are BQL. When the fitted model predicts almost no probability for observed BQL events, inspect the model and implementation rather than forcing an arbitrary replacement value. When the likelihood is numerically unstable or parameters are weakly identified, estimation diagnostics and sensitivity checks remain important.

Common mistake: “A BQL result is close enough to the LLOQ, so I can enter half the LLOQ as an exact observation in the M3 likelihood.” This is not M3. It replaces the observed event \(Y<L\) with an invented point value and a density at that point. M3 uses \(P(Y<L\mid\theta,\eta_i)\), retaining the uncertainty about the unobserved measurement.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A quantified measurement is recorded as 5.2 mg/L. What kind of likelihood contribution does it receive, and what kind does a BQL record with LLOQ 3.0 mg/L receive?
Under a normal additive-error model, the predicted mean is 5.0 mg/L, the residual standard deviation is 1.0 mg/L, and the LLOQ is 4.0 mg/L. Write the expression for the BQL probability.
Why must the conditional contributions for several records from one individual be multiplied before integrating over that individual's random effects?
Suppose a computed BQL probability is numerically returned as zero. Give one diagnostic question to ask before concluding that the model assigns exactly zero probability.

Solutions. (1) The quantified value receives the conditional density evaluated at 5.2 mg/L; the BQL record receives \(P(Y<3.0\text{ mg/L}\mid\theta,\eta_i)\), the conditional CDF at 3.0 mg/L. (2) The standardized threshold is \((4.0-5.0)/1.0=-1\), so the probability is \(\Phi(-1)\approx0.1587\). (3) The records share the same individual random effects. The population likelihood therefore averages their joint conditional likelihood over that shared effect; multiplying separately averaged contributions would generally give a different result. (4) Check whether the probability underflowed numerically, and whether the CDF or log-CDF was evaluated on the correct scale with the correct LLOQ. A numerical zero alone does not establish a theoretical zero probability.