Skills › Scientific Agent Skills › pkpd-modeling › fit_compartmental.py
fit_compartmental.py
567 lines · skills/pkpd-modeling/scripts/fit_compartmental.py at 5f2bc03
1#!/usr/bin/env python3
2"""Fit compartmental PK models to concentration-time data, with identifiability diagnostics.
4The fit itself is the easy part. What decides whether the result means anything
5is: was the residual error model right, is the extra compartment actually
6supported, and are the parameters identifiable from these data at all. This
7script reports all three, and flags the cases where a converged fit is still
8uninterpretable.
10 python3 fit_compartmental.py -i profile.csv --dose 100 --route iv-bolus --model 2cmt
11 python3 fit_compartmental.py -i profile.csv --dose 100 --route oral --compare 1cmt,2cmt
12 python3 fit_compartmental.py -i profile.csv --dose 500 --route iv-infusion --tinf 1 \
13 --model 2cmt --weight 1/y2
15Input columns default to ``time`` and ``conc``; an ``id`` column fits each
16subject separately. Parameters are estimated on the log scale, so they cannot
17go negative and their confidence intervals come out asymmetric, which is the
18honest shape for a clearance or a volume.
19"""
21from __future__ import annotations
23import argparse
24import math
25import sys
26from dataclasses import dataclass
27from typing import Callable, Sequence
29import numpy as np
31from _common import (
32 EXIT_INPUT,
33 InputError,
34 Report,
35 add_format_argument,
36 group_by,
37 main_wrapper,
38 parse_float,
39 read_table,
40 require_columns,
41)
42from _models import Dose, build_regimen, disposition, simulate_linear
44try:
45 from scipy.optimize import least_squares
46 from scipy.stats import f as f_dist
47 from scipy.stats import t as t_dist
48except ImportError as exc: # pragma: no cover
49 raise SystemExit("fit_compartmental.py needs scipy: uv pip install scipy") from exc
52MODELS = {
53 "1cmt": ["CL", "V1"],
54 "2cmt": ["CL", "V1", "Q2", "V2"],
55 "3cmt": ["CL", "V1", "Q2", "V2", "Q3", "V3"],
56}
58WEIGHTS = {
59 "uniform": "constant variance (additive error)",
60 "1/y": "variance proportional to the observation",
61 "1/y2": "constant CV (proportional error) - the usual PK default",
62 "1/yhat": "variance proportional to the prediction",
63 "1/yhat2": "constant CV against the prediction",
64}
67# ------------------------------------------------------------------- model
70@dataclass
71class FitSpec:
72 model: str
73 route: str
74 dose: float
75 tinf: float | None
76 interval: float | None
77 n_doses: int
78 absorption: str
79 estimate_lag: bool
81 @property
82 def parameter_names(self) -> list[str]:
83 names = list(MODELS[self.model])
84 if self.route == "oral":
85 names = [f"{n}/F" for n in names]
86 names.append("MTT" if self.absorption == "transit" else "ka")
87 if self.absorption == "transit":
88 names.append("Ntr")
89 if self.estimate_lag:
90 names.append("tlag")
91 return names
94def _regimen(spec: FitSpec) -> list[Dose]:
95 return build_regimen(
96 spec.dose,
97 interval=spec.interval,
98 n_doses=spec.n_doses,
99 duration=spec.tinf or 0.0,
100 route="oral" if spec.route == "oral" else "iv",
101 )
104def predictor(spec: FitSpec) -> Callable[[np.ndarray, np.ndarray], np.ndarray]:
105 """Return f(times, theta) -> predicted concentration."""
106 regimen = _regimen(spec)
107 n_disp = len(MODELS[spec.model])
109 def predict(times: np.ndarray, theta: np.ndarray) -> np.ndarray:
110 cl, v1 = theta[0], theta[1]
111 q = tuple(theta[2:n_disp:2])
112 vp = tuple(theta[3:n_disp:2])
113 disp = disposition(cl, v1, q, vp)
114 if spec.route != "oral":
115 return simulate_linear(times, regimen, disp)
116 extra = theta[n_disp:]
117 if spec.absorption == "transit":
118 mtt, ntr = extra[0], extra[1]
119 tlag = extra[2] if spec.estimate_lag else 0.0
120 from _models import conc_transit # local: only this branch needs it
122 total = np.zeros_like(times, dtype=float)
123 for dose in regimen:
124 total += conc_transit(np.maximum(times - dose.time - tlag, 0.0), dose.amount, mtt, ntr, disp)
125 return total
126 ka = extra[0]
127 tlag = extra[1] if spec.estimate_lag else 0.0
128 return simulate_linear(times, regimen, disp, ka=ka, f=1.0, tlag=tlag)
130 return predict
133# ------------------------------------------------------------------ fitting
136@dataclass
137class FitResult:
138 theta: np.ndarray
139 names: list[str]
140 se_log: np.ndarray
141 wssr: float
142 n: int
143 p: int
144 correlation: np.ndarray
145 condition: float
146 predicted: np.ndarray
147 residuals: np.ndarray
148 weighted_residuals: np.ndarray
149 sigma: float
150 success: bool
151 message: str
153 @property
154 def aic(self) -> float:
155 return self.n * math.log(max(self.wssr / self.n, np.finfo(float).tiny)) + 2 * self.p
157 @property
158 def bic(self) -> float:
159 return self.n * math.log(max(self.wssr / self.n, np.finfo(float).tiny)) + self.p * math.log(self.n)
162def _weights(y: np.ndarray, yhat: np.ndarray, scheme: str) -> np.ndarray:
163 floor = 1e-12
164 if scheme == "uniform":
165 return np.ones_like(y)
166 if scheme == "1/y":
167 return 1.0 / np.maximum(np.abs(y), floor)
168 if scheme == "1/y2":
169 return 1.0 / np.maximum(y**2, floor)
170 if scheme == "1/yhat":
171 return 1.0 / np.maximum(np.abs(yhat), floor)
172 if scheme == "1/yhat2":
173 return 1.0 / np.maximum(yhat**2, floor)
174 raise InputError(f"unknown weighting scheme {scheme!r}")
177def initial_estimates(spec: FitSpec, time: np.ndarray, conc: np.ndarray) -> np.ndarray:
178 """Heuristic starting values from the shape of the data.
180 Bad initial estimates are the usual cause of a 'model that will not fit'.
181 These come from the same quantities NCA would give: the terminal slope sets
182 the slow disposition, the AUC sets clearance, and the peak sets the central
183 volume.
184 """
185 positive = conc > 0
186 t, c = time[positive], conc[positive]
187 if len(t) < 3:
188 raise InputError("need at least 3 positive concentrations to fit")
189 tail = slice(max(len(t) - 3, 0), len(t))
190 slope = np.polyfit(t[tail], np.log(c[tail]), 1)[0]
191 lam = max(-slope, 1e-4)
192 auc = float(np.trapezoid(c, t)) + c[-1] / lam
193 cl = spec.dose / max(auc, 1e-9)
194 cmax = float(c.max())
195 v1 = spec.dose / cmax if spec.route == "iv-bolus" else cl / lam
196 v1 = max(v1, 1e-6)
198 theta = [cl, v1]
199 if spec.model in {"2cmt", "3cmt"}:
200 theta += [cl, v1 * 2.0]
201 if spec.model == "3cmt":
202 theta += [cl / 3.0, v1 * 5.0]
203 if spec.route == "oral":
204 tmax = float(t[int(np.argmax(c))])
205 if spec.absorption == "transit":
206 theta += [max(tmax * 0.7, 1e-3), 3.0]
207 else:
208 theta += [max(2.0 / max(tmax, 1e-3), 1e-3)]
209 if spec.estimate_lag:
210 theta += [max(float(t.min()) * 0.5, 1e-3)]
211 return np.asarray(theta, dtype=float)
214def fit_one(
215 spec: FitSpec,
216 time: np.ndarray,
217 conc: np.ndarray,
218 scheme: str,
219 starts: int = 5,
220 seed: int = 20260727,
221) -> FitResult:
222 predict = predictor(spec)
223 names = spec.parameter_names
224 p = len(names)
225 n = len(time)
226 if n <= p:
227 raise InputError(f"{n} observations cannot support {p} parameters")
229 base = initial_estimates(spec, time, conc)
230 rng = np.random.default_rng(seed)
232 def residual(log_theta: np.ndarray) -> np.ndarray:
233 theta = np.exp(log_theta)
234 try:
235 yhat = predict(time, theta)
236 except (ValueError, np.linalg.LinAlgError):
237 return np.full(n, 1e6)
238 if not np.all(np.isfinite(yhat)):
239 return np.full(n, 1e6)
240 w = _weights(conc, yhat, scheme)
241 return np.sqrt(w) * (conc - yhat)
243 best = None
244 for attempt in range(max(1, starts)):
245 start = base if attempt == 0 else base * np.exp(rng.normal(0.0, 0.5, size=p))
246 try:
247 candidate = least_squares(
248 residual,
249 np.log(np.maximum(start, 1e-12)),
250 method="trf",
251 x_scale="jac",
252 max_nfev=5000 * p,
253 )
254 except (ValueError, np.linalg.LinAlgError):
255 continue
256 if best is None or candidate.cost < best.cost:
257 best = candidate
258 if best is None: # pragma: no cover - every start failed
259 raise InputError("optimiser failed from every starting point; check dose, units and route")
261 theta = np.exp(best.x)
262 yhat = predict(time, theta)
263 weights = _weights(conc, yhat, scheme)
264 residuals = conc - yhat
265 weighted = np.sqrt(weights) * residuals
266 wssr = float(np.sum(weighted**2))
267 dof = n - p
268 sigma2 = wssr / dof
270 jac = best.jac
271 try:
272 # Covariance on the log scale from the Gauss-Newton approximation.
273 _, s, vt = np.linalg.svd(jac, full_matrices=False)
274 threshold = np.finfo(float).eps * max(jac.shape) * (s[0] if s.size else 0.0)
275 if np.count_nonzero(s > threshold) < p:
276 raise np.linalg.LinAlgError("rank-deficient sensitivity matrix; covariance is not identifiable")
277 s_inv = 1.0 / s
278 cov = (vt.T * s_inv**2) @ vt * sigma2
279 se_log = np.sqrt(np.maximum(np.diag(cov), 0.0))
280 outer = np.outer(se_log, se_log)
281 with np.errstate(divide="ignore", invalid="ignore"):
282 corr = np.where(outer > 0, cov / np.where(outer > 0, outer, 1.0), np.nan)
283 # NONMEM's condition number: largest over smallest eigenvalue of the
284 # *correlation* matrix of the estimates. Reported this way so the
285 # familiar ">1000 is ill-conditioned" rule of thumb actually applies;
286 # the condition number of the Jacobian is a different quantity on a
287 # different scale, and the two get confused constantly.
288 if np.all(np.isfinite(corr)):
289 eigenvalues = np.linalg.eigvalsh(corr)
290 smallest = float(eigenvalues.min())
291 condition = float(eigenvalues.max() / smallest) if smallest > 0 else float("inf")
292 else:
293 condition = float("inf")
294 except np.linalg.LinAlgError: # pragma: no cover
295 se_log = np.full(p, np.nan)
296 corr = np.full((p, p), np.nan)
297 condition = float("inf")
299 return FitResult(
300 theta=theta,
301 names=names,
302 se_log=se_log,
303 wssr=wssr,
304 n=n,
305 p=p,
306 correlation=corr,
307 condition=condition,
308 predicted=yhat,
309 residuals=residuals,
310 weighted_residuals=weighted,
311 sigma=math.sqrt(sigma2),
312 success=bool(best.success),
313 message=str(best.message),
314 )
317# -------------------------------------------------------------- derivations
320def secondary_parameters(spec: FitSpec, theta: np.ndarray) -> dict[str, float]:
321 n_disp = len(MODELS[spec.model])
322 cl, v1 = theta[0], theta[1]
323 q = tuple(theta[2:n_disp:2])
324 vp = tuple(theta[3:n_disp:2])
325 disp = disposition(cl, v1, q, vp)
326 suffix = "/F" if spec.route == "oral" else ""
327 out = {
328 f"Vss{suffix}": disp.vss,
329 "MRT_iv": disp.mrt_iv,
330 "terminal_t_half": disp.terminal_half_life,
331 }
332 for i, half in enumerate(disp.half_lives, start=1):
333 out[f"t_half_phase{i}"] = float(half)
334 if spec.route == "oral" and spec.absorption != "transit":
335 ka = theta[n_disp]
336 out["absorption_t_half"] = math.log(2.0) / ka
337 if ka < 1.0 / disp.terminal_half_life * math.log(2.0):
338 out["flip_flop_suspected"] = 1.0
339 if spec.interval:
340 from _models import steady_state_metrics
342 out.update({k: v for k, v in steady_state_metrics(disp, spec.dose, spec.interval).items()})
343 return out
346def runs_test_p(residuals: np.ndarray) -> float:
347 """Two-sided Wald-Wolfowitz runs test on residual signs."""
348 signs = np.sign(residuals)
349 signs = signs[signs != 0]
350 n = len(signs)
351 if n < 8:
352 return float("nan")
353 n_pos = int(np.sum(signs > 0))
354 n_neg = n - n_pos
355 if n_pos == 0 or n_neg == 0:
356 return 0.0
357 runs = 1 + int(np.sum(signs[1:] != signs[:-1]))
358 mean = 2.0 * n_pos * n_neg / n + 1.0
359 var = (mean - 1.0) * (mean - 2.0) / (n - 1.0)
360 if var <= 0:
361 return float("nan")
362 z = (runs - mean) / math.sqrt(var)
363 return float(2.0 * (1.0 - 0.5 * (1.0 + math.erf(abs(z) / math.sqrt(2.0)))))
366# --------------------------------------------------------------------- CLI
369def build_parser() -> argparse.ArgumentParser:
370 parser = argparse.ArgumentParser(
371 description="Fit 1/2/3-compartment models with identifiability and model-selection diagnostics.",
372 formatter_class=argparse.RawDescriptionHelpFormatter,
373 epilog=__doc__,
374 )
375 parser.add_argument("-i", "--input", required=True)
376 parser.add_argument("--dose", type=float, required=True)
377 parser.add_argument("--route", choices=("iv-bolus", "iv-infusion", "oral"), default="iv-bolus")
378 parser.add_argument("--model", choices=sorted(MODELS), default="1cmt")
379 parser.add_argument("--compare", help="comma-separated models to compare, e.g. 1cmt,2cmt,3cmt")
380 parser.add_argument("--tinf", type=float, help="infusion duration for --route iv-infusion")
381 parser.add_argument("--interval", type=float, help="dosing interval for a multiple-dose profile")
382 parser.add_argument("--n-doses", type=int, default=1)
383 parser.add_argument("--absorption", choices=("first-order", "transit"), default="first-order")
384 parser.add_argument("--lag", action="store_true", help="estimate an absorption lag time")
385 parser.add_argument("--weight", choices=sorted(WEIGHTS), default="1/y2")
386 parser.add_argument("--starts", type=int, default=5, help="random restarts (default: 5)")
387 parser.add_argument("--time-column", default="time")
388 parser.add_argument("--conc-column", default="conc")
389 parser.add_argument("--subject-column", default="id")
390 parser.add_argument("--max-rse", type=float, default=50.0, help="flag %%RSE above this (default: 50)")
391 parser.add_argument("--max-corr", type=float, default=0.95, help="flag |correlation| above this (default: 0.95)")
392 parser.add_argument("--predictions", action="store_true", help="also emit the observed/predicted table")
393 add_format_argument(parser)
394 return parser
397def _fit_and_report(spec: FitSpec, args, subject: str, time, conc, report: Report) -> FitResult:
398 result = fit_one(spec, time, conc, args.weight, starts=args.starts)
399 label = f"subject {subject}" if subject else "fit"
401 rows = []
402 crit = t_dist.ppf(0.975, max(result.n - result.p, 1))
403 for i, name in enumerate(result.names):
404 se = result.se_log[i]
405 rows.append(
406 {
407 "parameter": name,
408 "estimate": result.theta[i],
409 "rse_pct": 100.0 * se if math.isfinite(se) else float("nan"),
410 "ci95_low": result.theta[i] * math.exp(-crit * se) if math.isfinite(se) else float("nan"),
411 "ci95_high": result.theta[i] * float(np.exp(min(crit * se, 709.0))) if math.isfinite(se) else float("nan"),
412 }
413 )
414 if math.isfinite(se) and 100.0 * se > args.max_rse:
415 report.finding(
416 f"{label}: {name} has {100 * se:.0f}% RSE - not estimable from these data at this model size"
417 )
418 report.table(f"{label}: {spec.model} {spec.route} parameters", rows)
420 secondary = secondary_parameters(spec, result.theta)
421 report.table(f"{label}: secondary parameters", [{"parameter": k, "value": v} for k, v in secondary.items()])
422 if secondary.get("flip_flop_suspected"):
423 report.finding(
424 f"{label}: ka is slower than the terminal disposition rate - the fit is in the flip-flop "
425 "branch, where ka and k are numerically exchangeable. Without IV data the assignment is a "
426 "modelling assumption, not an estimate."
427 )
429 report.table(
430 f"{label}: fit statistics",
431 [
432 {
433 "statistic": "observations",
434 "value": result.n,
435 },
436 {"statistic": "parameters", "value": result.p},
437 {"statistic": "weighted SSR", "value": result.wssr},
438 {"statistic": "residual SD", "value": result.sigma},
439 {"statistic": "AIC", "value": result.aic},
440 {"statistic": "BIC", "value": result.bic},
441 {"statistic": "condition number", "value": result.condition},
442 {"statistic": "runs-test p (residual signs)", "value": runs_test_p(result.weighted_residuals)},
443 {"statistic": "converged", "value": result.success},
444 ],
445 )
447 with np.errstate(invalid="ignore"):
448 for i in range(result.p):
449 for j in range(i + 1, result.p):
450 rho = result.correlation[i, j]
451 if math.isfinite(rho) and abs(rho) > args.max_corr:
452 report.finding(
453 f"{label}: {result.names[i]} and {result.names[j]} are correlated at "
454 f"{rho:+.3f} - the data cannot separate them; consider a simpler model"
455 )
456 if not math.isfinite(result.condition) or result.condition > 1000:
457 report.finding(
458 f"{label}: condition number {result.condition:.0f} exceeds 1000 - the model is "
459 "ill-conditioned or rank deficient; standard errors do not establish identifiability"
460 )
461 runs_p = runs_test_p(result.weighted_residuals)
462 if math.isfinite(runs_p) and runs_p < 0.05:
463 report.finding(
464 f"{label}: residual signs are not random (runs test p = {runs_p:.4f}) - a structural "
465 "pattern needing review of structure, time ordering, serial dependence and error model"
466 )
467 if not result.success:
468 report.finding(f"{label}: optimiser did not report convergence ({result.message})")
470 if args.predictions:
471 report.table(
472 f"{label}: observed vs predicted",
473 [
474 {
475 "time": float(t),
476 "observed": float(o),
477 "predicted": float(p),
478 "residual": float(r),
479 "weighted_residual": float(w),
480 }
481 for t, o, p, r, w in zip(time, conc, result.predicted, result.residuals, result.weighted_residuals)
482 ],
483 )
484 return result
487def run(argv: Sequence[str] | None = None) -> int:
488 args = build_parser().parse_args(argv)
489 from _common import validate_numeric_args
490 validate_numeric_args(args)
491 if args.weight in {"1/yhat", "1/yhat2"}:
492 raise InputError("prediction-dependent WLS omits the likelihood variance term; use a full likelihood/error-model tool")
493 if args.dose <= 0 or args.n_doses < 1 or (args.interval is not None and args.interval <= 0):
494 raise InputError("dose, dose count and interval must be positive")
495 if args.route == "iv-infusion" and (args.tinf is None or args.tinf <= 0):
496 print("error: --route iv-infusion needs --tinf", file=sys.stderr)
497 return EXIT_INPUT
499 rows = read_table(args.input)
500 require_columns(rows, [args.time_column, args.conc_column], str(args.input))
501 subjects = group_by(rows, args.subject_column) if args.subject_column in rows[0] else {"": rows}
503 report = Report()
504 report.note(f"weighting: {args.weight} ({WEIGHTS[args.weight]})")
505 report.note("parameters estimated on the log scale; intervals are local Gauss-Newton approximations, not profile likelihood intervals")
506 report.note("Observed-concentration weights can bias estimates; compare error models. AIC/BIC here are weighted-SSR scores, conditional on fixed weights. Compartment F p-values are exploratory: boundary/nonidentifiability invalidate a general nominal F-test.")
507 if args.route == "oral":
508 report.note(
509 "extravascular data alone identify CL/F and V/F, never CL and V separately. "
510 "Bioavailability requires an intravenous reference."
511 )
513 models = [m.strip() for m in args.compare.split(",")] if args.compare else [args.model]
514 for name in models:
515 if name not in MODELS:
516 raise InputError(f"unknown model {name!r}; choose from {', '.join(sorted(MODELS))}")
518 for subject, subject_rows in subjects.items():
519 time = np.asarray([parse_float(r[args.time_column], "time") for r in subject_rows], dtype=float)
520 conc = np.asarray([parse_float(r[args.conc_column], "conc") for r in subject_rows], dtype=float)
521 order = np.argsort(time)
522 time, conc = time[order], conc[order]
523 keep = conc > 0
524 if np.sum(keep) < len(conc):
525 report.note(f"subject {subject or '1'}: dropped {int(np.sum(~keep))} non-positive concentration(s)")
526 time, conc = time[keep], conc[keep]
528 fits = {}
529 for name in models:
530 spec = FitSpec(
531 model=name,
532 route=args.route,
533 dose=args.dose,
534 tinf=args.tinf,
535 interval=args.interval,
536 n_doses=args.n_doses,
537 absorption=args.absorption,
538 estimate_lag=args.lag,
539 )
540 fits[name] = _fit_and_report(spec, args, subject, time, conc, report)
542 if len(models) > 1:
543 comparison = []
544 ordered = sorted(fits.items(), key=lambda kv: kv[1].p)
545 for index, (name, fit) in enumerate(ordered):
546 entry = {"model": name, "parameters": fit.p, "wssr": fit.wssr, "aic": fit.aic, "bic": fit.bic}
547 if index > 0:
548 simpler_name, simpler = ordered[index - 1]
549 d_p = fit.p - simpler.p
550 if d_p > 0 and fit.wssr > 0 and fit.n - fit.p > 0:
551 f_stat = ((simpler.wssr - fit.wssr) / d_p) / (fit.wssr / (fit.n - fit.p))
552 entry["f_vs_simpler"] = f_stat
553 entry["f_p_value"] = float(1.0 - f_dist.cdf(f_stat, d_p, fit.n - fit.p)) if f_stat > 0 else 1.0
554 entry["compared_with"] = simpler_name
555 comparison.append(entry)
556 report.table(f"subject {subject or '1'}: model comparison", comparison)
557 best = min(fits.items(), key=lambda kv: kv[1].aic)[0]
558 report.note(
559 f"subject {subject or '1'}: lowest AIC is {best}. AIC and BIC are comparable here only "
560 "because every candidate was fitted to the same observations with the same weighting."
561 )
563 return report.emit(args.format)
566if __name__ == "__main__":
567 raise SystemExit(main_wrapper(run))