Skills › Scientific Agent Skills › pkpd-modeling › bioequivalence.py

bioequivalence.py

420 lines · skills/pkpd-modeling/scripts/bioequivalence.py at 5f2bc03

Open on GitHub

1#!/usr/bin/env python3

2"""Conventional average bioequivalence and numerical sample-size planning.

4The CLI supports complete RT/TR 2x2 crossover and independent parallel data.

5It rejects replicate/reference-scaled analysis, which requires a design-specific

6model. Retained scalar scaling functions are arithmetic aids only.

8 python3 bioequivalence.py -i be.csv --design 2x2 --metric auc

9 python3 bioequivalence.py --power --cv 0.30 --gmr 0.95 --target-power 0.80

11Input: subject,treatment(T/R),value, plus sequence and period for crossover.

12"""

14from __future__ import annotations

16import argparse

17import math

18from dataclasses import dataclass

19from typing import Sequence

21import numpy as np

23from _common import (

24 InputError,

25 Report,

26 add_format_argument,

27 main_wrapper,

28 parse_float,

29 read_table,

30 require_columns,

31)

33try:

34 from scipy.stats import chi2, norm

35 from scipy.stats import t as t_dist

36except ImportError as exc: # pragma: no cover

37 raise SystemExit("bioequivalence.py needs scipy: uv pip install scipy") from exc

40# The FDA regulatory constant: ln(1.25) / 0.25, the point at which the scaled

41# criterion coincides with the conventional limits at sigma_w0 = 0.25.

42THETA_FDA = math.log(1.25) / 0.25

43# EMA's ABEL widening coefficient, and the CVwR at which widening is capped.

44K_ABEL = 0.760

45CV_ABEL_CAP = 0.50

46SWR_ABEL_CAP = math.sqrt(math.log(1.0 + CV_ABEL_CAP**2))

47CV_SCALING_THRESHOLD = 0.30

50def cv_from_s2(s2: float) -> float:

51 return math.sqrt(math.exp(s2) - 1.0)

54def s2_from_cv(cv: float) -> float:

55 return math.log(1.0 + cv**2)

58# ------------------------------------------------------------------ average BE

61@dataclass

62class AverageBE:

63 gmr: float

64 ci_low: float

65 ci_high: float

66 df: int

67 se: float

68 cv_within: float | None

69 n_subjects: int

70 method: str

73def crossover_2x2(records: list[dict], design: str = "2x2") -> AverageBE:

74 """Crossover analysis via within-subject period differences.

76 Averaging the sequence-specific mean differences cancels the period effect

77 exactly, which is the whole reason a crossover is run. The estimate and its

78 standard error are identical to what a sequence/subject(sequence)/period/

79 treatment ANOVA gives on balanced or unbalanced data, without building the

80 design matrix.

82 Replicate designs are rejected; repeated treatments need a different model.

83 """

84 if design != "2x2":

85 raise InputError("replicate designs need a design-specific period/sequence model; use FDA Appendix G or validated replicateBE tooling")

86 by_subject: dict[str, dict[str, list[float]]] = {}

87 sequences: dict[str, str] = {}

88 for row in records:

89 by_subject.setdefault(row["subject"], {}).setdefault(row["treatment"], []).append(row["logvalue"])

90 if row.get("sequence"):

91 sequences[row["subject"]] = row["sequence"]

93 for subject, treatments in by_subject.items():

94 if set(treatments) != {"T", "R"} or any(len(v) != 1 for v in treatments.values()):

95 raise InputError(f"subject {subject}: 2x2 analysis requires exactly one T and one R; handle incomplete records explicitly")

96 if sequences.get(subject) not in {"RT", "TR"}:

97 raise InputError("2x2 analysis requires RT/TR sequence and periods 1/2")

98 subject_rows = [r for r in records if r["subject"] == subject]

99 if {str(r.get("period")) for r in subject_rows} != {"1", "2"}:

100 raise InputError(f"subject {subject}: expected periods 1 and 2")

101 if any(r.get("sequence") != sequences[subject] or r["treatment"] != sequences[subject][int(r["period"])-1] for r in subject_rows):

102 raise InputError(f"subject {subject}: treatment, period and sequence disagree")

103 if set(sequences.values()) != {"RT", "TR"}:

104 raise InputError("both RT and TR sequences are required to separate period and treatment")

105 differences: dict[str, list[float]] = {}

106 for subject, treatments in by_subject.items():

107 if "T" not in treatments or "R" not in treatments:

108 continue # incomplete subject: dropped, as in the standard analysis

