# 10,000 simulated randomised experiments per row: 12 mice per group, between-mouse SD 12 mmHg,
# baseline-week 4 correlation 0.7. Share of runs with p < 0.05 for each analysis.
import numpy as np
from scipy import stats
rng = np.random.default_rng(3)
def sim(effect, runs=10_000, n=12, r=0.7, sd=12.0):
    hits = np.zeros(4)
    g = np.r_[np.zeros(n), np.ones(n)]
    for _ in range(runs):
        base = 120 + rng.normal(0, sd, 2 * n)
        wk4 = 120 + r * (base - 120) + rng.normal(0, sd * np.sqrt(1 - r * r), 2 * n) + effect * g
        ch, pc = wk4 - base, 100 * (wk4 - base) / base
        for k, y in enumerate((wk4, ch, pc)):
            hits[k] += stats.ttest_ind(y[g == 1], y[g == 0], equal_var=False).pvalue < 0.05
        X = np.c_[np.ones(2 * n), g, base]
        beta, res, *_ = np.linalg.lstsq(X, wk4, rcond=None)
        se = np.sqrt(res[0] / (2 * n - 3) * np.linalg.inv(X.T @ X)[1, 1])
        hits[3] += 2 * stats.t.sf(abs(beta[1] / se), 2 * n - 3) < 0.05
    return hits / runs
for effect, label in [(0.0, "no true effect (false positives)"), (-10.0, "true effect -10 mmHg (power)")]:
    w, c, p, a = sim(effect)
    print(f"{label}: week 4 only {w:.1%}, change {c:.1%}, percent change {p:.1%}, ANCOVA {a:.1%}")