BLOQ (M3 method)

This example shows how to fit a model to data that contains observations below the assay lower limit of quantification (BLOQ) using Beal’s M3 method. Instead of dropping BLOQ rows — which biases terminal-phase parameter estimates — each below-LOQ observation contributes P(y < LLOQ | θ, η) = Φ((LLOQ − f)/√V) to the likelihood. Above-LOQ rows use the mirrored upper-tail contribution, P(y > ULOQ | θ, η) = Φ((f − ULOQ)/√V).

Dataset (data/warfarin_bloq.csv)

This is the warfarin dataset with an added CENS column. Ten late-time observations that originally fell below an assay LLOQ of 2.0 µg/mL have been marked CENS=1 and their DV cell set to the LLOQ value:

ID,TIME,DV,EVID,AMT,CMT,RATE,MDV,CENS
1,0,.,1,100,1,0,1,0
1,0.5,5.3653,0,.,1,0,0,0
...
1,96,2.5019,0,.,1,0,0,0
1,120,2,0,.,1,0,0,1

CENS convention: when CENS=1, the row is left-censored and DV carries the LLOQ value — not the true (unobserved) concentration. When CENS=-1, the row is right-censored and DV carries the ULOQ value. This matches the Monolix/nlmixr2 polarity. Rows with no CENS column, or CENS=0, are treated as ordinary quantified observations.

Model File (examples/warfarin_bloq.ferx)

# One-compartment oral PK model (warfarin) with M3 BLOQ likelihood.

[parameters]
  theta TVCL(0.2, 0.001, 10.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.5, 0.01, 50.0)

  omega ETA_CL ~ 0.09
  omega ETA_V  ~ 0.04
  omega ETA_KA ~ 0.30

  sigma PROP_ERR ~ 0.02

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  DV ~ proportional(PROP_ERR)

[fit_options]
  method       = focei
  maxiter      = 300
  covariance   = true
  bloq_method  = m3

The only two lines that differ from a plain fit are bloq_method = m3 in [fit_options] and the added CENS column in the CSV. Nothing else in the model changes.

Running

ferx examples/warfarin_bloq.ferx --data data/warfarin_bloq.csv

Expected Results

--- Objective Function ---
OFV:  -217.1841

--- THETA Estimates ---
TVCL                 0.132811
TVV                  7.732354
TVKA                 0.811661

--- OMEGA Estimates ---
  OMEGA(1,1) = 0.029021  (CV% = 17.0)
  OMEGA(2,2) = 0.009557  (CV% = 9.8)
  OMEGA(3,3) = 0.338106  (CV% = 58.1)

--- SIGMA Estimates ---
  SIGMA(1) = 0.010761

These match a NONMEM 7.5.1 LAPLACE M3 reference (tests/nonmem/warfarin_bloq.{ctl,lst}) to ~4 significant figures on the structural and variance parameters (TVCL 0.132801, TVV 7.73139, TVKA 0.809824, PROP 0.010760, ω 0.028849 / 0.009544 / 0.335772). The NONMEM objective uses the F_FLAG likelihood convention for censored rows, which carries a different additive constant than ferx’s M3 term, so the OFV is not directly comparable (the MLE is).

The diagnostic table (warfarin_bloq-sdtab.csv) gains a CENS column, and the IWRES / CWRES cells for censored rows are written as empty (a weighted Gaussian residual is undefined when the observed value is censored). This blanking follows bloq_method, not the column: under the default drop the same row is fitted as an ordinary observation at its DV and keeps its residuals, so the diagnostics always describe the likelihood that was actually fitted (#1499).

Method Notes

  • Activation requires both pieces: the CENS column in the data file and bloq_method = m3 in [fit_options]. Without the option, CENS=1 (below LLOQ) and CENS=-1 (above ULOQ) rows are treated as ordinary observations at the DV limit value, which biases the fit.
  • FOCE and FOCEI give different M3 optima. method = foce keeps a consistent Sheiner–Beal objective: censored rows leave the linearized marginal and re-enter as −log Φ((LLOQ − f̂)/√R⁰) (and the ULOQ mirror for CENS=-1) with the population (η=0) variance, matching NONMEM METHOD=1 LAPLACE without INTER. method = focei evaluates the censored term at the conditional variance (NONMEM with INTER). Both have exact analytic gradients. On warfarin BLOQ the two land at meaningfully different TVKA (FOCE ≈ 0.71, FOCEI ≈ 0.81), exactly as the corresponding NONMEM runs do — pick the method to match your reference fit. (Earlier versions silently promoted censored subjects to FOCEI under method = foce.)
  • Gauss-Newton caveat. With method = gn or gn_hybrid, the BHHH information-matrix approximation degrades as the BLOQ fraction grows (each censored row carries less Fisher information than its Gaussian counterpart). A warning is emitted; for >20% censoring prefer focei.
  • SAEM optimizes θ/σ with M3 in the M-step directly, so no special handling is required beyond setting bloq_method = m3.