Comparing M1–M4 on a Declining Concentration Profile

Why Compare the Methods on One Profile?

A below-quantification-limit (BQL) record is evidence about the concentration: the assay result indicates that the observation fell below the lower limit of quantification (LLOQ). The four methods considered here handle that evidence differently. M1 deletes the record, M2 replaces it by half the LLOQ, M3 uses the probability of falling below the LLOQ, and M4 replaces it by the LLOQ. On a declining concentration profile, those choices can change how much the terminal observations influence the estimated elimination rate.

The calculations below isolate that difference in a small hypothetical example. They are not a general ranking of methods. In particular, the fitted values depend on the model, residual error, observations, and chosen limit. A single profile can make the likelihood mechanics visible, but it cannot establish which approach will perform best in another analysis.

Model and Assumptions for the Comparison

Definition: Let the typical concentration at time \(t\) after a bolus dose be \(C(t)=20e^{-kt}\) mg/L, where \(k>0\) is the elimination rate constant in \(\mathrm{h}^{-1}\). For this illustrative comparison, each observation has an independent additive normal residual error with known standard deviation \(\sigma=1\) mg/L: \(Y(t)\mid C(t)\sim N(C(t),1)\). A BQL result at an LLOQ \(L\) means the observed value is known only to satisfy \(Y(t)<L\). The profile has a quantified value \(Y(1)=10\) mg/L and a BQL record at \(t=3\) h with \(L=2.5\) mg/L.

This is a deliberately simple, one-profile calculation: it has no interindividual random effects, no estimation of \(\sigma\), and no additional dose or compartment parameters. The dose and volume are represented by the fixed initial concentration of 20 mg/L. Conditional independence of the residual errors means that the contributions from the two times multiply in the likelihood. These assumptions let us compare the handling of one BQL record without attributing differences to changes in other model components.

Under the stated normal error model, a quantified value \(y\) at time \(t\) contributes the normal density \(f(y\mid C(t),1)\). The BQL record contributes the probability \(P(Y(3)<2.5\mid C(3))\) when it is treated as censored. For the methods that substitute a number, the substituted value is instead treated as if it were an observed concentration in the normal density. That change is precisely why substitution is not the same likelihood statement as censoring.

Likelihood Contributions Under M1–M4

Write \(C_1=20e^{-k}\) and \(C_3=20e^{-3k}\). The quantified observation contributes \(f(10\mid C_1,1)\) whenever it is retained. The methods differ in what they do with the record at three hours.

Method	Contribution from the BQL record	Meaning of that contribution
M1	No contribution	The BQL record is deleted; only the quantified observation informs the fit.
M2	\(f(1.25\mid C_3,1)\)	The BQL value is replaced by \(L/2=1.25\) mg/L and treated as quantified.
M3	\(\Phi(2.5-C_3)\)	The event \(Y(3)<2.5\) mg/L contributes its probability under the normal observation model.
M4	\(f(2.5\mid C_3,1)\)	The BQL value is replaced by \(L=2.5\) mg/L and treated as quantified.

Here \(\Phi\) is the standard normal cumulative distribution function. The argument is \(2.5-C_3\) because the residual standard deviation is 1 mg/L. More generally, the M3 contribution is \(\Phi((L-C_3)/\sigma)\). The three density-based methods use the same quantified contribution at one hour, but M1 has no terminal-record contribution and M2 and M4 use densities at invented replacement values.

Worked Example: Writing the Four Likelihoods

For the hypothetical data above, write the likelihood as a function of \(k\), omitting no method-specific factors. Let \(\phi(y-C)\) denote the normal density at \(y\) with mean \(C\) and standard deviation 1.

The one-hour quantified observation contributes \(\phi(10-C_1)\). Thus the four likelihoods are

\begin{aligned} L_{\mathrm{M1}}(k)&=\phi(10-C_1),\\ L_{\mathrm{M2}}(k)&=\phi(10-C_1)\phi(1.25-C_3),\\ L_{\mathrm{M3}}(k)&=\phi(10-C_1)\Phi(2.5-C_3),\\ L_{\mathrm{M4}}(k)&=\phi(10-C_1)\phi(2.5-C_3). \end{aligned}

