# How often each analysis gives p < 0.05: 5,000 simulated studies per row of 6 + 6 or 12 + 12 dams, every pup weighed
# at weaning, with no treatment effect and with treated pups 1 g (about 10%) lighter. Python 3.13, NumPy, SciPy. Seed 20261112.
# A first run used a 0.5 g effect (power 12.5% and 22.7% for litter means with 6 and 12 dams per group); the effect
# was changed to 1 g after that run so that differences in power between the analyses are visible.
import numpy as np
from scipy import stats
from model import study
rng = np.random.default_rng(20261112); R = 5_000

def one(dams, effect):
    rows = study(rng, dams, effect)
    tr = np.array([r[0] == "treated" for r in rows]); lit = np.array([r[1] for r in rows]); size = np.array([r[2] for r in rows])
    pup = np.array([r[3] for r in rows]); w = np.array([r[5] for r in rows])
    pa = stats.ttest_ind(w[tr], w[~tr], equal_var=False).pvalue
    L, inv = np.unique(lit, return_inverse=True); lm = np.bincount(inv, w) / np.bincount(inv)
    lt = np.bincount(inv, tr) > 0; ls = np.bincount(inv, size) / np.bincount(inv)
    pb = stats.ttest_ind(lm[lt], lm[~lt], equal_var=False).pvalue
    X = np.column_stack([np.ones(len(L)), lt, ls]).astype(float); b, *_ = np.linalg.lstsq(X, lm, rcond=None)
    res = lm - X @ b; df = len(L) - 3; se = np.sqrt(res @ res / df * np.linalg.inv(X.T @ X)[1, 1])
    pc = 2 * stats.t.sf(abs(b[1] / se), df)
    first = pup == 1
    pe = stats.ttest_ind(w[first & tr], w[first & ~tr], equal_var=False).pvalue
    return np.array([pa, pb, pc, pe]) < 0.05

for dams in (6, 12):
    for effect, lab in ((0.0, "No effect (false positives)"), (-1.0, "Treated pups 1 g lighter (power)")):
        r = np.array([one(dams, effect) for _ in range(R)]).mean(axis=0)
        print(f"{lab}, {dams} + {dams} dams: A pups as n {r[0]:.1%}; B litter means {r[1]:.1%}; C litter means adjusted for litter size {r[2]:.1%}; E one pup per litter {r[3]:.1%}")