Skills › Scientific Agent Skills › analytical-method-validation › check_response.py
check_response.py
251 lines · skills/analytical-method-validation/scripts/check_response.py at 5f2bc03
1#!/usr/bin/env python3
2"""Evaluate a calibration response (linearity) the way ICH Q2(R2) 3.2.2 asks.
4Reports what the guideline asks to be reported -- slope, intercept, coefficient
5of determination, and an analysis of the deviation of points from the line --
6and adds the diagnostics that actually detect an unsuitable model: a lack-of-fit
7F test against pure error, a runs test on residual signs, back-calculated
8relative error per level, and a heteroscedasticity check that tells you whether
9weighting is needed.
11 python3 check_response.py --input calib.csv
12 python3 check_response.py --input calib.csv --weight 1/x2 --levels-required 5
13 python3 check_response.py --input calib.csv --max-back-calc-error 5 --format json
15Input columns: `level` (nominal concentration or %) and `response` (signal).
16An optional `replicate` column is ignored -- replicates are simply repeated rows
17at the same level, which is what enables the lack-of-fit test.
19Exit codes: 0 no findings, 1 findings raised, 2 bad input.
20"""
22from __future__ import annotations
24import argparse
25import math
26import sys
28sys.path.insert(0, str(__import__("pathlib").Path(__file__).resolve().parent))
30from _common import ( # noqa: E402
31 EXIT_FINDINGS,
32 EXIT_OK,
33 InputError,
34 add_common_args,
35 emit,
36 finding,
37 finite_float,
38 fit_linear,
39 heteroscedasticity,
40 lack_of_fit,
41 mean,
42 note,
43 parse_rows,
44 read_input,
45 require_columns,
46 run_cli,
47 runs_test,
48 to_float,
49)
51WEIGHT_SCHEMES = {
52 "none": lambda x: 1.0,
53 "1/x": lambda x: 1.0 / abs(x) if x != 0 else 0.0,
54 "1/x2": lambda x: 1.0 / (x * x) if x != 0 else 0.0,
55}
58def build_weights(xs: list[float], scheme: str) -> list[float] | None:
59 if scheme == "none":
60 return None
61 fn = WEIGHT_SCHEMES[scheme]
62 if any(x == 0 for x in xs):
63 raise InputError(f"weighting {scheme} needs non-zero levels; a zero level was supplied")
64 return [fn(x) for x in xs]
67def main() -> int:
68 parser = argparse.ArgumentParser(description="Check a calibration response for linearity.")
69 parser.add_argument("--input", "-i", help="CSV/TSV/JSON file, or '-' for stdin")
70 parser.add_argument("--weight", choices=sorted(WEIGHT_SCHEMES), default="none",
71 help="calibration weighting (default: none)")
72 parser.add_argument("--levels-required", type=int, default=5,
73 help="minimum distinct levels expected (Q2(R2) recommends 5)")
74 parser.add_argument("--max-back-calc-error", type=finite_float, default=None,
75 help="flag any level whose mean back-calculated error exceeds this %%")
76 parser.add_argument("--alpha", type=finite_float, default=0.05,
77 help="significance level for lack-of-fit and runs tests")
78 parser.add_argument("--through-origin-tolerance", type=finite_float, default=None,
79 help="flag when the intercept exceeds this %% of the response at the "
80 "highest level (a proxy for the y-intercept significance check)")
81 add_common_args(parser)
82 args = parser.parse_args()
84 if not 0 < args.alpha < 1:
85 raise InputError("--alpha must be between 0 and 1")
86 if args.levels_required < 3:
87 raise InputError("--levels-required must be at least 3")
88 for value in (args.max_back_calc_error, args.through_origin_tolerance):
89 if value is not None and value < 0:
90 raise InputError("error/intercept limits must be non-negative")
92 rows = parse_rows(read_input(args.input), args.input)
93 require_columns(rows, ["level", "response"])
94 xs = [to_float(r["level"], "level", i) for i, r in enumerate(rows)]
95 ys = [to_float(r["response"], "response", i) for i, r in enumerate(rows)]
97 weights = build_weights(xs, args.weight)
98 fit = fit_linear(xs, ys, weights)
99 distinct = sorted({round(x, 12) for x in xs})
101 findings: list[str] = []
103 # Q2(R2) 3.2.2.1: a minimum of five concentrations is recommended.
104 if len(distinct) < args.levels_required:
105 findings.append(
106 f"{len(distinct)} distinct levels; Q2(R2) 3.2.2.1 recommends at least "
107 f"{args.levels_required} appropriately distributed across the range"
108 )
110 lof = lack_of_fit(xs, ys, fit)
111 runs = runs_test([r for _, r in sorted(zip(xs, fit.residuals), key=lambda pair: pair[0])])
112 het = heteroscedasticity(xs, fit.residuals)
114 if lof.get("applicable") and lof["p_value"] < args.alpha:
115 findings.append(
116 f"lack-of-fit F={lof['f_statistic']:.3f} on {lof['df_lack_of_fit']}/"
117 f"{lof['df_pure_error']} df, p={lof['p_value']:.4g}: the straight line does not "
118 "describe the data beyond replicate scatter"
119 )
120 if math.isfinite(runs.get("p_value", float("nan"))) and runs["p_value"] < args.alpha:
121 findings.append(
122 f"residual signs are non-random (runs={runs['runs']}, expected "
123 f"{runs['expected_runs']:.1f}, p={runs['p_value']:.4g}): inspect the residual plot "
124 "for curvature"
125 )
126 if het.get("applicable") and het["variance_ratio_high_over_low"] > 10 and args.weight == "none":
127 findings.append(
128 f"residual variance is {het['variance_ratio_high_over_low']:.1f}x larger in the top "
129 "third of the range than the bottom, and the fit is unweighted: back-calculated "
130 "results at the low end may be imprecise. Investigate a justified variance model"
131 )
133 # Back-calculated relative error per level -- the practical test of the model.
134 level_rows = []
135 by_level: dict[float, list[float]] = {}
136 for x, y in zip(xs, ys):
137 by_level.setdefault(round(x, 12), []).append(y)
138 for level in distinct:
139 responses = by_level[level]
140 back = [
141 (r - fit.intercept) / fit.slope if fit.slope != 0 else float("nan")
142 for r in responses
143 ]
144 mean_back = mean(back)
145 rel_err = 100.0 * (mean_back - level) / level if level != 0 else float("nan")
146 level_rows.append(
147 {
148 "level": level,
149 "n": len(responses),
150 "mean_response": mean(responses),
151 "mean_back_calculated": mean_back,
152 "relative_error_pct": rel_err,
153 }
154 )
155 if (
156 args.max_back_calc_error is not None
157 and math.isfinite(rel_err)
158 and abs(rel_err) > args.max_back_calc_error
159 ):
160 findings.append(
161 f"level {level:g}: back-calculated mean deviates {rel_err:+.2f}% from nominal, "
162 f"outside the stated +/-{args.max_back_calc_error:g}%"
163 )
165 if args.through_origin_tolerance is not None:
166 top_response = fit.predict(max(distinct))
167 if top_response != 0:
168 pct = 100.0 * abs(fit.intercept) / abs(top_response)
169 if pct > args.through_origin_tolerance:
170 findings.append(
171 f"intercept is {pct:.2f}% of the response at the highest level, above the "
172 f"stated {args.through_origin_tolerance:g}% tolerance"
173 )
175 slope_lo, slope_hi = fit.slope_ci()
176 int_lo, int_hi = fit.intercept_ci()
177 summary = [
178 {"statistic": "n points", "value": fit.n},
179 {"statistic": "distinct levels", "value": len(distinct)},
180 {"statistic": "weighting", "value": args.weight},
181 {"statistic": "slope", "value": fit.slope},
182 {"statistic": "slope 95% CI", "value": f"{slope_lo:.6g} to {slope_hi:.6g}"},
183 {"statistic": "intercept", "value": fit.intercept},
184 {"statistic": "intercept 95% CI", "value": f"{int_lo:.6g} to {int_hi:.6g}"},
185 {"statistic": "intercept CI includes 0", "value": int_lo <= 0.0 <= int_hi},
186 {"statistic": "coefficient of determination (r2)", "value": fit.r_squared},
187 {"statistic": "correlation coefficient (r)", "value": fit.r},
188 {"statistic": "residual SD", "value": fit.residual_sd},
189 {"statistic": "residual df", "value": fit.df},
190 ]
191 if lof.get("applicable"):
192 summary += [
193 {"statistic": "lack-of-fit F", "value": lof["f_statistic"]},
194 {"statistic": "lack-of-fit p", "value": lof["p_value"]},
195 ]
196 else:
197 summary.append(
198 {"statistic": "lack-of-fit test", "value": f"not run ({lof.get('reason', '')})"}
199 )
200 summary += [
201 {"statistic": "runs test p", "value": runs.get("p_value", float("nan"))},
202 {
203 "statistic": "residual SD ratio (high/low third)",
204 "value": het.get("sd_ratio", float("nan")),
205 },
206 ]
208 if args.format == "json":
209 emit(
210 [
211 {
212 "summary": {r["statistic"]: r["value"] for r in summary},
213 "levels": level_rows,
214 "findings": findings,
215 }
216 ],
217 "json",
218 )
219 else:
220 emit(summary, args.format)
221 print()
222 emit(level_rows, args.format)
224 note(
225 "Q2(R2) 3.2.2.1 asks for the plot, r or r-squared, intercept, slope, and an analysis of "
226 "the deviation of points from the line"
227 )
228 note(
229 "r-squared alone does not demonstrate linearity: it rises with range and is insensitive "
230 "to curvature. The lack-of-fit test and the residual pattern are the evidence"
231 )
232 if args.weight != "none" and lof.get("applicable"):
233 note(
234 f"weighted lack-of-fit uses {args.weight} in both residual and pure-error sums; "
235 "its F reference assumes independent normal errors with variance proportional "
236 "to inverse weight. The weights need independent scientific justification"
237 )
238 if not lof.get("applicable"):
239 note(
240 f"lack-of-fit was not assessed: {lof.get('reason', 'unavailable')}"
241 )
242 for f in findings:
243 finding(f)
244 if not findings:
245 note("no findings against the checks that were run")
246 note("this tool does not decide that the calibration model is acceptable")
247 return EXIT_FINDINGS if findings else EXIT_OK
250if __name__ == "__main__":
251 run_cli(main)