For M3, the second factor is a probability, not a density: it is \(P(Y(3)<2.5\mid C_3)\). For M2 and M4, the second factor is a density because each replacement is treated as though it were an exact observed value. M1 has only one factor because the record is removed. These are distinct likelihood constructions, not four ways of writing the same likelihood.

A Score Derivation for M3

The M3 log-likelihood, up to constants that do not depend on \(k\), is

\ell_{\mathrm{M3}}(k)=-\frac{1}{2}(10-C_1)^2+\log\Phi(2.5-C_3).

Since \(dC_1/dk=-C_1\) and \(dC_3/dk=-3C_3\), differentiating the quantified term gives \(C_1(C_1-10)\). For the censored term, set \(z=2.5-C_3\). Then \(dz/dk=3C_3\), and the chain rule gives the derivative of \(\log\Phi(z)\) as \(3C_3\phi(z)/\Phi(z)\). Therefore the M3 score is

\frac{d\ell_{\mathrm{M3}}}{dk} =C_1(C_1-10)+3C_3\frac{\phi(2.5-C_3)}{\Phi(2.5-C_3)}.

The second term is positive: as \(k\) increases, the predicted terminal concentration falls and a BQL outcome becomes more probable. The quantified observation pulls in the opposite direction when \(C_1<10\), because increasing \(k\) then moves its predicted concentration farther below 10 mg/L. The M3 maximum balances these two score contributions.

Worked Fits and Terminal Predictions

Worked Example: M1 and M4 on the Same Constructed Profile

Under M1, only the quantified observation at one hour is used. The normal likelihood is maximized when its mean matches the observed value, so \(C_1=10\) mg/L. Substituting \(20e^{-k}=10\) gives

\hat{k}_{\mathrm{M1}}=\log(20/10)=\log 2=0.6931\ \mathrm{h}^{-1}, \qquad C_3=20e^{-3\log 2}=2.500\ \mathrm{mg/L}.

For M4, the substituted terminal value is 2.5 mg/L. At the same rate, both fitted means match the two values used by M4: \(C_1=10\) and \(C_3=2.5\). Thus the M4 likelihood is also maximized at \(k=\log 2\) for these constructed data, with a terminal prediction of 2.500 mg/L.

This equality is a feature of the example, not a general property of M1 and M4. The two constructed points happen to lie exactly on one exponential curve. With different quantified values, times, or substituted values, M4 will generally yield a different estimate from M1. M1's terminal prediction here is an extrapolation from the one quantified point and the fixed initial concentration; it is not direct use of the BQL record.

Worked Example: M2 With Half-LLOQ Substitution

M2 treats the one-hour value as 10 mg/L and the three-hour value as 1.25 mg/L. Up to constants, its log-likelihood is minus one half the sum of the two squared residuals. Differentiating with respect to \(k\) and setting the score to zero gives

$$ C_1(C_1-10)+3C_3(C_3-1.25)=0. $$

Solving this one-dimensional equation numerically gives \(\hat{k}_{\mathrm{M2}}\approx0.7515\ \mathrm{h}^{-1}\). Substitution into the profile gives

C_1=20e^{-0.7515}\approx9.433\ \mathrm{mg/L}, \qquad C_3=20e^{-3(0.7515)}\approx2.099\ \mathrm{mg/L}.

The fitted terminal mean is above the substituted value of 1.25 mg/L. The optimizer compromises between the one-hour observation and the substituted terminal point; it does not have enough freedom to match both exactly with one rate constant and the fixed initial concentration.

Worked Example: M3 Uses the BQL Event

For M3, solve the derived score equation

C_1(C_1-10)+3C_3\frac{\phi(2.5-C_3)}{\Phi(2.5-C_3)}=0.

The numerical maximum is approximately \(\hat{k}_{\mathrm{M3}}=0.7364\ \mathrm{h}^{-1}\), giving \(C_3=20e^{-3(0.7364)}\approx2.196\) mg/L. At this fitted mean, the probability of a BQL result is

P(Y(3)<2.5\mid C_3=2.196) =\Phi\!\left(\frac{2.5-2.196}{1}\right) =\Phi(0.304)\approx0.6194.

The predicted mean is below the LLOQ, but the observed BQL event does not assert that the mean itself must be below the LLOQ. It says the noisy observation was below the threshold. The probability calculation reflects both the mean and the residual standard deviation.

