Skills › Scientific Agent Skills › analytical-method-validation › check_accuracy_precision.py
check_accuracy_precision.py
359 lines · skills/analytical-method-validation/scripts/check_accuracy_precision.py at 5f2bc03
1#!/usr/bin/env python3
2"""Evaluate accuracy and precision per ICH Q2(R2) 3.3.
4Accuracy is reported as mean recovery with a confidence interval, which is what
5Q2(R2) 3.3.1.4 asks for -- a bare mean is not sufficient. Precision is
6decomposed into repeatability and intermediate precision by a one-way
7random-effects model, because pooling all results into a single standard
8deviation understates the day-to-day variability the procedure will actually
9show in routine use.
11 python3 check_accuracy_precision.py --input ap.csv
12 python3 check_accuracy_precision.py -i ap.csv --accuracy-limit 2 --rsd-limit 2
13 python3 check_accuracy_precision.py -i ap.csv --design-check assay --format json
15Input columns:
16 level nominal / added concentration (groups the accuracy analysis)
17 measured measured or recovered concentration
18 group optional: day, analyst, instrument or run -- the intermediate
19 precision factor. Without it only repeatability is estimated.
21Exit codes: 0 no findings, 1 findings raised, 2 bad input.
22"""
24from __future__ import annotations
26import argparse
27import math
28import sys
30sys.path.insert(0, str(__import__("pathlib").Path(__file__).resolve().parent))
32from _common import ( # noqa: E402
33 EXIT_FINDINGS,
34 EXIT_OK,
35 InputError,
36 add_common_args,
37 emit,
38 finding,
39 finite_float,
40 mean,
41 note,
42 one_way_components,
43 parse_rows,
44 read_input,
45 require_columns,
46 rsd_percent,
47 run_cli,
48 sample_sd,
49 sd_confidence_interval,
50 t_ppf,
51 to_float,
52)
54# ICH Q2(R2) 3.3.2.1 minima, used only to comment on the design. Either option
55# is sufficient on its own.
56DESIGN_MINIMA = {
57 "assay": {
58 "range_determinations": 9,
59 "range_levels": 3,
60 "single_level_determinations": 6,
61 "note": "9 determinations across the range (3x3), or 6 at 100% of test concentration",
62 },
63 "impurity": {
64 "range_determinations": 9,
65 "range_levels": 3,
66 "single_level_determinations": 6,
67 "note": "9 determinations across the range (3x3), or 6 at 100% of test concentration",
68 },
69}
72def main() -> int:
73 parser = argparse.ArgumentParser(
74 description="Check accuracy and precision from validation data."
75 )
76 parser.add_argument("--input", "-i", help="CSV/TSV/JSON file, or '-' for stdin")
77 parser.add_argument("--accuracy-limit", type=finite_float, default=None,
78 help="flag a level whose mean recovery deviates more than this %% "
79 "from 100%%")
80 parser.add_argument("--rsd-limit", type=finite_float, default=None,
81 help="flag repeatability or intermediate precision %%RSD above this")
82 parser.add_argument("--ci-level", type=finite_float, default=0.95,
83 help="confidence level for accuracy intervals (default 0.95)")
84 parser.add_argument("--sd-ci-level", type=finite_float, default=0.90,
85 help="confidence level for SD intervals (default 0.90)")
86 parser.add_argument("--design-check", choices=sorted(DESIGN_MINIMA), default=None,
87 help="comment on the design against Q2(R2) recommended minima")
88 parser.add_argument("--test-concentration", type=finite_float,
89 help="nominal level representing 100%% of test concentration for option (b)")
90 parser.add_argument("--require-ci-within-limit", action="store_true",
91 help="require the whole accuracy CI inside the limit, not just the mean")
92 add_common_args(parser)
93 args = parser.parse_args()
95 if not 0.5 <= args.ci_level < 1.0:
96 raise InputError("--ci-level must be in [0.5, 1)")
98 if not 0.5 <= args.sd_ci_level < 1.0:
99 raise InputError("--sd-ci-level must be in [0.5, 1)")
100 for name in ("accuracy_limit", "rsd_limit", "test_concentration"):
101 value = getattr(args, name)
102 if value is not None and (not math.isfinite(value) or value <= 0):
103 raise InputError(f"--{name.replace('_', '-')} must be finite and > 0")
104 if args.require_ci_within_limit and args.accuracy_limit is None:
105 raise InputError("--require-ci-within-limit needs --accuracy-limit")
107 rows = parse_rows(read_input(args.input), args.input)
108 require_columns(rows, ["level", "measured"])
109 has_group = any("group" in row for row in rows)
110 if has_group and any(not (r.get("group") or "").strip() for r in rows):
111 raise InputError("group must be present and non-empty on every row")
113 records = []
114 for i, r in enumerate(rows):
115 records.append(
116 {
117 "level": to_float(r["level"], "level", i),
118 "measured": to_float(r["measured"], "measured", i),
119 "group": (r.get("group") or "all").strip() or "all",
120 }
121 )
123 findings: list[str] = []
125 # ---------------- Accuracy ----------------
126 levels = sorted({rec["level"] for rec in records})
127 accuracy_rows = []
128 all_recoveries: list[float] = []
129 for level in levels:
130 vals = [rec["measured"] for rec in records if rec["level"] == level]
131 if level == 0:
132 raise InputError("a nominal level of 0 cannot be used for recovery")
133 recoveries = [100.0 * v / level for v in vals]
134 all_recoveries.extend(recoveries)
135 independent = recoveries
136 interval_basis = "individual independent preparations"
137 if has_group:
138 accuracy_groups: dict[str, list[float]] = {}
139 for rec in records:
140 if rec["level"] == level:
141 accuracy_groups.setdefault(rec["group"], []).append(100.0 * rec["measured"] / level)
142 if len(accuracy_groups) > 1:
143 independent = [mean(values) for values in accuracy_groups.values()]
144 interval_basis = "equally weighted independent group means"
145 m = mean(independent)
146 n = len(recoveries)
147 n_independent = len(independent)
148 sd = sample_sd(recoveries)
149 if n_independent >= 2:
150 half = t_ppf(0.5 + args.ci_level / 2.0, n_independent - 1) * sample_sd(independent) / math.sqrt(n_independent)
151 lo, hi = m - half, m + half
152 else:
153 sd, lo, hi = float("nan"), float("nan"), float("nan")
154 accuracy_rows.append(
155 {
156 "level": level,
157 "n": n,
158 "n_independent_for_ci": n_independent,
159 "interval_basis": interval_basis,
160 "mean_measured": m * level / 100.0,
161 "mean_recovery_pct": m,
162 "bias_pct": m - 100.0,
163 "sd_recovery_pct": sd,
164 f"ci{int(args.ci_level * 100)}_low": lo,
165 f"ci{int(args.ci_level * 100)}_high": hi,
166 }
167 )
168 if args.accuracy_limit is not None:
169 if abs(m - 100.0) > args.accuracy_limit:
170 findings.append(
171 f"level {level:g}: mean recovery {m:.2f}% is {m - 100.0:+.2f}% from nominal, "
172 f"outside +/-{args.accuracy_limit:g}%"
173 )
174 elif args.require_ci_within_limit:
175 if not math.isfinite(lo) or lo < 100.0 - args.accuracy_limit or hi > 100.0 + args.accuracy_limit:
176 findings.append(
177 f"level {level:g}: mean recovery {m:.2f}% is inside "
178 f"+/-{args.accuracy_limit:g}% but its "
179 f"{int(args.ci_level * 100)}% CI ({lo:.2f}, {hi:.2f}) is not -- the data "
180 "do not demonstrate accuracy at this limit"
181 )
183 # ---------------- Precision ----------------
184 # Precision is estimated WITHIN each concentration level. Pooling levels
185 # together would let the 80/100/120 spread masquerade as imprecision.
186 pct = int(args.sd_ci_level * 100)
187 precision_rows: list[dict] = []
189 def precision_block(label: str, groups: dict[str, list[float]]):
190 """Append rows for one concentration level."""
191 try:
192 comp = one_way_components(groups)
193 except InputError as exc:
194 vals = [v for vs in groups.values() for v in vs]
195 if len(vals) < 2:
196 note(f"{label}: too few values for a precision estimate")
197 if args.rsd_limit is not None:
198 findings.append(f"{label}: too few determinations to assess precision")
199 return None
200 if len(groups) > 1:
201 findings.append(f"{label}: precision components unidentifiable ({exc})")
202 return None
203 sd = sample_sd(vals)
204 base = abs(mean(vals))
205 rsd = 100.0 * sd / base if base else float("nan")
206 lo, hi = sd_confidence_interval(sd, len(vals) - 1, args.sd_ci_level)
207 precision_rows.append(
208 {"level": label, "component": "repeatability only", "sd": sd, "rsd_pct": rsd,
209 "df": len(vals) - 1, f"ci{pct}_low_sd": lo, f"ci{pct}_high_sd": hi}
210 )
211 note(f"{label}: intermediate precision not estimated ({exc})")
212 if args.rsd_limit is not None and (not math.isfinite(rsd) or rsd > args.rsd_limit):
213 findings.append(f"{label}: repeatability RSD exceeds the limit or is undefined")
214 return None
216 base = abs(comp.grand_mean)
218 def as_rsd(sd: float) -> float:
219 return 100.0 * sd / base if base else float("nan")
221 sr_lo, sr_hi = sd_confidence_interval(
222 comp.sd_repeatability, comp.df_within, args.sd_ci_level
223 )
224 sat_df = comp.satterthwaite_df()
225 si_lo, si_hi = sd_confidence_interval(comp.sd_intermediate, sat_df, args.sd_ci_level)
226 precision_rows.extend(
227 [
228 {"level": label, "component": "repeatability (within group)",
229 "sd": comp.sd_repeatability, "rsd_pct": as_rsd(comp.sd_repeatability),
230 "df": comp.df_within, f"ci{pct}_low_sd": sr_lo, f"ci{pct}_high_sd": sr_hi},
231 {"level": label, "component": "between-group", "sd": comp.sd_between,
232 "rsd_pct": as_rsd(comp.sd_between), "df": comp.df_between,
233 f"ci{pct}_low_sd": float("nan"), f"ci{pct}_high_sd": float("nan")},
234 {"level": label, "component": "intermediate precision (total)",
235 "sd": comp.sd_intermediate, "rsd_pct": as_rsd(comp.sd_intermediate),
236 "df": sat_df, f"ci{pct}_low_sd": si_lo, f"ci{pct}_high_sd": si_hi},
237 ]
238 )
239 if comp.sd_between == 0.0:
240 note(
241 f"{label}: between-group variance estimated as zero (MS_between <= MS_within); "
242 "the groups are indistinguishable at this level"
243 )
244 if not comp.balanced:
245 note(f"{label}: unbalanced design, effective group size {comp.n_effective:.3f}")
246 if args.rsd_limit is not None:
247 for name, value in (
248 ("repeatability", as_rsd(comp.sd_repeatability)),
249 ("intermediate precision", as_rsd(comp.sd_intermediate)),
250 ):
251 if not math.isfinite(value) or value > args.rsd_limit:
252 findings.append(
253 f"{label}: {name} {value:.3f}% RSD exceeds the stated "
254 f"{args.rsd_limit:g}% limit"
255 )
256 return comp
258 if has_group:
259 for level in levels:
260 groups: dict[str, list[float]] = {}
261 for rec in records:
262 if rec["level"] == level:
263 groups.setdefault(rec["group"], []).append(rec["measured"])
264 precision_block(f"{level:g}", groups)
266 # Normalized recoveries are not pooled: level-specific bias would inflate precision.
267 else:
268 for level in levels:
269 vals = [rec["measured"] for rec in records if rec["level"] == level]
270 if len(vals) < 2:
271 if args.rsd_limit is not None:
272 findings.append(f"level {level:g}: too few determinations to assess precision")
273 continue
274 sd = sample_sd(vals)
275 lo, hi = sd_confidence_interval(sd, len(vals) - 1, args.sd_ci_level)
276 rsd = 100.0 * sd / abs(mean(vals)) if mean(vals) else float("nan")
277 precision_rows.append(
278 {"level": f"{level:g}", "component": "repeatability (no group column)",
279 "sd": sd, "rsd_pct": rsd, "df": len(vals) - 1,
280 f"ci{pct}_low_sd": lo, f"ci{pct}_high_sd": hi}
281 )
282 if args.rsd_limit is not None and (not math.isfinite(rsd) or rsd > args.rsd_limit):
283 findings.append(
284 f"{level:g}: repeatability {rsd:.3f}% RSD exceeds {args.rsd_limit:g}%"
285 )
286 note(
287 "no `group` column, so only repeatability was estimated. Q2(R2) 3.3.2.2 expects "
288 "intermediate precision from different days, analysts or equipment"
289 )
291 naive = rsd_percent([rec["measured"] for rec in records])
292 if len(levels) > 1:
293 note(
294 f"a single SD over every result regardless of level would report {naive:.3f}% RSD, "
295 "which is a range effect and not precision. Precision is reported per level below"
296 )
298 # ---------------- Design commentary ----------------
299 if args.design_check:
300 note("design screen checks counts only; confirm independent preparations, coverage of "
301 "the reportable range, and repeatability conditions in the protocol")
302 # Q2(R2) 3.3.2.1 offers two alternatives, and either one is sufficient:
303 # (a) >=9 determinations covering the reportable range (e.g. 3 levels x 3), or
304 # (b) >=6 determinations at 100% of the test concentration.
305 # Flag only when neither holds. Requiring 9 unconditionally would raise a
306 # finding against a design the guideline explicitly permits.
307 spec = DESIGN_MINIMA[args.design_check]
308 option_a = len(records) >= spec["range_determinations"] and len(levels) >= spec["range_levels"]
309 per_level = {lv: sum(1 for r in records if r["level"] == lv) for lv in levels}
310 best_single = per_level.get(args.test_concentration, 0)
311 option_b = best_single >= spec["single_level_determinations"]
312 if option_a:
313 note(
314 f"count screen matches Q2(R2) 3.3.2.1 option (a): {len(records)} determinations "
315 f"across {len(levels)} levels"
316 )
317 elif option_b:
318 note(
319 f"count screen matches Q2(R2) 3.3.2.1 option (b): {best_single} determinations at "
320 "the declared 100% test concentration. Option (b) gives no information about precision "
321 "across the range"
322 )
323 else:
324 findings.append(
325 f"repeatability design meets neither Q2(R2) 3.3.2.1 option: "
326 f"{len(records)} determinations across {len(levels)} level(s), with at most "
327 f"{best_single} at the declared --test-concentration. Option (a) needs "
328 f"{spec['range_determinations']} across at least {spec['range_levels']} levels; "
329 f"option (b) needs {spec['single_level_determinations']} at 100% of the test "
330 "concentration"
331 )
333 if args.format == "json":
334 emit([{"accuracy": accuracy_rows, "precision": precision_rows,
335 "overall_mean_recovery_pct": mean(all_recoveries),
336 "findings": findings}], "json")
337 else:
338 emit(accuracy_rows, args.format)
339 print()
340 emit(precision_rows, args.format)
342 note(
343 "Q2(R2) 3.3.1.4: report accuracy as mean percent recovery, or the difference from the "
344 "accepted true value, with a 100(1-alpha)% confidence interval"
345 )
346 note(
347 "Q2(R2) 3.3.2.2: intermediate precision covers days, environmental conditions, analysts "
348 "and equipment; the effects need not be studied individually"
349 )
350 for f in findings:
351 finding(f)
352 if not findings:
353 note("no findings against the checks that were run")
354 note("this tool does not decide that accuracy or precision is acceptable")
355 return EXIT_FINDINGS if findings else EXIT_OK
358if __name__ == "__main__":
359 run_cli(main)