# Two treatments x two sexes, 5 mice per cell. Python 3.13, pandas, SciPy, statsmodels.
import numpy as np, pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm

d = pd.read_csv("grip.csv")
m = d.groupby(["sex", "treatment"]).grip_g.agg(["mean", "count"])
print("Cell means (g):")
for sex in ["female", "male"]:
    for trt in ["vehicle", "drug"]:
        print(f"  {sex:6s} {trt:7s} {m.loc[(sex, trt), 'mean']:.2f} (n = {m.loc[(sex, trt), 'count']})")

def g(sex, trt): return d.grip_g[(d.sex == sex) & (d.treatment == trt)]
print("Separate t-tests in each sex (the reading error):")
for sex in ["male", "female"]:
    t = stats.ttest_ind(g(sex, "drug"), g(sex, "vehicle"))
    print(f"  {sex}: difference {g(sex, 'drug').mean() - g(sex, 'vehicle').mean():.2f} g, t = {t.statistic:.2f}, p = {t.pvalue:.4f}")

t = stats.ttest_ind(d.grip_g[d.treatment == "drug"], d.grip_g[d.treatment == "vehicle"])
print(f"Pooled t-test ignoring sex: t = {t.statistic:.2f}, p = {t.pvalue:.4f}")

a = anova_lm(smf.ols("grip_g ~ C(treatment) * C(sex)", d).fit(), typ=2)   # Type II (all types agree here)
print("Two-way ANOVA (treatment x sex):")
for term, lab in [("C(treatment)", "treatment"), ("C(sex)", "sex"), ("C(treatment):C(sex)", "treatment x sex")]:
    print(f"  {lab:15s} F(1, {int(a.loc['Residual', 'df'])}) = {a.loc[term, 'F']:.2f}, p = {a.loc[term, 'PR(>F)']:.4f}")

# Contrasts on the four cell means, using the model's residual variance
n = 5; dfe = len(d) - 4; mse = a.loc["Residual", "sum_sq"] / a.loc["Residual", "df"]; tq = stats.t.ppf(0.975, dfe)
em = {(s, k): g(s, k).mean() for s in ["female", "male"] for k in ["vehicle", "drug"]}
eff = {s: em[(s, "drug")] - em[(s, "vehicle")] for s in ["female", "male"]}
rows = [("treatment effect, averaged over sexes", (eff["female"] + eff["male"]) / 2, np.sqrt(mse / n)),
        ("treatment effect in females", eff["female"], np.sqrt(2 * mse / n)),
        ("treatment effect in males", eff["male"], np.sqrt(2 * mse / n)),
        ("difference between sexes in the effect (male - female)", eff["male"] - eff["female"], np.sqrt(4 * mse / n))]
print("Effects with 95% CI (g):")
for lab, est, se in rows:
    print(f"  {lab}: {est:.2f} ({est - tq * se:.2f} to {est + tq * se:.2f})")