1"""Pilot power analysis (simulation) for the scenario-weighted paired contrast of acceptance probabilities.2Model: scenario random effect on the logit scale (sd = tau), arm effect delta (logit), R repeats, J scenarios; test = scenario-cluster3bootstrap percentile CI (approximated here by the paired t-test on scenario-level differences, which is conservative-equivalent at J = 10).4Usage: python power_simulation.py --out ../../analysis/tables/power_pilot.csv"""5from __future__ import annotations6import argparse, itertools, numpy as np, pandas as pd7from scipy import stats89def simulate(J, R, p_base, delta_pp, tau, n_sim=4000, seed=7, alpha=0.05):10 rng = np.random.default_rng(seed); rej = 011 logit = lambda p: np.log(p / (1 - p)); expit = lambda x: 1 / (1 + np.exp(-x))12 b0 = logit(p_base); d = logit(min(0.999, p_base + delta_pp / 100)) - b013 for _ in range(n_sim):14 u = rng.normal(0, tau, J)15 p0 = expit(b0 + u); p1 = expit(b0 + d + u)16 y0 = rng.binomial(R, p0) / R; y1 = rng.binomial(R, p1) / R17 diff = y1 - y018 if diff.std(ddof=1) == 0: rej += int(diff.mean() > 0); continue19 t = diff.mean() / (diff.std(ddof=1) / np.sqrt(J)); pval = 2 * (1 - stats.t.cdf(abs(t), J - 1))20 rej += int(pval < alpha and diff.mean() > 0)21 return rej / n_sim2223if __name__ == "__main__":24 ap = argparse.ArgumentParser(); ap.add_argument("--out", required=True); a = ap.parse_args()25 rows = []26 for J, R, p, dpp, tau in itertools.product([10, 20, 40], [2, 3, 5], [0.3, 0.5], [15, 25, 35], [0.8]):27 rows.append({"J_scenarios": J, "R_repeats": R, "p_base": p, "delta_pp": dpp, "tau_logit": tau, "power": simulate(J, R, p, dpp, tau)})28 df = pd.DataFrame(rows); df.to_csv(a.out, index=False); print(df.to_string())29