# 10,000 simulated studies like make_data.py (10 mice per genotype, ages 4-20 weeks, the SAME true slope -0.8 per
# week, noise SD 6). How often is a "difference between genotypes" claimed (alpha 0.05) by
# (a) "significant in one genotype, not in the other" and (b) the interaction test for a difference in slopes?
import numpy as np
from scipy import stats
rng = np.random.default_rng(9)
runs, rule, inter = 10_000, 0, 0
for _ in range(runs):
    ages = [np.linspace(4, 20, 10) + rng.normal(0, 0.5, 10) for _ in range(2)]
    ys = [60 - 0.8 * (a - 4) + rng.normal(0, 6, 10) for a in ages]
    p = [stats.linregress(a, y).pvalue for a, y in zip(ages, ys)]
    rule += (p[0] < 0.05) != (p[1] < 0.05)
    x = np.concatenate(ages); y = np.concatenate(ys); g = np.r_[np.zeros(10), np.ones(10)]
    X = np.c_[np.ones(20), x, g, x * g]
    beta, res, *_ = np.linalg.lstsq(X, y, rcond=None)
    se = np.sqrt(res[0] / 16 * np.linalg.inv(X.T @ X)[3, 3])
    inter += 2 * stats.t.sf(abs(beta[3] / se), 16) < 0.05
print(f"runs = {runs}; difference claimed with identical true slopes: 'significant in one, not the other' {rule / runs:.1%}, interaction test {inter / runs:.1%}")