109 diff = float(np.mean(treatments["T"])) - float(np.mean(treatments["R"]))

110 sequence = sequences.get(subject, "?")

111 differences.setdefault(sequence, []).append(diff)

112

113 total = sum(len(v) for v in differences.values())

114 if total < 3:

115 raise InputError("fewer than 3 subjects completed both treatments")

116

117 means = {seq: float(np.mean(v)) for seq, v in differences.items()}

118 estimate = float(np.mean(list(means.values())))

119 pooled_num = sum(float(np.var(v, ddof=1)) * (len(v) - 1) for v in differences.values() if len(v) > 1)

120 pooled_den = sum(len(v) - 1 for v in differences.values())

121 s2d = pooled_num / pooled_den

122 se = math.sqrt(s2d / 4.0 * sum(1.0 / len(v) for v in differences.values()))

123 df = pooled_den

124 s2w = s2d / 2.0

125 method = "2x2 crossover (period effect removed)"

126

127 crit = float(t_dist.ppf(0.95, df))

128 return AverageBE(

129 gmr=math.exp(estimate),

130 ci_low=math.exp(estimate - crit * se),

131 ci_high=math.exp(estimate + crit * se),

132 df=df,

133 se=se,

134 cv_within=cv_from_s2(s2w) if s2w > 0 else None,

135 n_subjects=total,

136 method=method,

137 )

138

139

140def parallel_design(records: list[dict]) -> AverageBE:

141 if len({r["subject"] for r in records}) != len(records):

142 raise InputError("parallel analysis requires one independent record per subject")

143 test = np.asarray([r["logvalue"] for r in records if r["treatment"] == "T"])

144 ref = np.asarray([r["logvalue"] for r in records if r["treatment"] == "R"])

145 if len(test) < 2 or len(ref) < 2:

146 raise InputError("a parallel design needs at least 2 subjects per arm")

147 n1, n2 = len(test), len(ref)

148 pooled = ((n1 - 1) * test.var(ddof=1) + (n2 - 1) * ref.var(ddof=1)) / (n1 + n2 - 2)

149 se = math.sqrt(pooled * (1.0 / n1 + 1.0 / n2))

150 df = n1 + n2 - 2

151 estimate = float(test.mean() - ref.mean())

152 crit = float(t_dist.ppf(0.95, df))

153 return AverageBE(

154 gmr=math.exp(estimate),

155 ci_low=math.exp(estimate - crit * se),

156 ci_high=math.exp(estimate + crit * se),

157 df=df,

158 se=se,

159 cv_within=cv_from_s2(float(pooled)), # total, not within-subject

160 n_subjects=n1 + n2,

161 method="parallel (pooled variance; the CV shown is total, not within-subject)",

162 )

163

164

165# -------------------------------------------------------------- scaled BE

166

167

168@dataclass

169class ReferenceVariability:

170 s2wr: float

171 df: int

172 n_subjects: int

173

174 @property

175 def swr(self) -> float:

176 return math.sqrt(self.s2wr)

177

178 @property

179 def cvwr(self) -> float:

180 return cv_from_s2(self.s2wr)

181

182

183def reference_variability(records: list[dict]) -> ReferenceVariability:

184 """Descriptive unadjusted replicate variance; NOT design-adjusted BE variance.

185

186 Period/sequence effects contaminate this quantity. Do not supply it to a

187 regulatory scaled analysis; the CLI rejects that unsupported workflow.

188 """

189 numerator = 0.0

190 df = 0

191 subjects = 0

192 for subject in {r["subject"] for r in records}:

193 values = [r["logvalue"] for r in records if r["subject"] == subject and r["treatment"] == "R"]

194 if len(values) < 2:

195 continue

196 numerator += float(np.var(values, ddof=1)) * (len(values) - 1)

197 df += len(values) - 1

198 subjects += 1

199 if df == 0:

200 raise InputError(

201 "no subject received the reference more than once. Reference-scaling requires a replicate "

202 "design (partial replicate RRT/RTR/TRR or full replicate RTRT/TRTR) - it cannot be applied "

203 "to a 2x2 study whatever the observed variability."

204 )

205 return ReferenceVariability(s2wr=numerator / df, df=df, n_subjects=subjects)

206

207

208def abel_limits(rv: ReferenceVariability) -> tuple[float, float, bool]:

209 """EMA widened acceptance limits. Returns (low, high, widened)."""

210 if rv.cvwr <= CV_SCALING_THRESHOLD:

211 return 0.80, 1.25, False

212 swr = min(rv.swr, SWR_ABEL_CAP)

