Worked M3 Analysis of a Sparse PK Profile

Why work through a sparse profile?

A BQL record does not report a concentration value; it reports that the concentration fell below the assay’s lower limit of quantification (LLOQ). In a sparse pharmacokinetic (PK) profile, a small number of quantified observations may therefore carry much of the numerical information, while one or more BQL records still constrain the plausible concentration-time curve. The M3 method uses both types of information without pretending a BQL record is an exact concentration.

This worked analysis makes the likelihood calculation explicit, then compares an M3 fit with a simple M4 fit. The calculations use a deliberately small hypothetical data set so each step can be checked. The model is simplified for arithmetic; a population PK analysis generally estimates structural and random-effects parameters using all subjects’ records.

Definition: Let \(Y_i\) be a concentration measurement with model-predicted mean \(\mu_i\), residual standard deviation \(\sigma\), and LLOQ \(L_i\). Under a normal residual model, a quantified observation \(y_i\) contributes its normal density evaluated at \(y_i\). A BQL observation contributes the probability \(P(Y_i<L_i)=\Phi((L_i-\mu_i)/\sigma)\), where \(\Phi\) is the standard normal cumulative distribution function. M3 combines these contributions in one likelihood.

Likelihood contributions under a stated error model

For the worked calculations, assume that, conditional on the model parameters, each observation has a normal distribution with mean \(\mu_i\) and common additive residual standard deviation \(\sigma\). The LLOQ is known, the BQL indicator is recorded correctly, and observations are conditionally independent given the model. The normal model is on the concentration scale, so \(\sigma\) and \(L_i\) have concentration units. These assumptions are not automatic: a different residual distribution, an assay-specific LLOQ, or a dependence structure requires the corresponding likelihood formulation.

For a quantified observation, the likelihood contribution is a density, not a probability of observing one exact real-valued number. Its units are inverse concentration. For a BQL observation, the contribution is a probability and has no units. Writing the log-likelihood is convenient because independent contributions add.

f(y_i\mid\mu_i,\sigma)=\frac{1}{\sqrt{2\pi}\sigma} \exp\!\left[-\frac{(y_i-\mu_i)^2}{2\sigma^2}\right], \qquad P(Y_i<L_i\mid\mu_i,\sigma) =\Phi\!\left(\frac{L_i-\mu_i}{\sigma}\right).

The strict or non-strict inequality at the LLOQ makes no difference under a continuous model, because the probability of exactly \(L_i\) is zero. For a profile with quantified index set \(Q\) and BQL index set \(B\), the conditional log-likelihood is

\ell = \sum_{i\in Q} \left[ -\log(\sqrt{2\pi}\sigma) -\frac{(y_i-\mu_i)^2}{2\sigma^2} \right] + \sum_{i\in B} \log\Phi\!\left(\frac{L_i-\mu_i}{\sigma}\right).

In a population model, a subject’s random effect is shared across that subject’s observations. The subject likelihood is therefore obtained by multiplying the conditional contributions across records and then integrating over the random-effects distribution. It is generally incorrect to integrate each record separately and multiply those marginal record contributions; doing so discards the dependence induced by the shared random effect.

L_j = \int \left[ \prod_{i\in Q_j} f(y_{ij}\mid\eta_j) \prod_{i\in B_j} P(Y_{ij}<L_{ij}\mid\eta_j) \right] g(\eta_j)\,d\eta_j.

Here \(g(\eta_j)\) is the random-effects density, and the conditional mean and possibly the residual scale depend on \(\eta_j\). The examples below have no random effects, so their likelihood calculations are conditional and do not require this integration.

Worked Example: Calculating One Quantified and One BQL Contribution

Suppose a hypothetical sample has LLOQ \(L=1.0\) mg/L. At one time, the model predicts \(\mu=2.0\) mg/L and the quantified measurement is \(y=2.2\) mg/L. At another time, the predicted mean is \(\mu=1.6\) mg/L and the record is BQL. Let the additive residual standard deviation be \(\sigma=0.4\) mg/L.

