# Daily clinical scores, vehicle vs drug, 10 mice each, days 7-28. Python 3.13, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("scores.csv")                       # columns: mouse, group, day, score; one row per mouse per day
per = d.groupby(["group", "mouse"]).agg(cumulative=("score", "sum"), maximum=("score", "max")).reset_index()
onset = d[d.score >= 1].groupby(["group", "mouse"]).day.min()  # onset: first day with a score of 1 or more (0.5 does not count)
per["onset"] = [onset.get((r.group, r.mouse)) for r in per.itertuples()]
per["onset"] = per.onset.astype(float); per["diseased"] = per.onset.notna().astype(int)
per["onset_time"] = per.onset.fillna(d.day.max())   # mice without onset are censored at the last day
groups = ["vehicle", "drug"]; v, g = (per[per.group == k] for k in groups)

def mw(a, b): return stats.mannwhitneyu(b, a, method="asymptotic", use_continuity=True).pvalue
print(f"Incidence (score >= 1 on any day): vehicle {v.diseased.sum()}/10, drug {g.diseased.sum()}/10, "
      f"Fisher p = {stats.fisher_exact([[v.diseased.sum(), 10 - v.diseased.sum()], [g.diseased.sum(), 10 - g.diseased.sum()]]).pvalue:.4f}")
print(f"Cumulative score (sum of daily scores): median vehicle {v.cumulative.median():.1f}, drug {g.cumulative.median():.1f}; "
      f"Mann-Whitney p = {mw(v.cumulative, g.cumulative):.4f}")
q = lambda x: f"{x.quantile(0.25):.1f} to {x.quantile(0.75):.1f}"
print(f"Cumulative score, interquartile range: vehicle {q(v.cumulative)}, drug {q(g.cumulative)}")
hl = np.median(np.subtract.outer(g.cumulative.to_numpy(), v.cumulative.to_numpy()))
print(f"Cumulative score, Hodges-Lehmann shift (drug - vehicle, median of all pairwise differences): {hl:.1f}")
print(f"Maximum score: median vehicle {v.maximum.median():.2f}, drug {g.maximum.median():.2f}; Mann-Whitney p = {mw(v.maximum, g.maximum):.4f}")

t, e = per.onset_time.to_numpy(), per.diseased.to_numpy(); grp = (per.group == "drug").to_numpy()
times = np.unique(t[e == 1]); o1 = e1 = var = 0.0
for tt in times:                                     # log-rank test for time to onset
    r = t >= tt; r1 = (r & grp).sum(); rn = r.sum(); dd = (t == tt) & (e == 1); dn = dd.sum()
    o1 += (dd & grp).sum(); e1 += dn * r1 / rn
    if rn > 1: var += dn * (r1 / rn) * (1 - r1 / rn) * (rn - dn) / (rn - 1)
chi = (o1 - e1) ** 2 / var
print(f"Day of onset: median vehicle {v.onset.median():.1f}, drug {g.onset.median():.1f} (diseased mice); "
      f"log-rank chi-square = {chi:.2f}, p = {stats.chi2.sf(chi, 1):.4f}")

days = sorted(d.day.unique()); ps = []
for day in days:
    x = d[d.day == day]; a, b = x.score[x.group == "vehicle"], x.score[x.group == "drug"]
    ps.append(1.0 if len(set(a) | set(b)) == 1 else mw(a, b))   # all scores equal (e.g. all 0): no test
print(f"Separate Mann-Whitney test on each of {len(days)} days: p < 0.05 on {sum(p < 0.05 for p in ps)} days "
      f"(days {', '.join(str(dd) for dd, p in zip(days, ps) if p < 0.05)})")