Skills › Scientific Agent Skills › pkpd-modeling › check_popk_dataset.py
check_popk_dataset.py
418 lines · skills/pkpd-modeling/scripts/check_popk_dataset.py at 5f2bc03
1#!/usr/bin/env python3
2"""Check common NONMEM-style population PK dataset conventions.
4This is a partial sanity checker, not an NM-TRAN emulator. Verify actual event
5ordering, numeric/censored DV conventions, dose history and the model's input
6mapping. Advanced MDV/SS, pure PD and multi-analyte datasets need tailored checks.
8 python3 check_popk_dataset.py -i nmdata.csv
9 python3 check_popk_dataset.py -i nmdata.csv --covariates WT,CRCL --time-varying WT
11Exit 1 means a finding at or above --fail-on; --strict also fails on warnings.
12"""
14from __future__ import annotations
16import argparse
17import math
18from collections import Counter
19from typing import Sequence
21from _common import (
22 Report,
23 add_format_argument,
24 main_wrapper,
25 read_table,
26)
28SEVERITY_ORDER = {"note": 0, "warning": 1, "error": 2}
30# NONMEM data item conventions. EVID 3 resets the system, 4 resets and doses.
31EVID_MEANINGS = {
32 "0": "observation",
33 "1": "dose",
34 "2": "other-type event (no dose, no observation)",
35 "3": "reset",
36 "4": "reset and dose",
37}
39MISSING = {"", ".", "na", "n/a", "nan", "null", "none"}
42def _num(value: str) -> float | None:
43 text = (value or "").strip()
44 if text.lower() in MISSING:
45 return None
46 try:
47 value = float(text)
48 return value if math.isfinite(value) else None
49 except ValueError:
50 return None
53def _is_missing(value: str) -> bool:
54 return (value or "").strip().lower() in MISSING
57class Findings:
58 def __init__(self) -> None:
59 self.items: list[dict[str, object]] = []
61 def add(self, severity: str, check: str, detail: str, rows: Sequence[int] = ()) -> None:
62 listed = list(rows)[:8]
63 self.items.append(
64 {
65 "severity": severity,
66 "check": check,
67 "detail": detail,
68 "n_records": len(rows),
69 "example_rows": ", ".join(str(r) for r in listed) + ("..." if len(rows) > 8 else ""),
70 }
71 )
73 def max_severity(self) -> int:
74 return max((SEVERITY_ORDER[i["severity"]] for i in self.items), default=-1)
77# ------------------------------------------------------------------ checks
80def check_dataset(rows: list[dict[str, str]], args: argparse.Namespace) -> tuple[Findings, dict]:
81 found = Findings()
82 columns = list(rows[0].keys())
83 upper = {c.lower(): c for c in columns}
85 def col(name: str) -> str | None:
86 return upper.get(name.lower())
88 id_col = col(args.id_column)
89 time_col = col(args.time_column)
90 dv_col = col(args.dv_column)
91 amt_col = col(args.amt_column)
92 evid_col, mdv_col = col("evid"), col("mdv")
93 rate_col, ss_col, ii_col, addl_col = col("rate"), col("ss"), col("ii"), col("addl")
95 for required, name in ((id_col, args.id_column), (time_col, args.time_column), (dv_col, args.dv_column)):
96 if required is None:
97 found.add("error", "required column", f"no '{name}' column; NONMEM cannot build a dataset without it")
98 if amt_col is None:
99 found.add("error", "required column", "no AMT column: without it no dose is ever administered")
100 if evid_col is None:
101 found.add(
102 "warning",
103 "EVID absent",
104 "no EVID column. NONMEM then infers events from AMT, which works only if every dose "
105 "record has AMT>0 and every observation has AMT=0 or missing. Supplying EVID is safer.",
106 )
107 if not id_col or not time_col:
108 return found, {}
110 # ---- column-name hazards
111 long_names = [c for c in columns if len(c) > 20]
112 if long_names:
113 found.add("note", "column names", f"unusually long column name(s): {', '.join(long_names[:4])}")
114 non_alnum = [c for c in columns if not c.replace("_", "").isalnum()]
115 if non_alnum:
116 found.add(
117 "warning",
118 "column names",
119 f"column name(s) with characters NM-TRAN will reject in $INPUT: {', '.join(non_alnum[:6])}",
120 )
122 # ---- per-record checks
123 bad_dv_text, neg_time, neg_conc, bad_evid = [], [], [], []
124 dose_without_amt, obs_with_amt, obs_missing_dv, mdv_conflict = [], [], [], []
125 rate_issues, ss_missing_ii, addl_missing_ii, zero_ii = [], [], [], []
127 subjects: dict[str, list[tuple[int, dict[str, str]]]] = {}
128 for index, row in enumerate(rows, start=2): # header is line 1
129 subject = (row.get(id_col) or "").strip()
130 subjects.setdefault(subject, []).append((index, row))
132 evid = (row.get(evid_col) or "0").strip() if evid_col else None
133 if evid_col and evid not in EVID_MEANINGS and not _is_missing(evid):
134 bad_evid.append(index)
136 time = _num(row.get(time_col, ""))
137 if time is None:
138 found.add("error", "TIME not numeric", f"row {index}: TIME is {row.get(time_col)!r}", [index])
139 elif time < 0 and not args.allow_negative_time:
140 neg_time.append(index)
142 dv_raw = row.get(dv_col, "") if dv_col else ""
143 amt = _num(row.get(amt_col, "")) if amt_col else None
144 is_dose = (evid in {"1", "4"}) if evid_col else bool(amt and amt > 0)
145 is_obs = (evid == "0") if evid_col else not is_dose
147 if is_obs:
148 if not _is_missing(dv_raw) and _num(dv_raw) is None:
149 bad_dv_text.append(index)
150 value = _num(dv_raw)
151 if value is not None and value < 0:
152 neg_conc.append(index)
153 mdv = (row.get(mdv_col) or "0").strip() if mdv_col else "0"
154 if _is_missing(dv_raw) and mdv != "1":
155 obs_missing_dv.append(index)
156 if not _is_missing(dv_raw) and mdv == "1":
157 mdv_conflict.append(index)
158 if amt not in (None, 0.0):
159 obs_with_amt.append(index)
160 if is_dose and (amt is None or amt <= 0) and not (ss_col and (_num(row.get(ss_col, "")) or 0) > 0 and rate_col and (_num(row.get(rate_col, "")) or 0) > 0):
161 dose_without_amt.append(index)
163 if rate_col and not _is_missing(row.get(rate_col, "")):
164 rate = _num(row.get(rate_col, ""))
165 if rate is None or (rate < 0 and rate not in (-1.0, -2.0)):
166 rate_issues.append(index)
167 if ss_col and (row.get(ss_col) or "0").strip() not in {"0", "", "."}:
168 ii = _num(row.get(ii_col, "")) if ii_col else None
169 constant_infusion = amt in (None, 0.0) and (_num(row.get(rate_col, "")) or 0) > 0
170 if (ii is None or ii <= 0) and not constant_infusion:
171 ss_missing_ii.append(index)
172 if addl_col and (_num(row.get(addl_col, "")) or 0) > 0:
173 ii = _num(row.get(ii_col, "")) if ii_col else None
174 if ii is None:
175 addl_missing_ii.append(index)
176 elif ii <= 0:
177 zero_ii.append(index)
179 if bad_dv_text:
180 found.add(
181 "error",
182 "non-numeric DV",
183 "DV contains nonnumeric or nonfinite text (for example BLQ). Do not rely on parsing or "
184 "automatic recoding to represent censoring. Recode to a numeric "
185 "value plus a BLQ flag column and choose an M-method.",
186 bad_dv_text,
187 )
188 if neg_time:
189 found.add("error", "negative TIME", "TIME before zero without --allow-negative-time", neg_time)
190 if neg_conc:
191 found.add("warning", "negative DV", "negative concentrations in observation records", neg_conc)
192 if bad_evid:
193 found.add("error", "invalid EVID", f"EVID outside {sorted(EVID_MEANINGS)}", bad_evid)
194 if dose_without_amt:
195 found.add("error", "dose without AMT", "EVID marks a dose but AMT is missing or not positive", dose_without_amt)
196 if obs_with_amt:
197 found.add(
198 "warning",
199 "observation carries AMT",
200 "records marked as observations also carry a non-zero AMT; if EVID is dropped from $INPUT "
201 "these become doses",
202 obs_with_amt,
203 )
204 if obs_missing_dv:
205 found.add("error", "missing DV with MDV=0", "observation record has no DV but is not flagged MDV=1", obs_missing_dv)
206 if mdv_conflict:
207 found.add("warning", "MDV=1 with a DV present", "the observation will be ignored by the estimation", mdv_conflict)
208 if rate_issues:
209 found.add("error", "invalid RATE", "RATE must be 0 (bolus), >0 (known infusion rate), -1 (modelled rate), or -2 (modelled duration)", rate_issues)
210 if ss_missing_ii:
211 found.add("error", "SS without II", "a steady-state record needs a positive II", ss_missing_ii)
212 if addl_missing_ii:
213 found.add("error", "ADDL without II", "ADDL needs positive II; invalid events must be corrected before NM-TRAN/PREDPP", addl_missing_ii)
214 if zero_ii:
215 found.add("error", "ADDL with II<=0", "II must be positive for ADDL to place additional doses", zero_ii)
217 seen_ids, previous_id = set(), None
218 for row in rows:
219 sid = row.get(id_col, "")
220 if sid != previous_id and sid in seen_ids:
221 found.add("error", "noncontiguous ID", f"subject {sid} occurs in separate blocks")
222 seen_ids.add(sid)
223 previous_id = sid
224 if _num(sid) is None:
225 found.add("error", "ID not numeric", "NONMEM ID must be numeric")
227 # ---- per-subject checks
228 no_obs, no_dose, unsorted, duplicated, first_not_dose = [], [], [], [], []
229 obs_counts, dose_counts = [], []
230 for subject, records in subjects.items():
231 if subject == "":
232 found.add("error", "blank ID", "records with an empty ID column", [r[0] for r in records])
233 continue
234 times = []
235 n_obs = n_dose = 0
236 for index, row in records:
237 time = _num(row.get(time_col, ""))
238 times.append((time, index))
239 evid = (row.get(evid_col) or "0").strip() if evid_col else None
240 amt = _num(row.get(amt_col, "")) if amt_col else None
241 if (evid in {"1", "4"}) if evid_col else bool(amt and amt > 0):
242 n_dose += 1
243 elif (evid == "0") if evid_col else True:
244 if not _is_missing(row.get(dv_col, "")) and (not mdv_col or _num(row.get(mdv_col, "0")) in (0, None)):
245 n_obs += 1
246 obs_counts.append(n_obs)
247 dose_counts.append(n_dose)
248 if n_obs == 0:
249 no_obs.append(subject)
250 if n_dose == 0:
251 no_dose.append(subject)
252 numeric = [t for t, _ in times if t is not None]
253 previous_time = None
254 decreased = False
255 for _, row in records:
256 current_time = _num(row.get(time_col, ""))
257 reset = evid_col and row.get(evid_col) in {"3", "4"}
258 if current_time is not None and previous_time is not None and current_time < previous_time and not reset:
259 decreased = True
260 if current_time is not None:
261 previous_time = current_time
262 if decreased:
263 unsorted.append(subject)
264 counts = Counter(numeric)
265 if any(v > 1 for v in counts.values()):
266 duplicated.append(subject)
267 if records:
268 first_evid = (records[0][1].get(evid_col) or "0").strip() if evid_col else None
269 first_amt = _num(records[0][1].get(amt_col, "")) if amt_col else None
270 is_dose_first = (first_evid in {"1", "4"}) if evid_col else bool(first_amt and first_amt > 0)
271 if not is_dose_first:
272 first_not_dose.append(subject)
274 if no_obs:
275 found.add(
276 "error",
277 "subject with no observations",
278 f"{len(no_obs)} subject(s) contribute no DV: {', '.join(no_obs[:6])}. Their ETAs are drawn "
279 "entirely from the prior and they inflate the apparent N of the analysis.",
280 )
281 if no_dose:
282 found.add("error", "subject with no dose", f"{len(no_dose)} subject(s) have observations but no dose: {', '.join(no_dose[:6])}")
283 if unsorted:
284 found.add("error", "TIME not sorted", f"{len(unsorted)} subject(s) have out-of-order TIME: {', '.join(unsorted[:6])}")
285 if duplicated:
286 found.add(
287 "warning",
288 "duplicate TIME within a subject",
289 f"{len(duplicated)} subject(s) have records sharing a timestamp: {', '.join(duplicated[:6])}. "
290 "NONMEM applies them in file order, so a dose and an observation at the same time give a "
291 "pre-dose or post-dose value depending purely on which row came first.",
292 )
293 if first_not_dose:
294 found.add(
295 "note",
296 "first record is not a dose",
297 f"{len(first_not_dose)} subject(s) start with a non-dose record: {', '.join(first_not_dose[:6])}. "
298 "Fine for a pre-dose baseline; wrong if a dose is missing.",
299 )
301 # ---- covariates
302 for name in args.covariates:
303 cov_col = col(name)
304 if cov_col is None:
305 found.add("error", "covariate absent", f"--covariates names '{name}' but there is no such column")
306 continue
307 missing_rows, non_numeric = [], []
308 per_subject_values: dict[str, set[str]] = {}
309 for index, row in enumerate(rows, start=2):
310 raw = row.get(cov_col, "")
311 subject = (row.get(id_col) or "").strip()
312 per_subject_values.setdefault(subject, set()).add(raw.strip())
313 if _is_missing(raw):
314 missing_rows.append(index)
315 elif _num(raw) is None:
316 non_numeric.append(index)
317 if missing_rows:
318 found.add(
319 "error",
320 f"covariate {name} missing",
321 f"{len(missing_rows)} record(s) have no value. A blank or '.' is read as 0, not as "
322 f"missing, so those records enter the covariate model with {name} = 0 - a 0 kg "
323 "body weight, a 0 mL/min creatinine clearance - rather than being excluded.",
324 missing_rows,
325 )
326 if non_numeric:
327 found.add("error", f"covariate {name} non-numeric", "values NM-TRAN cannot read as numbers", non_numeric)
328 varying = [s for s, values in per_subject_values.items() if len(values) > 1]
329 if varying and name.lower() not in {t.lower() for t in args.time_varying}:
330 found.add(
331 "warning",
332 f"covariate {name} varies within subject",
333 f"{len(varying)} subject(s) have more than one value: {', '.join(varying[:6])}. "
334 "If this is intended, list it in --time-varying; if not, the model will use whichever "
335 "value is on the record being evaluated.",
336 )
338 summary = {
339 "records": len(rows),
340 "subjects": len(subjects),
341 "observations": sum(obs_counts),
342 "dose_records": sum(dose_counts),
343 "median_observations_per_subject": sorted(obs_counts)[len(obs_counts) // 2] if obs_counts else 0,
344 "subjects_with_one_observation": sum(1 for c in obs_counts if c == 1),
345 "columns": len(columns),
346 }
347 if summary["subjects"] and summary["observations"] / max(summary["subjects"], 1) < 2:
348 found.add(
349 "note",
350 "sparse sampling",
351 "fewer than two observations per subject on average. A structural model with an absorption "
352 "phase is unlikely to be identifiable; consider fixing parameters from a richer study.",
353 )
354 return found, summary
357# --------------------------------------------------------------------- CLI
360def build_parser() -> argparse.ArgumentParser:
361 parser = argparse.ArgumentParser(
362 description="Validate a population PK dataset against NONMEM/nlmixr2 data conventions.",
363 formatter_class=argparse.RawDescriptionHelpFormatter,
364 epilog=__doc__,
365 )
366 parser.add_argument("-i", "--input", required=True)
367 parser.add_argument("--id-column", default="ID")
368 parser.add_argument("--time-column", default="TIME")
369 parser.add_argument("--dv-column", default="DV")
370 parser.add_argument("--amt-column", default="AMT")
371 parser.add_argument("--covariates", default="", help="comma-separated covariate columns to check")
372 parser.add_argument("--time-varying", default="", help="covariates that are expected to change within a subject")
373 parser.add_argument("--allow-negative-time", action="store_true")
374 parser.add_argument(
375 "--fail-on",
376 choices=("note", "warning", "error"),
377 default="error",
378 help="lowest severity that makes the exit code 1 (default: error)",
379 )
380 parser.add_argument("--strict", action="store_true", help="shorthand for --fail-on warning")
381 add_format_argument(parser)
382 return parser
385def run(argv: Sequence[str] | None = None) -> int:
386 args = build_parser().parse_args(argv)
387 from _common import validate_numeric_args
388 validate_numeric_args(args)
389 args.covariates = [c.strip() for c in args.covariates.split(",") if c.strip()]
390 args.time_varying = [c.strip() for c in args.time_varying.split(",") if c.strip()]
391 if args.strict:
392 args.fail_on = "warning"
394 rows = read_table(args.input)
395 found, summary = check_dataset(rows, args)
397 report = Report()
398 for key, value in summary.items():
399 report.scalar(key, value)
400 if found.items:
401 ordered = sorted(found.items, key=lambda i: -SEVERITY_ORDER[str(i["severity"])])
402 report.table("dataset findings", ordered)
403 else:
404 report.table("dataset findings", [{"severity": "none", "check": "all checks passed", "detail": "", "n_records": 0, "example_rows": ""}])
406 report.note(
407 "This checks data conventions, not pharmacology. A dataset that passes can still be wrong: "
408 "verify units, time origin, and that DV is the quantity the model predicts."
409 )
410 threshold = SEVERITY_ORDER[args.fail_on]
411 for item in found.items:
412 if SEVERITY_ORDER[str(item["severity"])] >= threshold:
413 report.finding(f"[{item['severity']}] {item['check']}: {item['detail']}")
414 return report.emit(args.format)
417if __name__ == "__main__":
418 raise SystemExit(main_wrapper(run))