213 return math.exp(-K_ABEL * swr), math.exp(K_ABEL * swr), True

214

215

216def rsabe_bound(estimate: float, se: float, df_point: int, rv: ReferenceVariability) -> dict[str, float]:

217 """Scalar FDA 2026 Appendix G bound; inputs must come from its design model.

218

219 This does not establish the HVD applicability threshold, point-estimate

220 constraint or validity of externally supplied variance/degrees of freedom.

221 """

222 if not all(math.isfinite(v) for v in (estimate, se, df_point, rv.s2wr, rv.df)) or se < 0 or df_point <= 0 or rv.s2wr < 0 or rv.df <= 0:

223 raise InputError("scaled bound requires finite contrast, non-negative SE/variance and positive degrees of freedom")

224 e_point = estimate**2 - se**2

225 e_bound = (abs(estimate) + float(t_dist.ppf(0.95, df_point)) * se) ** 2

226 h_point = -(THETA_FDA**2) * rv.s2wr

227 h_bound = -(THETA_FDA**2) * rv.s2wr * rv.df / float(chi2.ppf(0.95, rv.df))

228 upper = e_point + h_point + math.sqrt((e_bound - e_point) ** 2 + (h_bound - h_point) ** 2)

229 return {

230 "criterion_point_estimate": e_point + h_point,

231 "criterion_95_upper_bound": upper,

232 "passes_scaled_criterion": float(upper <= 0.0),

233 }

234

235

236# ------------------------------------------------------------- power / N

237

238

239def tost_power(n_total: int, cv: float, gmr: float, design: str = "2x2", limits: tuple[float, float] = (0.80, 1.25)) -> float:

240 """Approximate TOST power by a finite quantile grid over estimated variance.

241

242 Equal allocation and normal log-endpoints are assumed. This is not an

243 exact Owen-Q routine; verify borderline planning decisions independently.

244 """

245 if design not in {"2x2", "parallel"}:

246 raise InputError("replicate power needs a specified sequence design and variance model; use PowerTOST")

247 if not math.isfinite(cv) or cv <= 0 or not math.isfinite(gmr) or gmr <= 0:

248 raise InputError("CV and GMR must be positive and finite")

249 if n_total < 4 or n_total % 2:

250 raise InputError("power calculation requires an even total N >= 4 for equal allocation")

251 sigma = math.sqrt(s2_from_cv(cv))

252 delta = math.log(gmr)

253 theta_low, theta_high = math.log(limits[0]), math.log(limits[1])

254

255 if design == "parallel":

256 df = n_total - 2

257 factor = math.sqrt(2.0 / (n_total / 2.0)) # equal arms

258 elif design == "2x2":

259 df = n_total - 2

260 factor = math.sqrt(2.0 / n_total)

261 if df < 1:

262 return 0.0

263

264 crit = float(t_dist.ppf(0.95, df))

265 # s^2 * df / sigma^2 ~ chi2_df

266 grid = np.linspace(1e-6, 1 - 1e-6, 2001)

267 chi_values = chi2.ppf(grid, df)

268 s_values = np.sqrt(chi_values / df) * sigma

269 se_values = s_values * factor

270 upper = (theta_high - crit * se_values - delta) / (sigma * factor)

271 lower = (theta_low + crit * se_values - delta) / (sigma * factor)

272 conditional = np.clip(norm.cdf(upper) - norm.cdf(lower), 0.0, 1.0)

273 return float(np.mean(conditional))

274

275

276def sample_size(cv: float, gmr: float, target: float, design: str = "2x2", limits: tuple[float, float] = (0.80, 1.25)) -> tuple[int, float]:

277 if not 0 < target < 1:

278 raise InputError("target power must be in (0,1)")

279 step = 2

280 for n in range(4, 5002, step):

281 power = tost_power(n, cv, gmr, design, limits)

282 if power >= target:

283 return n, power

284 raise InputError("no sample size below 5000 reaches the target power; the GMR is too far from 1")

285

286

287# --------------------------------------------------------------------- CLI

288

289

290def load_records(args: argparse.Namespace) -> list[dict]:

291 rows = read_table(args.input)

292 require_columns(rows, [args.subject_column, args.treatment_column, args.value_column], str(args.input))

293 records = []

294 for index, row in enumerate(rows, start=1):

295 treatment = row[args.treatment_column].strip().upper()

296 if treatment in {"T", "TEST"}:

297 treatment = "T"

298 elif treatment in {"R", "REF", "REFERENCE"}:

299 treatment = "R"

