# 10,000 simulated experiments per scenario. Grip strength SD 12 g; males 20 g stronger at baseline.
# Factorial design: 5 mice per sex per treatment (20 mice). Single-sex comparison: 10 males per treatment (20 mice).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
rng = np.random.default_rng(20261111); R = 10_000; n = 5; sd = 12; base = {"f": 100, "m": 120}

def sim(eff_f, eff_m):
    c = {(s, k): rng.normal(base[s] + (e if k == "d" else 0), sd, (R, n))
         for s, e in [("f", eff_f), ("m", eff_m)] for k in ["v", "d"]}
    mean = {key: x.mean(1) for key, x in c.items()}
    ss = sum(((x - x.mean(1, keepdims=True)) ** 2).sum(1) for x in c.values())
    dfe = 4 * n - 4; mse = ss / dfe
    ef, em = mean["f", "d"] - mean["f", "v"], mean["m", "d"] - mean["m", "v"]
    p_trt = 2 * stats.t.sf(abs(((ef + em) / 2) / np.sqrt(mse / n)), dfe)            # factorial: treatment main effect
    p_int = 2 * stats.t.sf(abs((em - ef) / np.sqrt(4 * mse / n)), dfe)              # factorial: treatment x sex
    p_sep = {s: stats.ttest_ind(c[s, "d"], c[s, "v"], axis=1).pvalue for s in ["f", "m"]}
    pooled = stats.ttest_ind(np.hstack([c["f", "d"], c["m", "d"]]), np.hstack([c["f", "v"], c["m", "v"]]), axis=1).pvalue
    one_only = ((p_sep["f"] < 0.05) != (p_sep["m"] < 0.05)).mean()
    return (p_trt < 0.05).mean(), (pooled < 0.05).mean(), (p_int < 0.05).mean(), one_only

def single_sex(eff):
    v, d_ = rng.normal(120, sd, (R, 2 * n)), rng.normal(120 + eff, sd, (R, 2 * n))
    return (stats.ttest_ind(d_, v, axis=1).pvalue < 0.05).mean()

for label, ef, em in [("No treatment effect", 0, 0), ("Same effect in both sexes (12 g)", 12, 12),
                      ("Effect in males only (24 g; average 12 g)", 0, 24)]:
    pt, pp, pi, oo = sim(ef, em)
    print(f"{label}: factorial treatment p<0.05 {pt:.1%}; pooled t-test ignoring sex {pp:.1%}; "
          f"interaction p<0.05 {pi:.1%}; significant in exactly one sex {oo:.1%}")
print(f"Males only, 10 per treatment, effect 12 g: t-test p<0.05 {single_sex(12):.1%}")