Quantified record. The residual is \(2.2-2.0=0.2\) mg/L. Its density contribution is

f(2.2\mid 2.0,0.4) = \frac{1}{\sqrt{2\pi}(0.4)} \exp\!\left[-\frac{(0.2)^2}{2(0.4)^2}\right] = 0.9974\exp(-0.125) \approx 0.8802\;(\mathrm{mg/L})^{-1}.

BQL record. Its standardized LLOQ is \(z=(1.0-1.6)/0.4=-1.5\), so the contribution is

P(Y<1.0\mid\mu=1.6,\sigma=0.4) = \Phi(-1.5) \approx 0.06681, \qquad \log\{\Phi(-1.5)\}\approx -2.706.

The BQL contribution says that, under this predicted mean and residual variability, about 6.7% of measurements would fall below 1.0 mg/L. It does not assign the unobserved concentration a value of zero, half the LLOQ, or the LLOQ.

Fitting a one-parameter M3 model

To make the fit reproducible by hand, consider a hypothetical sparse profile whose structural predictions have a known shape and an unknown concentration scale \(\theta\). At four sampling times, let the predictions be \(\mu_i=\theta f_i\), where the fixed factors are \(f=(1,0.5,0.25,0.125)\). Use additive normal error with \(\sigma=1.0\) concentration unit and a common LLOQ \(L=1.0\) concentration unit. The records are \(y_1=4.0\) and \(y_2=2.0\), both quantified; records 3 and 4 are BQL. For this example, assume \(\theta>0\).

The model shape and error scale are treated as fixed; only \(\theta\) is estimated. This is a teaching calculation, not a claim that a four-record profile can reliably estimate a full population PK model. The quantified observations contribute normal densities, and the two BQL records contribute probabilities.

Worked Example: Maximizing the M3 Likelihood

For the quantified records, the derivative of the log-likelihood with respect to \(\theta\) is \(\sum f_i(y_i-\theta f_i)/\sigma^2\). For a BQL record, define \(z_i=(L-\theta f_i)/\sigma\). Since \(dz_i/d\theta=-f_i/\sigma\), differentiating \(\log\Phi(z_i)\) gives \(-f_i\phi(z_i)/[\sigma\Phi(z_i)]\), where \(\phi\) is the standard normal density.

With \(\sigma=1\), setting the total derivative to zero gives the score equation

0 = (4-\theta)+0.5(2-0.5\theta) -0.25\frac{\phi(1-0.25\theta)}{\Phi(1-0.25\theta)} -0.125\frac{\phi(1-0.125\theta)}{\Phi(1-0.125\theta)}.

The quantified part simplifies to \(5-1.25\theta\). Solving the full equation numerically gives \(\hat{\theta}_{M3}\approx3.797\) concentration units. At this value, the two BQL standardized LLOQs are approximately \(z_3=1-0.25(3.797)=0.05075\) and \(z_4=1-0.125(3.797)=0.5254\). Thus their probabilities are approximately \(\Phi(0.05075)=0.5202\) and \(\Phi(0.5254)=0.7003\).

For a numerical check on the second probability, at \(z=0.5268\), \(\Phi(z)\approx0.7008\) and \(\ln\{\Phi(z)\}\approx-0.3555\). These rounded values are consistent with the standard normal CDF and its logarithm. The M3 estimate is a maximum, not merely a score-equation solution: the normal log-density is concave in its mean and \(\log\Phi(z)\) is concave in \(z\); with positive \(f_i\) and quantified records, the total log-likelihood is strictly concave in \(\theta\) here, so the root is the unique maximizer.

The calculation also shows how a BQL record can constrain the fit. At \(\hat{\theta}_{M3}\), the predictions at the last two times are about \(0.9493\) and \(0.4746\) concentration units. Each is below the LLOQ as a predicted mean, but neither is treated as an exact observed concentration.

Comparing the M3 estimate with an M4 fit

