Skills › Scientific Agent Skills › pkpd-modeling › exposure_response.py
exposure_response.py
353 lines · skills/pkpd-modeling/scripts/exposure_response.py at 5f2bc03
1#!/usr/bin/env python3
2"""Exposure-response analysis: Emax, logistic, concentration-QTc, and exposure quartiles.
4Exposure-response is where dose selection is actually decided, and where the
5most consequential statistical mistakes are made. Two dominate. The first is
6fitting Emax to data that never approached the plateau, which produces an Emax
7and EC50 that are individually meaningless but jointly reproduce the observed
8slope. The second is reading an exposure-response relationship causally when
9exposure is itself a consequence of the patient's condition — sicker patients
10clear drug differently, so a flat or inverted E-R curve can be confounding
11rather than pharmacology.
13 python3 exposure_response.py --emax -i er.csv
14 python3 exposure_response.py --logistic -i er.csv --response-column responder
15 python3 exposure_response.py --cqtc -i qt.csv --cmax 250
16 python3 exposure_response.py --quartiles -i er.csv --response-column response
18Default columns are ``exposure`` and ``response``.
19"""
21from __future__ import annotations
23import argparse
24import math
25from typing import Sequence
27import numpy as np
29from _common import (
30 InputError,
31 Report,
32 add_format_argument,
33 main_wrapper,
34 parse_float,
35 read_table,
36 require_columns,
37)
39try:
40 from scipy.optimize import least_squares, minimize
41 from scipy.stats import norm
42 from scipy.stats import t as t_dist
43except ImportError as exc: # pragma: no cover
44 raise SystemExit("exposure_response.py needs scipy: uv pip install scipy") from exc
46# ICH E14's threshold of regulatory concern: an upper bound of the two-sided
47# 90% confidence interval for placebo-corrected change-from-baseline QTc above
48# 10 ms.
49QTC_THRESHOLD_MS = 10.0
52def fit_emax(exposure: np.ndarray, response: np.ndarray, sigmoid: bool) -> dict:
53 """Fit E0 + Emax*C^h/(EC50^h + C^h) by least squares on the log of positive parameters."""
54 if len(exposure) <= (4 if sigmoid else 3) or np.ptp(exposure) <= 0:
55 raise InputError("Emax needs more observations than parameters and varied exposures")
56 e0_init = float(np.min(response))
57 emax_init = float(np.max(response) - np.min(response)) or 1.0
58 ec50_init = float(np.median(exposure[exposure > 0])) if np.any(exposure > 0) else 1.0
60 def unpack(theta: np.ndarray) -> tuple[float, float, float, float]:
61 e0 = theta[0]
62 emax_value = theta[1]
63 ec50 = math.exp(theta[2])
64 hill = math.exp(theta[3]) if sigmoid else 1.0
65 return e0, emax_value, ec50, hill
67 def predict(theta: np.ndarray, x: np.ndarray) -> np.ndarray:
68 e0, emax_value, ec50, hill = unpack(theta)
69 # Logistic evaluation avoids overflow of C**hill and EC50**hill.
70 z = hill * (np.log(np.maximum(x, np.finfo(float).tiny)) - math.log(ec50))
71 fraction = np.exp(-np.logaddexp(0.0, -z))
72 return e0 + emax_value * np.where(x > 0, fraction, 0.0)
74 start = [e0_init, emax_init, math.log(max(ec50_init, 1e-9))] + ([0.0] if sigmoid else [])
75 result = least_squares(lambda th: predict(th, exposure) - response, start, max_nfev=20000)
76 e0, emax_value, ec50, hill = unpack(result.x)
78 n, p = len(exposure), len(result.x)
79 residuals = response - predict(result.x, exposure)
80 ssr = float(np.sum(residuals**2))
81 sigma2 = ssr / max(n - p, 1)
82 try:
83 _, s, vt = np.linalg.svd(result.jac, full_matrices=False)
84 threshold = np.finfo(float).eps * max(result.jac.shape) * s[0]
85 if np.count_nonzero(s > threshold) < p:
86 raise np.linalg.LinAlgError("rank-deficient Emax sensitivity")
87 s_inv = 1.0 / s
88 cov = (vt.T * s_inv**2) @ vt * sigma2
89 se = np.sqrt(np.maximum(np.diag(cov), 0.0))
90 except np.linalg.LinAlgError: # pragma: no cover
91 se = np.full(p, np.nan)
93 max_observed = float(np.max(exposure))
94 return {
95 "converged": bool(result.success),
96 "identifiable_local_covariance": bool(np.all(np.isfinite(se))),
97 "e0": e0,
98 "emax": emax_value,
99 "ec50": ec50,
100 "hill": hill,
101 "se_e0": float(se[0]),
102 "se_emax": float(se[1]),
103 "rse_ec50_pct": 100.0 * float(se[2]),
104 "rse_hill_pct": 100.0 * float(se[3]) if sigmoid else float("nan"),
105 "residual_sd": math.sqrt(sigma2),
106 "n": n,
107 "max_observed_exposure": max_observed,
108 "ec50_over_max_exposure": ec50 / max_observed if max_observed else float("nan"),
109 "predicted_effect_at_max_exposure": float(predict(result.x, np.array([max_observed]))[0]),
110 "fraction_of_emax_reached": float(
111 (predict(result.x, np.array([max_observed]))[0] - e0) / emax_value
112 )
113 if emax_value
114 else float("nan"),
115 }
118def fit_logistic(exposure: np.ndarray, response: np.ndarray) -> dict:
119 """Maximum-likelihood logistic regression of a binary response on exposure."""
120 unique = set(np.unique(response).tolist())
121 if not unique <= {0.0, 1.0}:
122 raise InputError("logistic exposure-response needs a 0/1 response column")
124 if unique != {0.0, 1.0} or np.ptp(exposure) <= 0:
125 raise InputError("logistic regression needs both outcomes and varied exposures")
126 zero, one = exposure[response == 0], exposure[response == 1]
127 if zero.max() <= one.min() or one.max() <= zero.min():
128 raise InputError("complete or quasi-complete separation: unpenalized logistic slope is not finite")
130 def neg_loglik(theta: np.ndarray) -> float:
131 eta = theta[0] + theta[1] * exposure
132 # log(1+exp(eta)) computed stably
133 return float(np.sum(np.logaddexp(0.0, eta) - response * eta))
135 result = minimize(neg_loglik, np.array([0.0, 0.0]), method="BFGS")
136 intercept, slope = result.x
137 try:
138 cov = result.hess_inv
139 se = np.sqrt(np.maximum(np.diag(np.asarray(cov)), 0.0))
140 except (ValueError, TypeError): # pragma: no cover
141 se = np.array([float("nan"), float("nan")])
143 ec50 = -intercept / slope if slope else float("nan")
144 return {
145 "converged": bool(result.success),
146 "intercept": float(intercept),
147 "slope": float(slope),
148 "se_slope": float(se[1]),
149 "slope_z": float(slope / se[1]) if se[1] else float("nan"),
150 "slope_p": float(2 * (1 - norm.cdf(abs(slope / se[1])))) if se[1] else float("nan"),
151 "odds_ratio_per_unit": float(math.exp(slope)),
152 "exposure_at_50pct_probability": float(ec50),
153 "n": len(exposure),
154 "n_responders": int(response.sum()),
155 "neg_loglik": float(result.fun),
156 }
159def fit_cqtc(conc: np.ndarray, delta_qtc: np.ndarray, cmax: float | None) -> dict:
160 """Linear concentration-QTc model with the two-sided 90% CI of the predicted effect.
162 ICH E14 asks whether the **upper bound of the two-sided 90% confidence
163 interval** for placebo-corrected change-from-baseline QTc is below 10 ms at
164 the clinically relevant exposure. Reporting the point estimate, or a 95%
165 interval, answers a different question than the guideline asks.
166 """
167 n = len(conc)
168 if n < 4:
169 raise InputError("concentration-QTc analysis needs at least 4 observations")
170 if np.ptp(conc) <= 0:
171 raise InputError("concentration-QTc needs varied concentrations")
172 design = np.column_stack([np.ones(n), conc])
173 coefficients, *_ = np.linalg.lstsq(design, delta_qtc, rcond=None)
174 fitted = design @ coefficients
175 residuals = delta_qtc - fitted
176 dof = n - 2
177 sigma2 = float(residuals @ residuals) / dof
178 xtx_inv = np.linalg.inv(design.T @ design)
179 se = np.sqrt(np.diag(xtx_inv) * sigma2)
180 crit = float(t_dist.ppf(0.95, dof)) # two-sided 90%
182 out = {
183 "intercept_ms": float(coefficients[0]),
184 "slope_ms_per_conc": float(coefficients[1]),
185 "se_slope": float(se[1]),
186 "slope_ci90_low": float(coefficients[1] - crit * se[1]),
187 "slope_ci90_high": float(coefficients[1] + crit * se[1]),
188 "residual_sd_ms": math.sqrt(sigma2),
189 "n": n,
190 "degrees_of_freedom": dof,
191 }
192 if cmax is not None:
193 x0 = np.array([1.0, cmax])
194 prediction = float(x0 @ coefficients)
195 se_mean = math.sqrt(float(x0 @ xtx_inv @ x0) * sigma2)
196 out.update(
197 {
198 "exposure_evaluated": cmax,
199 "predicted_delta_delta_qtc_ms": prediction,
200 "ci90_low_ms": prediction - crit * se_mean,
201 "ci90_high_ms": prediction + crit * se_mean,
202 "upper_bound_exceeds_10ms": float((prediction + crit * se_mean) >= QTC_THRESHOLD_MS),
203 }
204 )
205 if cmax > float(np.max(conc)) or cmax < float(np.min(conc)):
206 out["extrapolated_beyond_observed"] = 1.0
207 return out
210def quartile_summary(exposure: np.ndarray, response: np.ndarray, n_bins: int) -> list[dict]:
211 edges = np.quantile(exposure, np.linspace(0, 1, n_bins + 1))
212 edges[-1] = np.nextafter(edges[-1], np.inf)
213 rows = []
214 for i in range(n_bins):
215 mask = (exposure >= edges[i]) & (exposure < edges[i + 1])
216 if not np.any(mask):
217 continue
218 values = response[mask]
219 rows.append(
220 {
221 "bin": i + 1,
222 "exposure_range": f"{edges[i]:.4g} - {edges[i + 1]:.4g}",
223 "n": int(mask.sum()),
224 "median_exposure": float(np.median(exposure[mask])),
225 "mean_response": float(values.mean()),
226 "sd_response": float(values.std(ddof=1)) if len(values) > 1 else float("nan"),
227 "responder_fraction": float(values.mean()) if set(np.unique(values).tolist()) <= {0.0, 1.0} else float("nan"),
228 }
229 )
230 return rows
233# --------------------------------------------------------------------- CLI
236def build_parser() -> argparse.ArgumentParser:
237 parser = argparse.ArgumentParser(
238 description="Exposure-response modelling: Emax, logistic, concentration-QTc, quartiles.",
239 formatter_class=argparse.RawDescriptionHelpFormatter,
240 epilog=__doc__,
241 )
242 parser.add_argument("--emax", action="store_true")
243 parser.add_argument("--sigmoid", action="store_true", help="estimate the Hill coefficient too")
244 parser.add_argument("--logistic", action="store_true")
245 parser.add_argument("--cqtc", action="store_true")
246 parser.add_argument("--quartiles", action="store_true")
247 parser.add_argument("-i", "--input", required=True)
248 parser.add_argument("--exposure-column", default="exposure")
249 parser.add_argument("--response-column", default="response")
250 parser.add_argument("--cmax", type=float, help="exposure at which to evaluate the C-QTc prediction")
251 parser.add_argument("--bins", type=int, default=4)
252 add_format_argument(parser)
253 return parser
256def run(argv: Sequence[str] | None = None) -> int:
257 args = build_parser().parse_args(argv)
258 from _common import validate_numeric_args
259 validate_numeric_args(args)
260 if args.bins < 1 or (args.cmax is not None and args.cmax < 0):
261 raise InputError("--bins must be positive and --cmax non-negative")
262 modes = [m for m in ("emax", "logistic", "cqtc", "quartiles") if getattr(args, m)]
263 if len(modes) != 1:
264 raise InputError("choose exactly one of --emax, --logistic, --cqtc, --quartiles")
266 rows = read_table(args.input)
267 require_columns(rows, [args.exposure_column, args.response_column], str(args.input))
268 exposure = np.asarray([parse_float(r[args.exposure_column], "exposure") for r in rows], dtype=float)
269 response = np.asarray([parse_float(r[args.response_column], "response") for r in rows], dtype=float)
270 if np.any(exposure < 0):
271 raise InputError("exposure must not be negative")
273 report = Report()
274 report.note(
275 "Exposure-response is observational even inside a randomised trial: patients are randomised to "
276 "dose, not to exposure. Differences across exposure quantiles can reflect the covariates that "
277 "drive clearance rather than the drug."
278 )
280 if args.emax:
281 fit = fit_emax(exposure, response, args.sigmoid)
282 report.table("Emax model", [{"parameter": k, "value": v} for k, v in fit.items()])
283 if fit["fraction_of_emax_reached"] < 0.5:
284 report.finding(
285 f"the highest observed exposure reaches only {100 * fit['fraction_of_emax_reached']:.0f}% "
286 "of the estimated Emax. Emax and EC50 are then extrapolations, strongly correlated with "
287 "each other, and should not be quoted as independent estimates."
288 )
289 if fit["ec50_over_max_exposure"] > 1.0:
290 report.finding(
291 f"EC50 ({fit['ec50']:.4g}) exceeds the highest observed exposure "
292 f"({fit['max_observed_exposure']:.4g}); the plateau is entirely outside the data"
293 )
294 if math.isfinite(fit["rse_ec50_pct"]) and fit["rse_ec50_pct"] > 50:
295 report.finding(f"EC50 has {fit['rse_ec50_pct']:.0f}% relative standard error")
296 report.note("a linear-looking E-R relationship is the low-concentration limb of an Emax curve; the two are not distinguishable without data near the plateau")
298 elif args.logistic:
299 fit = fit_logistic(exposure, response)
300 report.table("logistic exposure-response", [{"parameter": k, "value": v} for k, v in fit.items()])
301 if fit["n_responders"] < 10 or fit["n"] - fit["n_responders"] < 10:
302 report.finding(
303 f"{fit['n_responders']} responders out of {fit['n']}; with fewer than about 10 events "
304 "per covariate the slope estimate is unstable and its confidence interval unreliable"
305 )
306 if math.isfinite(fit["slope_p"]) and fit["slope_p"] > 0.05:
307 report.note(f"the exposure slope is not statistically distinguishable from zero (p = {fit['slope_p']:.3f})")
309 elif args.cqtc:
310 fit = fit_cqtc(exposure, response, args.cmax)
311 report.table("concentration-QTc", [{"parameter": k, "value": v} for k, v in fit.items()])
312 report.note(
313 "ICH E14 evaluates the upper bound of the two-sided 90% confidence interval for "
314 "placebo-corrected change-from-baseline QTc against a 10 ms threshold. The response column "
315 "must already be that placebo-corrected change; this script does not compute it."
316 )
317 if fit.get("upper_bound_exceeds_10ms"):
318 report.finding(
319 f"the 90% upper bound at the evaluated exposure is {fit['ci90_high_ms']:.2f} ms, above "
320 "the 10 ms threshold of regulatory concern"
321 )
322 if fit.get("extrapolated_beyond_observed"):
323 report.finding(
324 f"the evaluated exposure ({fit['exposure_evaluated']:.4g}) is above the highest observed "
325 "concentration; the prediction is an extrapolation and the interval understates its uncertainty"
326 )
327 report.note(
328 "This is an ordinary linear model. A regulatory C-QTc analysis uses a mixed model with a "
329 "random intercept and slope per subject and a treatment-specific intercept; use this for "
330 "screening and exploration, not for submission."
331 )
333 else:
334 rows_out = quartile_summary(exposure, response, args.bins)
335 report.table(f"exposure {args.bins}-quantile summary", rows_out)
336 if len(rows_out) >= 2:
337 first, last = rows_out[0]["mean_response"], rows_out[-1]["mean_response"]
338 report.scalar("response_change_across_range", last - first)
339 report.note(
340 "Quantile summaries are descriptive. They do not adjust for the covariates that determine "
341 "exposure, and a monotone trend across quantiles is not by itself evidence of a causal "
342 "dose-response relationship."
343 )
345 if (args.emax or args.logistic) and not fit["converged"]:
346 report.finding("optimizer did not converge; estimates and Wald intervals are unreliable")
347 if args.emax and not fit["identifiable_local_covariance"]:
348 report.finding("rank-deficient Emax sensitivity; covariance unavailable, use profile likelihood or a simpler model")
349 return report.emit(args.format)
352if __name__ == "__main__":
353 raise SystemExit(main_wrapper(run))