# Litter effects: a treatment given to pregnant dams, pups weighed at weaning. Pups as n, litter means, and litter
# means adjusted for litter size. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("pups.csv")                          # one row per pup: treatment, litter, litter_size, pup, sex, weight_g
m = d.groupby(["treatment", "litter"]).agg(weight_g=("weight_g", "mean"), litter_size=("litter_size", "first")).reset_index()
for g in ("control", "treated"):
    x = m[m.treatment == g]
    print(f"{g}: {len(x)} litters, {int((d.treatment == g).sum())} pups, mean litter size {x.litter_size.mean():.1f}, mean of litter means {x.weight_g.mean():.2f} g")
t = stats.ttest_ind(d.weight_g[d.treatment == "treated"], d.weight_g[d.treatment == "control"], equal_var=False)
print(f"A. Pups as n ({int((d.treatment == 'treated').sum())} vs {int((d.treatment == 'control').sum())}), Welch t-test: difference {d.weight_g[d.treatment == 'treated'].mean() - d.weight_g[d.treatment == 'control'].mean():.2f} g, p = {t.pvalue:.2g}")
tr, co = m.weight_g[m.treatment == "treated"], m.weight_g[m.treatment == "control"]
t = stats.ttest_ind(tr, co, equal_var=False); ci = t.confidence_interval()
print(f"B. Litters as n ({len(tr)} vs {len(co)}), Welch t-test on litter means: difference {tr.mean() - co.mean():.2f} g (95% CI {ci.low:.2f} to {ci.high:.2f}), p = {t.pvalue:.2g}")
X = np.column_stack([np.ones(len(m)), m.treatment == "treated", m.litter_size]).astype(float); y = m.weight_g.to_numpy()
b, *_ = np.linalg.lstsq(X, y, rcond=None); res = y - X @ b; df = len(y) - 3
se = np.sqrt(res @ res / df * np.linalg.inv(X.T @ X)[1, 1]); q = stats.t.ppf(0.975, df)
print(f"C. Litters as n, litter means adjusted for litter size (linear model): difference {b[1]:.2f} g (95% CI {b[1] - q * se:.2f} to {b[1] + q * se:.2f}), p = {2 * stats.t.sf(abs(b[1] / se), df):.2g}")