In the simple M4 implementation used here, each BQL record is replaced by the LLOQ and then treated as an observed value in the usual residual likelihood. This substitution is not the M3 probability contribution. With normal errors of common known \(\sigma\), maximizing the resulting likelihood for \(\theta\) is weighted least squares. When the weights are equal, the estimate is the least-squares fit to the quantified values and the substituted LLOQ values.

Worked Example: M4 Substitution on the Same Profile

Replace the two BQL records with \(1.0\), the LLOQ. The response vector used by this M4 calculation is therefore \((4,2,1,1)\), with prediction factors \((1,0.5,0.25,0.125)\). Minimizing \(\sum_i(y_i-\theta f_i)^2\) gives

\hat{\theta}_{M4} = \frac{\sum_i f_i y_i}{\sum_i f_i^2} = \frac{1(4)+0.5(2)+0.25(1)+0.125(1)} {1^2+0.5^2+0.25^2+0.125^2} = \frac{5.375}{1.328125} \approx 4.047.

For this particular hypothetical profile, the M4 estimate is about \(0.250\) concentration units higher than the M3 estimate of \(3.797\). The substitution pulls the fit toward the LLOQ as though both BQL records had measured values exactly equal to 1.0. The direction and size of the difference are data- and model-dependent; M4 does not always produce a larger estimate than M3.

This comparison is meaningful only because the two fits use the same structural shape, quantified records, and fixed error scale. In a real analysis, changes in parameter estimates may also reflect jointly estimated residual or random-effects parameters, different records included in the data set, or a different likelihood implementation. The comparison should be made by checking the actual objective function and the predictions at the sampling times, not by assuming the method label alone determines the result.

Checking predicted BQL frequencies

A BQL-frequency diagnostic compares observed BQL indicators with the probabilities implied by the fitted model. For a record with predicted mean \(\hat{\mu}_i\), residual standard deviation \(\hat{\sigma}\), and LLOQ \(L_i\), the predicted probability of a BQL result is \(\hat{p}_i=\Phi((L_i-\hat{\mu}_i)/\hat{\sigma})\). Summing these probabilities gives the expected number of BQL records across a set of observations. This is an expectation, not a requirement that the observed count equal it exactly.

Worked Example: Interpreting BQL-Frequency Diagnostics

Use the M3 estimate \(\hat{\theta}=3.797\) and the four records from the sparse-profile example. The predicted means are approximately \(3.797,\;1.899,\;0.9493,\;0.4746\) concentration units. The corresponding BQL probabilities are approximately \(0.00258,\;0.1845,\;0.5202,\;0.7003\), respectively. Their sum is

E(N_{\mathrm{BQL}}) = 0.00258+0.1845+0.5202+0.7003 \approx1.408.

There are two observed BQL records, compared with an expected count of about 1.408 under the fitted model. The observed count is not by itself evidence of a failure: with only four observations, sampling variation is substantial. The pattern is also important. The model assigns very low BQL probability to the first time and high probability to the last, which is consistent with the profile’s two quantified early records and two BQL later records.

In an actual population analysis, examine observed and predicted BQL frequencies across clinically meaningful strata, such as sampling time or assay, and account for repeated records within subjects. A systematic excess or deficit in a time region may indicate a misspecified concentration-time shape, residual model, LLOQ coding, or random-effects distribution. A single aggregate frequency can hide such patterns.

Assumptions, diagnostics, and common failures

The likelihood calculation is only as reliable as the model and data coding behind it. Before interpreting an M3 fit, verify that the BQL flag and LLOQ correspond to the assay and record, the correct residual distribution is used on the correct scale, and the model prediction and residual standard deviation are evaluated consistently for quantified and BQL records. If different records have different LLOQs, use each record’s own \(L_i\). If the residual model is proportional or defined on a transformed concentration scale, calculate the probability of the BQL event on that model’s measurement scale rather than reusing the additive-normal formula without adjustment.

