Skills › Scientific Agent Skills › analytical-method-validation › compare_methods.py
compare_methods.py
218 lines · skills/analytical-method-validation/scripts/compare_methods.py at 5f2bc03
1#!/usr/bin/env python3
2"""Compare two analytical procedures for a transfer, bridging, or bias study.
4Uses the regressions that belong to method comparison -- Deming and
5Passing-Bablok, which allow error in both measurements -- rather than ordinary
6least squares, which assumes the comparative procedure is error-free and biases
7the slope toward zero. Adds Bland-Altman agreement and, most importantly, a TOST
8equivalence test.
10The default reflex at a method transfer is a t test, and "p > 0.05, no
11significant difference" is then written up as evidence of equivalence. It is not:
12failing to detect a difference is not the same as showing there is none, and with
13a small transfer dataset that outcome is close to guaranteed. TOST tests the
14hypothesis that matters -- that the true difference lies inside a pre-stated
15acceptance margin.
17 python3 compare_methods.py --input paired.csv --margin 2
18 python3 compare_methods.py -i paired.csv --margin 2 --relative --lambda 1.0
20Input columns: `reference` and `test` (paired results on the same samples).
21An optional `sample` column labels the rows.
23Exit codes: 0 equivalence demonstrated at the stated margin, 1 not demonstrated
24or other findings, 2 bad input.
25"""
27from __future__ import annotations
29import argparse
30import math
31import sys
33sys.path.insert(0, str(__import__("pathlib").Path(__file__).resolve().parent))
35from _common import ( # noqa: E402
36 EXIT_FINDINGS,
37 EXIT_OK,
38 InputError,
39 add_common_args,
40 bland_altman,
41 deming,
42 emit,
43 finding,
44 finite_float,
45 fit_linear,
46 mean,
47 note,
48 parse_rows,
49 passing_bablok,
50 read_input,
51 require_columns,
52 run_cli,
53 t_cdf,
54 to_float,
55 tost_paired,
56)
59def main() -> int:
60 parser = argparse.ArgumentParser(
61 description="Compare two procedures with the statistics method comparison requires."
62 )
63 parser.add_argument("--input", "-i", help="CSV/TSV/JSON with `reference` and `test`")
64 parser.add_argument("--margin", type=finite_float, required=True,
65 help="pre-stated equivalence margin for TOST, in the units of the "
66 "difference (absolute, or %% when --relative is used)")
67 parser.add_argument("--relative", action="store_true",
68 help="work in percent differences relative to the pair mean")
69 parser.add_argument("--lambda", dest="lambda_ratio", type=finite_float, default=1.0,
70 help="Deming error-variance ratio var(y error)/var(x error) -- test "
71 "procedure over comparative procedure. Default 1.0 means equal "
72 "precision. Estimate it as (SD of y replicates / SD of x "
73 "replicates) squared")
74 parser.add_argument("--alpha", type=finite_float, default=0.05,
75 help="one-sided alpha for TOST (default 0.05)")
76 parser.add_argument("--slope-tolerance", type=finite_float, default=None,
77 help="flag when the slope CI excludes 1 +/- this amount")
78 add_common_args(parser)
79 args = parser.parse_args()
81 if args.slope_tolerance is not None and args.slope_tolerance < 0:
82 raise InputError("--slope-tolerance must be non-negative")
83 rows = parse_rows(read_input(args.input), args.input)
84 require_columns(rows, ["reference", "test"])
85 ref = [to_float(r["reference"], "reference", i) for i, r in enumerate(rows)]
86 test = [to_float(r["test"], "test", i) for i, r in enumerate(rows)]
87 n = len(ref)
88 if n < 3:
89 raise InputError("method comparison needs at least 3 paired results")
91 findings: list[str] = []
93 # Agreement.
94 ba = bland_altman(ref, test, relative=args.relative)
95 unit = "%" if args.relative else "units"
97 # TOST on the differences that Bland-Altman used.
98 if args.relative:
99 diffs = [
100 100.0 * (t - r) / (0.5 * (r + t)) for r, t in zip(ref, test)
101 ]
102 else:
103 diffs = [t - r for r, t in zip(ref, test)]
104 tost = tost_paired(diffs, args.margin, args.alpha)
106 # The naive test, computed only to show what it does not establish.
107 d_mean = mean(diffs)
108 sd = ba["sd_differences"]
109 se = sd / math.sqrt(n) if sd > 0 else float("nan")
110 t_stat = d_mean / se if se and math.isfinite(se) and se > 0 else float("nan")
111 p_naive = (
112 2.0 * (1.0 - t_cdf(abs(t_stat), n - 1)) if math.isfinite(t_stat) else float("nan")
113 )
115 # Regressions.
116 ols = fit_linear(ref, test)
117 dem = deming(ref, test, args.lambda_ratio)
118 try:
119 pb = passing_bablok(ref, test)
120 except InputError as exc:
121 pb = None
122 note(f"Passing-Bablok not computed: {exc}")
124 summary = [
125 {"statistic": "n pairs", "value": n},
126 {"statistic": f"mean difference ({unit})", "value": ba["bias"]},
127 {"statistic": "difference 95% CI", "value":
128 f"{ba['bias_ci95'][0]:.6g} to {ba['bias_ci95'][1]:.6g}"},
129 {"statistic": f"SD of differences ({unit})", "value": sd},
130 {"statistic": "limits of agreement", "value":
131 f"{ba['loa_lower']:.6g} to {ba['loa_upper']:.6g}"},
132 {"statistic": "LoA 95% CI half-width", "value": ba["loa_ci95_halfwidth"]},
133 {"statistic": "--- equivalence ---", "value": ""},
134 {"statistic": "TOST margin", "value": args.margin},
135 {"statistic": "TOST p-value", "value": tost["p_value"]},
136 {"statistic": f"{100 * (1 - 2 * args.alpha):.0f}% CI (TOST)", "value":
137 f"{tost['ci_1_minus_2alpha'][0]:.6g} to {tost['ci_1_minus_2alpha'][1]:.6g}"},
138 {"statistic": "equivalent at stated margin", "value": tost["equivalent"]},
139 {"statistic": "--- for contrast only ---", "value": ""},
140 {"statistic": "paired t-test p (NOT equivalence)", "value": p_naive},
141 {"statistic": "--- regressions ---", "value": ""},
142 {"statistic": "OLS slope (biased here)", "value": ols.slope},
143 {"statistic": "Deming slope", "value": dem["slope"]},
144 {"statistic": "Deming slope 95% CI", "value":
145 f"{dem['slope_ci95'][0]:.6g} to {dem['slope_ci95'][1]:.6g}"},
146 {"statistic": "Deming intercept", "value": dem["intercept"]},
147 ]
148 if pb:
149 summary += [
150 {"statistic": "Passing-Bablok slope", "value": pb["slope"]},
151 {"statistic": "Passing-Bablok slope 95% CI", "value":
152 f"{pb['slope_ci95'][0]:.6g} to {pb['slope_ci95'][1]:.6g}"},
153 {"statistic": "Passing-Bablok intercept", "value": pb["intercept"]},
154 ]
156 prop = ba.get("proportional_bias")
157 if prop and math.isfinite(prop.get("p_value", float("nan"))):
158 summary.append(
159 {"statistic": "proportional bias p (difference vs mean)", "value": prop["p_value"]}
160 )
161 if prop["p_value"] < 0.05:
162 findings.append(
163 f"the difference trends with concentration (slope {prop['slope']:.4g}, "
164 f"p={prop['p_value']:.4g}): a single mean bias does not describe the "
165 "disagreement, and limits of agreement are misleading"
166 )
168 if not tost["equivalent"]:
169 findings.append(
170 f"equivalence NOT demonstrated at +/-{args.margin:g} {unit}: the "
171 f"{100 * (1 - 2 * args.alpha):.0f}% CI "
172 f"({tost['ci_1_minus_2alpha'][0]:.4g}, {tost['ci_1_minus_2alpha'][1]:.4g}) is not "
173 f"contained in the margin"
174 )
176 if args.slope_tolerance is not None:
177 lo, hi = dem["slope_ci95"]
178 target_lo, target_hi = 1.0 - args.slope_tolerance, 1.0 + args.slope_tolerance
179 if not (math.isfinite(lo) and math.isfinite(hi)) or lo < target_lo or hi > target_hi:
180 findings.append(
181 f"Deming slope 95% CI ({lo:.4g}, {hi:.4g}) is not contained in "
182 f"({target_lo:g}, {target_hi:g})"
183 )
185 if args.format == "json":
186 emit([{"summary": {r["statistic"]: r["value"] for r in summary
187 if not r["statistic"].startswith("---")},
188 "bland_altman": ba, "deming": dem,
189 "passing_bablok": pb, "tost": tost,
190 "ols_slope": ols.slope, "paired_t_p_value": p_naive,
191 "findings": findings}], "json")
192 else:
193 emit(summary, args.format)
195 if math.isfinite(p_naive) and p_naive > 0.05 and not tost["equivalent"]:
196 note(
197 f"the paired t-test gives p={p_naive:.4g}, which would often be written up as "
198 "'no significant difference'. TOST shows equivalence is NOT established at the "
199 "stated margin. Absence of a detected difference is not evidence of equivalence"
200 )
201 note(
202 "OLS treats reference values as fixed without error; appreciable reference error can "
203 "attenuate the slope. Choose a regression consistent with the error model"
204 )
205 note(
206 "the equivalence margin must be pre-stated from the specification or the analytical "
207 "target profile, never chosen after seeing the data"
208 )
209 for f in findings:
210 finding(f)
211 if not findings:
212 note("no findings against the checks that were run")
213 note("this tool does not decide that a transfer or comparison passes")
214 return EXIT_FINDINGS if findings else EXIT_OK
217if __name__ == "__main__":
218 run_cli(main)