# Planned n = 8 per group; 4 more per group added after seeing p = 0.08. Python 3.13, pandas, SciPy, statsmodels.
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

d = pd.read_csv("rotarod.csv")                       # columns: mouse, group, cohort, rotarod_s
for lab, x in [("cohort 1", d[d.cohort == 1]), ("both cohorts", d)]:
    s = x.groupby("group").rotarod_s.agg(["mean", "std"])
    print(f"Mean (SD), {lab}: control {s.loc['control', 'mean']:.2f} ({s.loc['control', 'std']:.2f}), drug {s.loc['drug', 'mean']:.2f} ({s.loc['drug', 'std']:.2f})")
def test(x, welch=False):                            # Student's t-test (equal variances) unless welch=True
    a, b = x.rotarod_s[x.group == "control"], x.rotarod_s[x.group == "drug"]
    t = stats.ttest_ind(b, a, equal_var=not welch)
    return len(a), b.mean() - a.mean(), t.statistic, t.pvalue
n, diff, t, p = test(d[d.cohort == 1])
print(f"Planned analysis, cohort 1 only: n = {n} per group, difference {diff:.2f} s, t = {t:.2f}, p = {p:.4f}")
n, diff, t, p = test(d)
print(f"After adding cohort 2 (unplanned): n = {n} per group, difference {diff:.2f} s, t = {t:.2f}, p = {p:.4f}")
print(f"Welch's t-test instead: cohort 1 p = {test(d[d.cohort == 1], True)[3]:.4f}; both cohorts p = {test(d, True)[3]:.4f}")
fit = smf.ols("rotarod_s ~ C(group) + C(cohort)", d).fit()
est, ci, pv = fit.params["C(group)[T.drug]"], fit.conf_int().loc["C(group)[T.drug]"], fit.pvalues["C(group)[T.drug]"]
print(f"Both cohorts, cohort as a block: difference {est:.2f} s (95% CI {ci[0]:.2f} to {ci[1]:.2f}), p = {pv:.4f}")
n, diff, t, p = test(d[d.cohort == 2])
print(f"Cohort 2 alone: n = {n} per group, difference {diff:.2f} s, p = {p:.4f}")