In population PK, BQL status is informative about the latent concentration but does not reveal it. The likelihood must preserve the shared subject-level random effects when combining records, with the conditional contributions integrated over their distribution. Under heavy censoring, the data may contain little information to distinguish parameters that produce similar probabilities of falling below the LLOQ. A numerically converged fit is not proof that all parameters are identified or that the assumed model is adequate.

A practical diagnostic sequence is to verify the contribution for a few records by hand, check whether fitted predictions are plausible at quantified and BQL times, and compare predicted with observed BQL frequencies across time and assay strata. Then inspect objective-function behavior and parameter stability under defensible changes to residual and structural assumptions. Frequency agreement is necessary evidence to consider, not a stand-alone validation: compensating model errors can yield a plausible total BQL count.

Check record classification.
Confirm which observations are quantified, which are BQL, and the LLOQ attached to each record. Do not silently remove BQL observations or treat a missing record as BQL.

Check likelihood scale and inputs.
For a quantified record, verify the density uses the recorded value. For a BQL record, verify the CDF argument is \((L_i-\mu_i)/\sigma\) under the stated additive-normal model.

Check fit and frequency patterns.
Compare predicted and observed BQL rates across time or assay groups, and investigate structured discrepancies rather than relying only on a pooled count.

A frequent implementation mistake is to add the normal density evaluated at the LLOQ for an M3 BQL record. That is not the probability of the observed event. The M3 record says only \(Y<L\), so its contribution is the integral of the density below \(L\), namely \(\Phi((L-\mu)/\sigma)\) under the normal model. Evaluating a density at \(L\) instead assigns the record an exact value and changes the method toward substitution-based handling.

M3 also has limitations. It relies on a valid model for the distribution of measurements and the censoring threshold; it cannot recover an unobserved concentration without assumptions. Sparse profiles and many BQL records may make structural parameters weakly informed, and misspecification can affect both estimates and predicted BQL rates. M4 can be simpler to implement, but its substitution treats censored values as exact at the LLOQ and can distort the likelihood. The choice should be justified by the measurement process and assessed with diagnostics and sensitivity analysis, not by a presumption that one method is universally adequate.

Key takeaway: Under a normal residual model, an M3 BQL record contributes a CDF probability, while a quantified record contributes a density. The fitted model should reproduce BQL patterns reasonably across the profile, but frequency agreement cannot establish that the structural, residual, and random-effects assumptions are correct.

Check Your Understanding

Answer each question, then compare with the solutions below.

For \(L=2.0\) mg/L, \(\mu=2.6\) mg/L, and \(\sigma=0.3\) mg/L, calculate the M3 BQL probability contribution.
In the four-record profile, why is the BQL contribution at the third time not the normal density evaluated at 1.0?
For the sparse-profile example, calculate the M4 estimate using the substituted responses \((4,2,1,1)\) and factors \((1,0.5,0.25,0.125)\).
The fitted model predicts 1.408 BQL records among four records, while two were observed. Does this difference alone show model failure? Explain.
Why must the conditional contributions for a subject generally be multiplied before integrating over that subject’s random effects?

Solutions. (1) The standardized LLOQ is \(z=(2.0-2.6)/0.3=-2.0\), so the BQL contribution is \(\Phi(-2.0)\approx0.02275\). (2) A BQL record identifies the event \(Y<1.0\), not an exact value \(Y=1.0\). Its likelihood contribution is the probability below the threshold, here \(\Phi((1.0-\mu)/\sigma)\), rather than a density at 1.0. (3) The numerator is \(4+1+0.25+0.125=5.375\), and the denominator is \(1+0.25+0.0625+0.015625=1.328125\), giving \(5.375/1.328125\approx4.047\) concentration units. (4) No. An expected count need not equal the observed count in a small sample; inspect patterns by time or assay and assess the model with other diagnostics. (5) Records from the same subject share the random effect, which induces dependence. The joint conditional likelihood must retain that shared effect before averaging over its distribution; separately averaging each record would remove the dependence.