Skills › Scientific Agent Skills › pkpd-modeling › nca.py
nca.py
609 lines · skills/pkpd-modeling/scripts/nca.py at 5f2bc03
1#!/usr/bin/env python3
2"""Non-compartmental analysis of concentration-time data.
4NCA looks arithmetically trivial and is not. Essentially every disagreement
5between two NCA results traces to one of four choices that are rarely written
6down: how lambda_z was selected, which trapezoidal rule was used, what happened
7to BLQ values, and whether AUCinf was based on observed or predicted Clast.
8This script makes all four explicit, reports the diagnostics that decide
9whether the result is reportable, and refuses to hide them behind a default.
11 python3 nca.py -i profile.csv --dose 100 --route extravascular
12 python3 nca.py -i profile.csv --dose-column dose --route iv-infusion --tinf 1
13 python3 nca.py -i ss.csv --dose 50 --route extravascular --tau 12
14 python3 nca.py -i profile.csv --dose 100 --partial-auc 0-24 --partial-auc 0-72
16Input is a table with subject, time, and concentration columns (default names
17``id``, ``time``, ``conc``). Concentrations below the limit of quantification
18may be given as the text ``BLQ`` / ``BQL`` / ``<LLOQ``, or as blanks, or as a
19separate 0/1 ``blq`` column.
20"""
22from __future__ import annotations
24import argparse
25import math
26import sys
27from dataclasses import dataclass
28from typing import Sequence
30import numpy as np
32from _common import (
33 EXIT_INPUT,
34 InputError,
35 Report,
36 add_format_argument,
37 main_wrapper,
38 parse_float,
39 read_table,
40 require_columns,
41)
43BLQ_TOKENS = {"blq", "bql", "<lloq", "lloq", "bloq", "nq", "<loq"}
45ROUTES = {
46 "iv-bolus": "intravenous bolus",
47 "iv-infusion": "intravenous infusion (zero order)",
48 "extravascular": "extravascular (oral, SC, IM, ...)",
49}
52# ----------------------------------------------------------------- parsing
55@dataclass
56class Profile:
57 subject: str
58 time: np.ndarray
59 conc: np.ndarray
60 is_blq: np.ndarray
61 dose: float
62 tau: float | None = None
63 tinf: float | None = None
66def _parse_conc(raw: str, lloq: float | None, rule: str) -> tuple[float | None, bool]:
67 """Return (value, is_blq). ``None`` means 'drop this record'."""
68 text = (raw or "").strip()
69 lowered = text.lower()
70 blq = lowered in BLQ_TOKENS or lowered.startswith("<") or lowered == ""
71 if blq:
72 if rule == "missing":
73 return None, True
74 if rule == "half-lloq":
75 if lloq is None:
76 raise InputError("--blq-rule half-lloq needs --lloq")
77 return lloq / 2.0, True
78 return 0.0, True
79 value = parse_float(text, "conc")
80 if value is None:
81 return None, True
82 if value < 0:
83 raise InputError("concentrations must be non-negative")
84 if lloq is not None and value < lloq:
85 if rule == "missing":
86 return None, True
87 return (lloq / 2.0 if rule == "half-lloq" else 0.0), True
88 return value, False
91def load_profiles(args: argparse.Namespace) -> list[Profile]:
92 rows = read_table(args.input)
93 require_columns(rows, [args.time_column, args.conc_column], str(args.input))
94 has_subject = args.subject_column in rows[0]
95 if args.dose is None and args.dose_column not in rows[0]:
96 raise InputError(
97 f"give a dose with --dose, or a per-subject '{args.dose_column}' column in the input"
98 )
100 buckets: dict[str, list[dict]] = {}
101 for row in rows:
102 subject = row[args.subject_column].strip() if has_subject else "1"
103 buckets.setdefault(subject or "1", []).append(row)
105 profiles = []
106 for subject, subject_rows in buckets.items():
107 times, concs, blqs = [], [], []
108 for index, row in enumerate(subject_rows):
109 time = parse_float(row.get(args.time_column), f"{args.time_column} (row {index + 1})")
110 flagged = str(row.get("blq", "")).strip() in {"1", "y", "yes", "true"}
111 value, is_blq = _parse_conc(row.get(args.conc_column, ""), args.lloq, args.blq_rule)
112 if flagged and not is_blq:
113 value, is_blq = _parse_conc("BLQ", args.lloq, args.blq_rule)
114 if value is None:
115 continue
116 times.append(time)
117 concs.append(value)
118 blqs.append(is_blq)
119 if not times:
120 raise InputError(f"subject {subject}: no usable concentration records")
122 order = np.argsort(np.asarray(times, dtype=float), kind="stable")
123 time_array = np.asarray(times, dtype=float)[order]
124 conc_array = np.asarray(concs, dtype=float)[order]
125 blq_array = np.asarray(blqs, dtype=bool)[order]
127 if args.dose is not None:
128 dose = args.dose
129 else:
130 dose = parse_float(subject_rows[0].get(args.dose_column), args.dose_column)
131 if dose is None or dose <= 0:
132 raise InputError(f"subject {subject}: dose must be positive")
134 tau = args.tau
135 if args.tau is None and "tau" in subject_rows[0]:
136 tau = parse_float(subject_rows[0].get("tau"), "tau", allow_missing=True)
137 tinf = args.tinf
138 if args.tinf is None and "tinf" in subject_rows[0]:
139 tinf = parse_float(subject_rows[0].get("tinf"), "tinf", allow_missing=True)
141 if np.any(time_array < 0) or np.any(np.diff(time_array) <= 0):
142 raise InputError(f"subject {subject}: times must be non-negative and distinct")
143 if tau is not None and (not math.isfinite(tau) or tau <= 0):
144 raise InputError("tau must be finite and positive")
145 if args.route == "iv-infusion" and (tinf is None or not math.isfinite(tinf) or tinf <= 0):
146 raise InputError("iv-infusion needs a positive --tinf or tinf column")
147 profiles.append(Profile(subject, time_array, conc_array, blq_array, float(dose), tau, tinf))
148 return profiles
151# ------------------------------------------------------------- trapezoidal
154def _segment(t0: float, t1: float, c0: float, c1: float, method: str) -> tuple[float, float]:
155 """AUC and AUMC over one interval. Returns (auc, aumc)."""
156 dt = t1 - t0
157 if dt <= 0:
158 return 0.0, 0.0
159 use_log = (
160 method in {"linup-logdown", "log"}
161 and c0 > 0
162 and c1 > 0
163 and (c1 < c0 or method == "log")
164 and not math.isclose(c0, c1)
165 )
166 if not use_log:
167 return (c0 + c1) / 2.0 * dt, (t0 * c0 + t1 * c1) / 2.0 * dt
168 k = math.log(c0 / c1) / dt
169 auc = (c0 - c1) / k
170 aumc = (t0 * c0 - t1 * c1) / k + (c0 - c1) / (k * k)
171 return auc, aumc
174def cumulative_auc(time: np.ndarray, conc: np.ndarray, method: str) -> tuple[np.ndarray, np.ndarray]:
175 auc = np.zeros_like(time)
176 aumc = np.zeros_like(time)
177 for i in range(1, len(time)):
178 d_auc, d_aumc = _segment(time[i - 1], time[i], conc[i - 1], conc[i], method)
179 auc[i] = auc[i - 1] + d_auc
180 aumc[i] = aumc[i - 1] + d_aumc
181 return auc, aumc
184def interpolate(time: np.ndarray, conc: np.ndarray, target: float, method: str) -> float:
185 """Concentration at an arbitrary time, consistent with the AUC rule."""
186 if target <= time[0]:
187 return float(conc[0])
188 if target >= time[-1]:
189 return float(conc[-1])
190 i = int(np.searchsorted(time, target))
191 t0, t1, c0, c1 = time[i - 1], time[i], conc[i - 1], conc[i]
192 use_log = method in {"linup-logdown", "log"} and c0 > 0 and c1 > 0 and (c1 < c0 or method == "log")
193 if not use_log:
194 return float(c0 + (c1 - c0) * (target - t0) / (t1 - t0))
195 k = math.log(c0 / c1) / (t1 - t0)
196 return float(c0 * math.exp(-k * (target - t0)))
199def partial_auc(time: np.ndarray, conc: np.ndarray, start: float, end: float, method: str) -> float:
200 if start < time[0] or end > time[-1] or end <= start:
201 raise InputError("partial AUC needs both boundaries within sampled times and end > start")
202 inner = [t for t in time if start < t < end]
203 knots = [start, *inner, end]
204 values = [interpolate(time, conc, t, method) for t in knots]
205 total = 0.0
206 for i in range(1, len(knots)):
207 total += _segment(knots[i - 1], knots[i], values[i - 1], values[i], method)[0]
208 return total
211# --------------------------------------------------------------- lambda_z
214@dataclass
215class LambdaZ:
216 lam: float | None = None
217 intercept: float | None = None
218 r2: float | None = None
219 r2_adj: float | None = None
220 n_points: int = 0
221 t_first: float | None = None
222 t_last: float | None = None
223 clast_pred: float | None = None
224 reason: str = ""
225 excluded_cmax: bool = True
227 @property
228 def half_life(self) -> float | None:
229 return math.log(2.0) / self.lam if self.lam and self.lam > 0 else None
232def _loglinear_fit(t: np.ndarray, c: np.ndarray) -> tuple[float, float, float]:
233 """Slope, intercept, r2 of ln(C) on time."""
234 if len(t) < 2 or np.ptp(t) <= 0:
235 raise InputError("lambda_z needs distinct sample times")
236 y = np.log(c)
237 tbar, ybar = t.mean(), y.mean()
238 sxx = float(np.sum((t - tbar) ** 2))
239 slope = float(np.sum((t - tbar) * (y - ybar)) / sxx)
240 intercept = float(ybar - slope * tbar)
241 fitted = intercept + slope * t
242 ss_res = float(np.sum((y - fitted) ** 2))
243 ss_tot = float(np.sum((y - ybar) ** 2))
244 r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
245 return slope, intercept, r2
248def estimate_lambda_z(
249 time: np.ndarray,
250 conc: np.ndarray,
251 is_blq: np.ndarray,
252 tmax: float,
253 min_points: int = 3,
254 manual: tuple[float, float] | None = None,
255) -> LambdaZ:
256 """Select the terminal window by best adjusted r-squared.
258 This helper starts from the last three quantifiable points, extends
259 backwards, and retains a longer window only when adjusted r-squared
260 improves by more than 0.0001. Unlike Phoenix Best Fit, it does not favor
261 the longer window when scores are within that tolerance. Neither plain
262 nor adjusted r-squared is monotonic as observations are added.
264 Points at or before Tmax are never eligible in this helper. Distribution
265 and absorption can still contaminate the tail; bias direction is not fixed.
266 """
267 if min_points < 3:
268 raise InputError("lambda_z requires at least 3 points")
269 eligible = (~is_blq) & (conc > 0) & (time > tmax)
270 idx = np.flatnonzero(eligible)
271 if manual is not None:
272 lo, hi = manual
273 idx = np.flatnonzero(eligible & (time >= lo) & (time <= hi))
274 if len(idx) < 3:
275 return LambdaZ(reason=f"manual window {lo}-{hi} contains {len(idx)} quantifiable points")
276 slope, intercept, r2 = _loglinear_fit(time[idx], conc[idx])
277 if slope >= 0:
278 return LambdaZ(reason="manual terminal slope is not negative", n_points=len(idx))
279 n = len(idx)
280 r2_adj = 1.0 - (1.0 - r2) * (n - 1) / (n - 2) if n > 2 else float("nan")
281 return LambdaZ(
282 lam=-slope,
283 intercept=intercept,
284 r2=r2,
285 r2_adj=r2_adj,
286 n_points=n,
287 t_first=float(time[idx][0]),
288 t_last=float(time[idx][-1]),
289 clast_pred=float(math.exp(intercept + slope * time[np.flatnonzero((~is_blq) & (conc > 0))[-1]])),
290 reason="manual window",
291 excluded_cmax=bool(np.all(time[idx] > tmax)),
292 )
294 if len(idx) < min_points:
295 return LambdaZ(reason=f"only {len(idx)} quantifiable points after Tmax; need {min_points}")
297 best: LambdaZ | None = None
298 for start in range(len(idx) - min_points, -1, -1):
299 window = idx[start:]
300 n = len(window)
301 slope, intercept, r2 = _loglinear_fit(time[window], conc[window])
302 if n <= 2 or not math.isfinite(r2) or slope >= 0:
303 continue
304 r2_adj = 1.0 - (1.0 - r2) * (n - 1) / (n - 2)
305 candidate = LambdaZ(
306 lam=-slope,
307 intercept=intercept,
308 r2=r2,
309 r2_adj=r2_adj,
310 n_points=n,
311 t_first=float(time[window][0]),
312 t_last=float(time[window][-1]),
313 clast_pred=float(math.exp(intercept + slope * time[window][-1])),
314 reason="best adjusted r2",
315 )
316 if best is None or (candidate.r2_adj or -math.inf) > (best.r2_adj or -math.inf) + 1e-4:
317 best = candidate
318 if best is None:
319 return LambdaZ(reason="no window of at least 3 points could be fitted")
320 if best.lam is not None and best.lam <= 0:
321 return LambdaZ(reason=f"terminal slope is not negative (lambda_z = {best.lam:.4g})", n_points=best.n_points)
322 return best
325# ------------------------------------------------------------------- NCA
328def analyse(profile: Profile, args: argparse.Namespace) -> tuple[dict, LambdaZ, list[str]]:
329 findings: list[str] = []
330 time, conc, is_blq = profile.time, profile.conc, profile.is_blq
331 quantifiable = (~is_blq) & (conc > 0)
333 if not np.any(quantifiable):
334 findings.append(f"subject {profile.subject}: every sample is BLQ; no parameters computed")
335 return {"id": profile.subject}, LambdaZ(reason="all BLQ"), findings
337 peak = int(np.argmax(np.where(quantifiable, conc, -np.inf)))
338 cmax, tmax = float(conc[peak]), float(time[peak])
339 last = int(np.flatnonzero(quantifiable)[-1])
340 clast, tlast = float(conc[last]), float(time[last])
342 if time[0] > 0:
343 findings.append(f"subject {profile.subject}: first sample is at {time[0]:g}; AUC starts there, not at dose time 0; CL and volume estimates omit early exposure")
344 auc_cum, aumc_cum = cumulative_auc(time, conc, args.auc_method)
345 auc_last, aumc_last = float(auc_cum[last]), float(aumc_cum[last])
347 lz = estimate_lambda_z(time, conc, is_blq, tmax, args.lambda_z_points, args.lambda_z_window)
349 row: dict[str, object] = {
350 "id": profile.subject,
351 "dose": profile.dose,
352 "n_obs": int(len(time)),
353 "n_blq": int(np.sum(is_blq)),
354 "cmax": cmax,
355 "tmax": tmax,
356 "clast": clast,
357 "tlast": tlast,
358 "auc_last": auc_last,
359 }
360 if profile.tau:
361 row["auc_tau"] = partial_auc(time, conc, 0.0, profile.tau, args.auc_method)
363 if lz.lam and profile.tau is None:
364 lam = lz.lam
365 row["lambda_z"] = lam
366 row["t_half"] = lz.half_life
367 row["r2_adj"] = lz.r2_adj
368 row["lambda_z_n"] = lz.n_points
369 auc_inf_obs = auc_last + clast / lam
370 auc_inf_pred = auc_last + (lz.clast_pred or clast) / lam
371 aumc_inf = aumc_last + tlast * clast / lam + clast / lam**2
372 row["auc_inf_obs"] = auc_inf_obs
373 row["auc_inf_pred"] = auc_inf_pred
374 row["pct_auc_extrap"] = 100.0 * (auc_inf_obs - auc_last) / auc_inf_obs
375 row["pct_aumc_extrap"] = 100.0 * (aumc_inf - aumc_last) / aumc_inf if aumc_inf > 0 else float("nan")
377 mrt = aumc_inf / auc_inf_obs
378 if args.route == "iv-infusion" and profile.tinf:
379 mrt -= profile.tinf / 2.0
380 row["mrt"] = mrt
382 clearance = profile.dose / auc_inf_obs
383 row["cl_f" if args.route == "extravascular" else "cl"] = clearance
384 row["vz_f" if args.route == "extravascular" else "vz"] = clearance / lam
385 if args.route != "extravascular":
386 row["vss"] = clearance * mrt
388 span = (lz.t_last - lz.t_first) / lz.half_life
389 row["span_ratio"] = span
391 prefix = f"subject {profile.subject}"
392 if row["pct_auc_extrap"] > args.max_extrap:
393 findings.append(
394 f"{prefix}: {row['pct_auc_extrap']:.1f}% of AUCinf is extrapolated "
395 f"(above {args.max_extrap:.0f}%); AUCinf is driven by the lambda_z fit, not by data"
396 )
397 if lz.r2_adj is not None and lz.r2_adj < args.min_r2_adj:
398 findings.append(
399 f"{prefix}: terminal-phase adjusted r2 is {lz.r2_adj:.4f} (below {args.min_r2_adj}); "
400 "half-life, Vz and AUCinf inherit that uncertainty"
401 )
402 if span < args.min_span:
403 findings.append(
404 f"{prefix}: lambda_z window spans {span:.2f} half-lives (below {args.min_span}); "
405 "the terminal phase may not have been reached"
406 )
407 if lz.n_points < 3:
408 findings.append(f"{prefix}: lambda_z estimated from {lz.n_points} points")
409 elif not lz.lam:
410 findings.append(f"subject {profile.subject}: lambda_z not estimable - {lz.reason}")
412 if profile.tau:
413 if np.any(time > profile.tau):
414 raise InputError("--tau input must contain only one interval [0, tau]")
415 findings.append(f"subject {profile.subject}: --tau assumes demonstrated steady state; AUCinf and single-dose CL/Vz/Vss are withheld")
416 tau = profile.tau
417 auc_tau = float(row["auc_tau"])
418 cmin = float(np.min(conc[(time >= 0) & (time <= tau) & quantifiable])) if np.any(quantifiable) else float("nan")
419 cavg = auc_tau / tau
420 row["cavg_ss"] = cavg
421 row["cmin_ss"] = cmin
422 row["ptf_pct"] = 100.0 * (cmax - cmin) / cavg if cavg else float("nan")
423 row["swing"] = (cmax - cmin) / cmin if cmin else float("nan")
424 row["cl_ss_f" if args.route == "extravascular" else "cl_ss"] = profile.dose / auc_tau
425 if lz.lam:
426 row["terminal_slope_accumulation_index"] = 1.0 / (-math.expm1(-lz.lam * tau))
427 findings.append("terminal-slope accumulation index is a monoexponential heuristic, not the observed AUC/peak accumulation ratio")
429 for start, end in args.partial_auc:
430 if end > tlast:
431 findings.append(
432 f"subject {profile.subject}: partial AUC {start:g}-{end:g} extends past the last "
433 f"quantifiable sample at {tlast:g}; integration uses the declared BLQ substitutions beyond Clast"
434 )
435 row[f"auc_{start:g}_{end:g}"] = partial_auc(time, conc, start, end, args.auc_method)
437 prefix = f"subject {profile.subject}"
438 if peak == len(time) - 1:
439 findings.append(f"{prefix}: Cmax is the last sample; the peak and the terminal phase are not characterised")
440 if peak == 0 and args.route == "extravascular":
441 findings.append(f"{prefix}: Cmax is the first sample; the true peak may precede the first draw")
442 if is_blq[:peak].any() and args.blq_rule == "zero":
443 row["leading_blq_set_to_zero"] = int(np.sum(is_blq[:peak]))
445 return row, lz, findings
448# --------------------------------------------------------------- summaries
451def summarise(rows: Sequence[dict], keys: Sequence[str]) -> list[dict]:
452 out = []
453 for key in keys:
454 values = [float(r[key]) for r in rows if isinstance(r.get(key), (int, float)) and math.isfinite(float(r[key]))]
455 if len(values) < 1:
456 continue
457 arr = np.asarray(values)
458 entry: dict[str, object] = {
459 "parameter": key,
460 "n": len(arr),
461 "mean": float(arr.mean()),
462 "sd": float(arr.std(ddof=1)) if len(arr) > 1 else float("nan"),
463 "cv_pct": float(100.0 * arr.std(ddof=1) / arr.mean()) if len(arr) > 1 and arr.mean() else float("nan"),
464 "median": float(np.median(arr)),
465 "min": float(arr.min()),
466 "max": float(arr.max()),
467 }
468 if np.all(arr > 0):
469 logs = np.log(arr)
470 entry["geo_mean"] = float(np.exp(logs.mean()))
471 entry["geo_cv_pct"] = (
472 float(100.0 * math.sqrt(math.exp(float(logs.var(ddof=1))) - 1.0)) if len(arr) > 1 else float("nan")
473 )
474 out.append(entry)
475 return out
478# --------------------------------------------------------------------- CLI
481def _partial_spec(text: str) -> tuple[float, float]:
482 try:
483 start, end = text.split("-", 1)
484 lo, hi = float(start), float(end)
485 except ValueError as exc:
486 raise argparse.ArgumentTypeError(f"partial AUC must look like 0-24, got {text!r}") from exc
487 if hi <= lo:
488 raise argparse.ArgumentTypeError(f"partial AUC end must exceed start, got {text!r}")
489 return lo, hi
492def _window_spec(text: str) -> tuple[float, float]:
493 return _partial_spec(text)
496def build_parser() -> argparse.ArgumentParser:
497 parser = argparse.ArgumentParser(
498 description="Non-compartmental analysis with explicit lambda_z, BLQ, and trapezoidal choices.",
499 formatter_class=argparse.RawDescriptionHelpFormatter,
500 epilog=__doc__,
501 )
502 parser.add_argument("-i", "--input", required=True, help="concentration-time table (CSV/TSV, or - for stdin)")
503 parser.add_argument("--dose", type=float, help="dose given to every subject")
504 parser.add_argument("--dose-column", default="dose", help="per-subject dose column (default: dose)")
505 parser.add_argument("--subject-column", default="id")
506 parser.add_argument("--time-column", default="time")
507 parser.add_argument("--conc-column", default="conc")
508 parser.add_argument(
509 "--route",
510 choices=sorted(ROUTES),
511 default="extravascular",
512 help="route; decides whether CL/Vz are apparent (/F) and whether Vss is reported",
513 )
514 parser.add_argument("--tinf", type=float, help="infusion duration, for --route iv-infusion")
515 parser.add_argument("--tau", type=float, help="dosing interval; requests steady-state parameters")
516 parser.add_argument(
517 "--auc-method",
518 choices=("linup-logdown", "linear", "log"),
519 default="linup-logdown",
520 help="trapezoidal rule (default: linup-logdown, the usual regulatory choice)",
521 )
522 parser.add_argument("--lloq", type=float, help="lower limit of quantification")
523 parser.add_argument(
524 "--blq-rule",
525 choices=("zero", "half-lloq", "missing"),
526 default="zero",
527 help="how BLQ samples enter the AUC (default: zero)",
528 )
529 parser.add_argument("--lambda-z-points", type=int, default=3, help="minimum points in the terminal fit (default: 3)")
530 parser.add_argument("--lambda-z-window", type=_window_spec, help="force the terminal window, e.g. 8-48")
531 parser.add_argument("--partial-auc", type=_partial_spec, action="append", default=[], help="e.g. --partial-auc 0-24")
532 parser.add_argument("--min-r2-adj", type=float, default=0.80, help="flag terminal fits below this (default: 0.80)")
533 parser.add_argument("--max-extrap", type=float, default=20.0, help="flag %%AUC extrapolated above this (default: 20)")
534 parser.add_argument("--min-span", type=float, default=2.0, help="flag lambda_z spans below this many half-lives (default: 2)")
535 parser.add_argument("--no-summary", action="store_true", help="per-subject table only")
536 add_format_argument(parser)
537 return parser
540def run(argv: Sequence[str] | None = None) -> int:
541 args = build_parser().parse_args(argv)
542 from _common import validate_numeric_args
543 validate_numeric_args(args)
544 if args.lloq is not None and args.lloq <= 0:
545 raise InputError("--lloq must be positive")
547 profiles = load_profiles(args)
548 report = Report()
549 report.note(f"route: {ROUTES[args.route]}")
550 report.note(f"trapezoidal rule: {args.auc_method}; BLQ rule: {args.blq_rule}")
551 report.note(
552 "AUCinf_obs uses the observed Clast; AUCinf_pred uses the value predicted by the "
553 "lambda_z fit. Report which one you used - they are not interchangeable."
554 )
556 rows, diagnostics, all_findings = [], [], []
557 for profile in profiles:
558 row, lz, findings = analyse(profile, args)
559 rows.append(row)
560 diagnostics.append(
561 {
562 "id": profile.subject,
563 "lambda_z": lz.lam,
564 "t_half": lz.half_life,
565 "n_points": lz.n_points,
566 "window_start": lz.t_first,
567 "window_end": lz.t_last,
568 "r2": lz.r2,
569 "r2_adj": lz.r2_adj,
570 "clast_pred": lz.clast_pred,
571 "basis": lz.reason,
572 }
573 )
574 all_findings.extend(findings)
576 report.table("per-subject parameters", rows)
577 report.table("lambda_z diagnostics", diagnostics)
579 if not args.no_summary and len(rows) > 1:
580 candidates = [
581 "cmax",
582 "tmax",
583 "auc_last",
584 "auc_inf_obs",
585 "t_half",
586 "cl_f",
587 "cl",
588 "vz_f",
589 "vz",
590 "vss",
591 "mrt",
592 "auc_tau",
593 "cavg_ss",
594 ]
595 keys = [k for k in candidates if any(k in r for r in rows)]
596 keys += [k for k in rows[0] if k.startswith("auc_") and k not in keys and k not in {"auc_last", "auc_tau"}]
597 report.table("summary statistics", summarise(rows, keys))
598 report.note(
599 "Exposure metrics (AUC, Cmax) are conventionally summarised as geometric mean with "
600 "geometric CV%; Tmax as median and range. Arithmetic statistics are shown alongside."
601 )
603 for finding in all_findings:
604 report.finding(finding)
605 return report.emit(args.format)
608if __name__ == "__main__":
609 raise SystemExit(main_wrapper(run))