Skills › Scientific Agent Skills › analytical-method-validation › check_detection_limits.py
check_detection_limits.py
326 lines · skills/analytical-method-validation/scripts/check_detection_limits.py at 5f2bc03
1#!/usr/bin/env python3
2"""Estimate DL and QL with supported approaches from ICH Q2(R2) 3.2.3.
4Supported estimates can disagree because they use different response-variation
5data. Reporting one number without naming the approach is the finding an
6assessor raises, so this script computes all of the applicable ones side by side
7and checks the answer against the reporting threshold it has to serve.
9 # sigma from the calibration line, slope from the same fit
10 python3 check_detection_limits.py --calibration calib.csv
12 # sigma from blank responses
13 python3 check_detection_limits.py --calibration calib.csv --blanks blanks.csv
15 # confirm an estimated QL with real data at that level
16 python3 check_detection_limits.py --calibration calib.csv \
17 --confirm-ql 0.05 --confirm-data ql_check.csv --reporting-threshold 0.05 \
18 --confirm-accuracy-limit 10 --confirm-rsd-limit 10
20Input:
21 --calibration CSV with `level` and `response` (low-range calibration curve)
22 --blanks CSV with `response` (blank measurements)
23 --intercepts CSV with `intercept` (independent low-range calibration curves)
24 --confirm-data CSV with `measured` (results at or near the claimed QL)
26Exit codes: 0 no findings, 1 findings raised, 2 bad input.
27"""
29from __future__ import annotations
31import argparse
32import math
33import sys
35sys.path.insert(0, str(__import__("pathlib").Path(__file__).resolve().parent))
37from _catalog import DL_QL_APPROACHES # noqa: E402
38from _common import ( # noqa: E402
39 EXIT_FINDINGS,
40 EXIT_OK,
41 InputError,
42 add_common_args,
43 emit,
44 finding,
45 finite_float,
46 fit_linear,
47 mean,
48 note,
49 parse_rows,
50 read_input,
51 require_columns,
52 rsd_percent,
53 run_cli,
54 sample_sd,
55 to_float,
56)
58DL_FACTOR = 3.3
59QL_FACTOR = 10.0
62def main() -> int:
63 parser = argparse.ArgumentParser(
64 description="Estimate detection and quantitation limits per ICH Q2(R2) 3.2.3."
65 )
66 parser.add_argument("--calibration", required=True,
67 help="CSV/TSV/JSON with `level` and `response`")
68 parser.add_argument("--blanks", help="CSV/TSV/JSON with `response` for blank samples")
69 parser.add_argument("--intercepts", help="CSV/TSV/JSON with `intercept` from independent low-range calibration curves")
70 parser.add_argument("--signal-to-noise", type=finite_float, default=None,
71 help="measured S/N at a stated concentration (use with --sn-level)")
72 parser.add_argument("--sn-level", type=finite_float, default=None,
73 help="the concentration at which --signal-to-noise was measured")
74 parser.add_argument("--confirm-ql", type=finite_float, default=None,
75 help="the QL being claimed, to be confirmed with --confirm-data")
76 parser.add_argument("--confirm-data", help="CSV/TSV/JSON with `measured` at/near the QL")
77 parser.add_argument("--confirm-accuracy-limit", type=finite_float, default=None,
78 help="pre-stated max %% bias for QL confirmation; no guideline default")
79 parser.add_argument("--confirm-rsd-limit", type=finite_float, default=None,
80 help="pre-stated max %%RSD for QL confirmation; no guideline default")
81 parser.add_argument("--reporting-threshold", type=finite_float, default=None,
82 help="impurity reporting threshold the QL must be at or below")
83 parser.add_argument("--weight", choices=("none", "1/x", "1/x2"), default="none")
84 add_common_args(parser)
85 args = parser.parse_args()
86 for name in ("signal_to_noise", "sn_level", "confirm_ql", "reporting_threshold",
87 "confirm_accuracy_limit", "confirm_rsd_limit"):
88 value = getattr(args, name)
89 if value is not None and (not math.isfinite(value) or value <= 0):
90 raise InputError(f"--{name.replace('_', '-')} must be finite and > 0")
91 if args.weight != "none":
92 raise InputError("DL/QL sigma-and-slope estimates require an unweighted low-range fit; "
93 "weighted residual scale is not a response SD")
94 if args.confirm_data and (args.confirm_accuracy_limit is None or args.confirm_rsd_limit is None):
95 raise InputError("QL confirmation needs pre-stated --confirm-accuracy-limit and --confirm-rsd-limit")
97 rows = parse_rows(read_input(args.calibration), args.calibration)
98 require_columns(rows, ["level", "response"])
99 xs = [to_float(r["level"], "level", i) for i, r in enumerate(rows)]
100 ys = [to_float(r["response"], "response", i) for i, r in enumerate(rows)]
101 fit = fit_linear(xs, ys)
102 slope = fit.slope
103 if slope == 0:
104 raise InputError("calibration slope is zero; DL/QL cannot be computed")
106 findings: list[str] = []
107 results: list[dict] = []
109 # Approach 3.2.3.3, sigma = residual SD of the regression line.
110 results.append(
111 {
112 "approach": "sd-and-slope (sigma = residual SD of regression)",
113 "sigma": fit.residual_sd,
114 "slope": slope,
115 "DL": DL_FACTOR * fit.residual_sd / abs(slope),
116 "QL": QL_FACTOR * fit.residual_sd / abs(slope),
117 "reference": "Q2(R2) 3.2.3.3",
118 }
119 )
121 # Q2 refers to SD across regression-line intercepts, not the standard error
122 # of the intercept from one line. Do not silently substitute the latter.
123 if args.intercepts:
124 irows = parse_rows(read_input(args.intercepts), args.intercepts)
125 require_columns(irows, ["intercept"])
126 intercepts = [to_float(r["intercept"], "intercept", i) for i, r in enumerate(irows)]
127 if len(intercepts) < 3:
128 raise InputError("intercept SD needs at least 3 independent calibration curves (tool minimum)")
129 sigma = sample_sd(intercepts)
130 results.append({"approach": "sd-and-slope (sigma = SD of independent y-intercepts)",
131 "sigma": sigma, "slope": slope,
132 "DL": DL_FACTOR * sigma / abs(slope),
133 "QL": QL_FACTOR * sigma / abs(slope),
134 "reference": "Q2(R2) 3.2.3.3"})
136 # Approach 3.2.3.3, sigma = SD of blank responses.
137 if args.blanks:
138 brows = parse_rows(read_input(args.blanks), args.blanks)
139 require_columns(brows, ["response"])
140 blanks = [to_float(r["response"], "response", i) for i, r in enumerate(brows)]
141 if len(blanks) < 3:
142 raise InputError("blank SD needs at least 3 blank measurements")
143 bsd = sample_sd(blanks)
144 results.append(
145 {
146 "approach": f"sd-and-slope (sigma = SD of {len(blanks)} blanks)",
147 "sigma": bsd,
148 "slope": slope,
149 "DL": DL_FACTOR * bsd / abs(slope),
150 "QL": QL_FACTOR * bsd / abs(slope),
151 "reference": "Q2(R2) 3.2.3.3",
152 }
153 )
154 note(f"blank mean response {mean(blanks):.6g}, SD {bsd:.6g}, n={len(blanks)}")
156 # Approach 3.2.3.2, signal-to-noise.
157 if args.signal_to_noise is not None:
158 if args.sn_level is None:
159 raise InputError("--signal-to-noise needs --sn-level")
160 if args.signal_to_noise <= 0 or args.sn_level <= 0:
161 raise InputError("--signal-to-noise and --sn-level must be > 0")
162 per_unit = args.signal_to_noise / args.sn_level
163 results.append(
164 {
165 "approach": f"signal-to-noise (S/N {args.signal_to_noise:g} at {args.sn_level:g})",
166 "sigma": float("nan"),
167 "slope": slope,
168 "DL": 3.0 / per_unit,
169 "QL": 10.0 / per_unit,
170 "reference": "Q2(R2) 3.2.3.2",
171 }
172 )
173 note(
174 "signal-to-noise DL uses the 3:1 ratio and QL the 10:1 ratio from Q2(R2) 3.2.3.2, "
175 "scaled linearly from the measured S/N. Linear scaling of noise is an assumption -- "
176 "confirm at the resulting level"
177 )
179 if any(result["QL"] <= 0 for result in results):
180 findings.append("zero estimated response variation cannot establish a physical detection/quantitation limit")
182 # Spread across approaches: the point of computing all of them.
183 qls = [r["QL"] for r in results if math.isfinite(r["QL"])]
184 spread_note = ""
185 if len(qls) >= 2:
186 ratio = max(qls) / min(qls) if min(qls) > 0 else float("inf")
187 spread_note = (
188 f"QL estimates span {min(qls):.6g} to {max(qls):.6g} ({ratio:.2f}x) across "
189 f"{len(qls)} approaches"
190 )
191 if ratio > 2.0:
192 findings.append(
193 spread_note
194 + ": name the approach used in the report, and confirm the claimed limit with "
195 "real data at that level (Q2(R2) 3.2.3.5)"
196 )
198 # QL confirmation with real data (Q2(R2) 3.2.3.4 / 3.2.3.5).
199 confirm_rows: list[dict] = []
200 if args.confirm_data:
201 if args.confirm_ql is None:
202 raise InputError("--confirm-data needs --confirm-ql")
203 crows = parse_rows(read_input(args.confirm_data), args.confirm_data)
204 require_columns(crows, ["measured"])
205 measured = [to_float(r["measured"], "measured", i) for i, r in enumerate(crows)]
206 if len(measured) < 3:
207 raise InputError("QL confirmation needs at least 3 determinations")
208 m = mean(measured)
209 bias = 100.0 * (m - args.confirm_ql) / args.confirm_ql
210 rsd = rsd_percent(measured)
211 confirm_rows = [
212 {"metric": "claimed QL", "value": args.confirm_ql},
213 {"metric": "n determinations", "value": len(measured)},
214 {"metric": "mean measured", "value": m},
215 {"metric": "bias vs claimed QL (%)", "value": bias},
216 {"metric": "RSD (%)", "value": rsd},
217 {"metric": "accuracy limit (%)", "value": args.confirm_accuracy_limit},
218 {"metric": "RSD limit (%)", "value": args.confirm_rsd_limit},
219 ]
220 if abs(bias) > args.confirm_accuracy_limit:
221 findings.append(
222 f"QL confirmation: bias {bias:+.2f}% at the claimed QL exceeds "
223 f"+/-{args.confirm_accuracy_limit:g}%"
224 )
225 if not math.isfinite(rsd):
226 findings.append("QL confirmation: RSD is undefined; precision was not demonstrated")
227 elif rsd > args.confirm_rsd_limit:
228 findings.append(
229 f"QL confirmation: {rsd:.2f}% RSD at the claimed QL exceeds "
230 f"{args.confirm_rsd_limit:g}%"
231 )
232 elif args.confirm_ql is not None:
233 findings.append(
234 "a QL was claimed but no confirmation data supplied. Q2(R2) 3.2.3.5 asks that an "
235 "estimated limit be validated by analysing samples at or near it"
236 )
238 # The requirement that actually gates an impurity method.
239 if args.reporting_threshold is not None:
240 # With no claimed QL, use the LARGEST estimate. Taking the smallest would
241 # let the check pass on the most flattering choice of sigma, which is the
242 # wrong direction to err on a compliance requirement.
243 conservative = max(qls) if qls else float("nan")
244 claimed = args.confirm_ql if args.confirm_ql is not None else conservative
245 if args.confirm_ql is None and qls:
246 note(
247 f"no --confirm-ql given, so the reporting-threshold check uses the most "
248 f"conservative estimate ({conservative:.6g}), not the most favourable "
249 f"({min(qls):.6g})"
250 )
251 if min(qls) <= args.reporting_threshold < conservative:
252 findings.append(
253 f"the QL estimates straddle the reporting threshold "
254 f"{args.reporting_threshold:.6g}: {min(qls):.6g} would pass and "
255 f"{conservative:.6g} would not. Whether this procedure meets Q2(R2) "
256 "3.2.3.5 depends on which approach is chosen, so choose it, justify it, "
257 "and confirm the claimed limit with data at that level"
258 )
259 if math.isfinite(claimed) and claimed > args.reporting_threshold:
260 findings.append(
261 f"QL {claimed:.6g} is above the reporting threshold "
262 f"{args.reporting_threshold:.6g}. Q2(R2) 3.2.3.5 requires the QL for an "
263 "impurity procedure to be at or below the reporting threshold"
264 )
265 elif math.isfinite(claimed):
266 ratio = args.reporting_threshold / claimed if claimed > 0 else float("inf")
267 if math.isclose(ratio, 1.0, rel_tol=1e-9):
268 note(
269 f"QL {claimed:.6g} sits exactly at the reporting threshold "
270 f"{args.reporting_threshold:.6g}, which satisfies Q2(R2) 3.2.3.5 with no "
271 "margin -- any drift puts the procedure out of compliance"
272 )
273 else:
274 note(
275 f"QL {claimed:.6g} is {ratio:.1f}x below the reporting threshold "
276 f"{args.reporting_threshold:.6g}"
277 )
278 if ratio >= 10:
279 note(
280 "the QL is roughly 10x or more below the reporting limit, so Q2(R2) 3.2.3.5 "
281 "allows the confirmatory validation to be omitted with justification"
282 )
284 if args.format == "json":
285 emit(
286 [
287 {
288 "calibration": {
289 "n": fit.n,
290 "slope": slope,
291 "intercept": fit.intercept,
292 "residual_sd": fit.residual_sd,
293 "weighting": args.weight,
294 },
295 "estimates": results,
296 "confirmation": confirm_rows,
297 "findings": findings,
298 }
299 ],
300 "json",
301 )
302 else:
303 emit(results, args.format)
304 if confirm_rows:
305 print()
306 emit(confirm_rows, args.format)
308 if not args.confirm_data:
309 note(
310 f"approach 'accuracy-precision' ({DL_QL_APPROACHES['accuracy-precision']['note']}) "
311 "validates the QL directly rather than estimating it -- supply --confirm-data "
312 "to use it"
313 )
314 if spread_note and not any(spread_note in f for f in findings):
315 note(spread_note)
316 note("Q2(R2) 3.2.3.5: report the limit AND the approach used to determine it")
317 for f in findings:
318 finding(f)
319 if not findings:
320 note("no findings against the checks that were run")
321 note("this tool does not decide that a detection or quantitation limit is acceptable")
322 return EXIT_FINDINGS if findings else EXIT_OK
325if __name__ == "__main__":
326 run_cli(main)