What the Comparison Shows—and Does Not Show

For these data, M1 and M4 give \(k\approx0.6931\ \mathrm{h}^{-1}\) and a terminal mean of 2.500 mg/L; M2 gives approximately \(0.7515\ \mathrm{h}^{-1}\) and 2.099 mg/L; M3 gives approximately \(0.7364\ \mathrm{h}^{-1}\) and 2.196 mg/L. M2 and M3 both yield a steeper decline than M1 in this example, but for different reasons: M2 treats 1.25 mg/L as an exact target, whereas M3 uses the probability of an observation below 2.5 mg/L.

Raw likelihood values or objective-function values from these four fits should not be compared as if they were computed from an identical set of observed data under a common likelihood. M1 omits a record, M2 and M4 replace it with a point value, and M3 uses an event probability. In particular, a density and a probability are different kinds of likelihood contributions. The calculations are useful for understanding how each rule changes the fit, not for selecting a method by simply choosing the largest raw likelihood across these altered data representations.

In a population NLME analysis, individual profiles contribute through a hierarchical model and the random-effects distribution; the effect of a BQL record is then conditional on that broader structure. The small calculation here is deliberately not a substitute for that analysis. It shows the local mechanism: handling the terminal record differently changes how the elimination rate is supported by the observed data.

Diagnostics, Pitfalls, and Limitations

For a practical fit, check the data status and the fitted behavior rather than assuming that a method's label guarantees correct implementation. Confirm the actual LLOQ and units at each record, verify that BQL status was not accidentally converted into an invented concentration for M3, and inspect fitted terminal predictions alongside the quantified data. A rapidly changing estimate when BQL records are included or excluded is a warning that the terminal behavior is sensitive to their handling.

Check the residual model and scale. The example uses additive normal error on the original concentration scale with known \(\sigma=1\) mg/L. Changing the residual model changes the BQL probability and density contributions. A probability computed with the wrong scale is not a valid M3 contribution.
Check numerical behavior. For M3, evaluate \(\Phi((L-C)/\sigma)\) and its log stably when the probability is very small. A numerical underflow or a boundary estimate can make the computed objective or score unreliable.
Do not equate a substituted value with the assay result. Under M2 and M4, the replacement is a modeling convention, not a measured concentration. Interpreting it as a true observation gives the terminal point more precise information than the BQL record supplied.
Do not infer general bias from one profile. The direction and size of a parameter shift depend on the full design, number and timing of BQL records, assay limits, residual variability, and model structure. The present example has one individual and a known residual standard deviation.

A further limitation is that the assumed initial concentration is fixed at 20 mg/L. If both the initial level and elimination rate were estimated, a single quantified point plus a terminal BQL event could provide limited information to distinguish them. Likewise, uncertainty in the residual error or a misspecified one-compartment profile would affect all fitted results. The comparisons therefore illustrate the consequences of the four handling rules under explicit assumptions; they do not remove the need to assess whether those assumptions suit the intended analysis.

Key takeaway: M1 discards the terminal record, M2 and M4 turn it into a point-valued pseudo-observation, and M3 retains its threshold information as a probability. The resulting terminal predictions are conditional on the model and error assumptions, and raw likelihood values across these different constructions are not directly comparable.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

Under the example's normal error model, what is the M3 contribution for a BQL observation with mean \(C_3\) and LLOQ 2.5 mg/L?
Why does the example give the same rate estimate under M1 and M4?
At the M3 fit, why is the predicted BQL probability not one even though the fitted mean is below the LLOQ?
Can the M2 and M3 raw likelihood values be compared directly to decide which handling method is better? Explain.

Solutions. (1) The contribution is \(\Phi((2.5-C_3)/1)\), the probability that the noisy observation falls below 2.5 mg/L. (2) The constructed values lie exactly on the same profile at \(k=\log 2\): \(C_1=10\) and \(C_3=2.5\); this coincidence is not general. (3) Residual error allows observations both above and below the mean, so a mean below the LLOQ gives a BQL probability between zero and one; here it is about 0.6194. (4) No. M2 uses a density at a substituted point, while M3 uses an event probability, so they do not represent the same observed-data likelihood.