300 else:

301 raise InputError(f"row {index}: treatment must be T or R, got {row[args.treatment_column]!r}")

302 value = parse_float(row[args.value_column], f"{args.value_column} (row {index})")

303 if value is None or value <= 0:

304 raise InputError(f"row {index}: value must be positive to log-transform, got {value!r}")

305 records.append(

306 {

307 "subject": row[args.subject_column].strip(),

308 "treatment": treatment,

309 "value": value,

310 "logvalue": math.log(value),

311 "sequence": (row.get("sequence") or "").strip().upper(),

312 "period": (row.get("period") or "").strip(),

313 }

314 )

315 return records

316

317

318def build_parser() -> argparse.ArgumentParser:

319 parser = argparse.ArgumentParser(

320 description="Average, reference-scaled, and prospective bioequivalence calculations.",

321 formatter_class=argparse.RawDescriptionHelpFormatter,

322 epilog=__doc__,

323 )

324 parser.add_argument("-i", "--input")

325 parser.add_argument("--design", choices=("2x2", "parallel", "replicate"), default="2x2")

326 parser.add_argument("--metric", default="", help="label for the metric being tested, used in output only")

327 parser.add_argument("--subject-column", default="subject")

328 parser.add_argument("--treatment-column", default="treatment")

329 parser.add_argument("--value-column", default="value")

330 parser.add_argument(

331 "--scaling",

332 choices=("none", "abel", "rsabe", "both"),

333 default="none",

334 help="unsupported: requires external design-specific analysis",

335 )

336 parser.add_argument("--limits", default="0.80,1.25", help="acceptance limits for average BE (default: 0.80,1.25)")

337 parser.add_argument("--nti", action="store_true", help="narrow therapeutic index: apply 90.00-111.11%% limits")

338 parser.add_argument("--power", action="store_true", help="prospective power / sample size instead of an analysis")

339 parser.add_argument("--cv", type=float, help="assumed within-subject CV (as a fraction) for --power")

340 parser.add_argument("--gmr", type=float, default=0.95, help="assumed true GMR for --power (default: 0.95)")

341 parser.add_argument("--target-power", type=float, default=0.80)

342 parser.add_argument("--n", type=int, help="compute power at this N instead of solving for N")

343 add_format_argument(parser)

344 return parser

345

346

347def run(argv: Sequence[str] | None = None) -> int:

348 args = build_parser().parse_args(argv)

349 from _common import validate_numeric_args

350 validate_numeric_args(args)

351 report = Report()

352

353 limits = (0.90, 1.1111) if args.nti else tuple(parse_float(x, "acceptance limit") for x in args.limits.split(","))

354 if len(limits) != 2 or not all(math.isfinite(x) and x > 0 for x in limits) or limits[0] >= limits[1]:

355 raise InputError(f"--limits must be low,high with low < high; got {args.limits!r}")

356

357 if args.design == "replicate" or args.scaling != "none":

358 raise InputError("replicate/reference-scaling requires a design-specific validated analysis (FDA 2026 Appendix G, replicateBE); this CLI supports 2x2 and parallel ABE only")

359 if args.power:

360 if args.cv is None:

361 raise InputError("--power needs --cv")

362 if args.n:

363 power = tost_power(args.n, args.cv, args.gmr, args.design, limits)

364 report.scalar("n_total", args.n)

365 report.scalar("power", power)

366 else:

367 n, power = sample_size(args.cv, args.gmr, args.target_power, args.design, limits)

368 report.scalar("n_total_required", n)

369 report.scalar("achieved_power", power)

370 report.scalar("assumed_cv_within", args.cv)

371 report.scalar("assumed_gmr", args.gmr)

372 report.scalar("acceptance_limits", f"{limits[0]:.4f}-{limits[1]:.4f}")

373 report.note(

374 "Power is approximated by numerical integration over the sampling distribution of the estimated "

375 "standard deviation; the normal approximation overstates it at these sample sizes."

376 )

377 report.note(

378 "Planning depends jointly on GMR, CV, allocation and dropout. Assess plausible GMR/CV scenarios and inflate for dropout; the returned N is evaluable subjects."

379 )

380 return report.emit(args.format)

381

382 if not args.input:

383 raise InputError("give -i INPUT for an analysis, or --power for a sample-size calculation")

384

385 records = load_records(args)

386 label = args.metric or args.value_column

387

388 result = parallel_design(records) if args.design == "parallel" else crossover_2x2(records, args.design)

389 passes = limits[0] <= result.ci_low and result.ci_high <= limits[1]

390