# Histopathology scores (0-4), treated vs control: fields as n, a t-test on mouse scores, a rank test on mouse scores,
# a proportional-odds (ordinal) model, and the mean of the field scores per mouse. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats, optimize
d = pd.read_csv("scores.csv")                        # one row per field: group, mouse, field, score (0-4)
m = d.groupby(["group", "mouse"]).score.median().reset_index()        # one score per mouse: the median of its 9 fields
for g, x in m.groupby("group"):
    counts = ", ".join(f"{k}: {int((x.score == k).sum())}" for k in range(5))
    print(f"{g}: mouse median scores {counts}; median of mice {x.score.median():g}")
a, b = d.score[d.group == "treated"], d.score[d.group == "control"]
u = stats.mannwhitneyu(a, b, alternative="two-sided", method="asymptotic", use_continuity=True)
print(f"A. Fields as n ({len(a)} vs {len(b)}), Mann-Whitney test: p = {u.pvalue:.2g}")
mt, mc = m.score[m.group == "treated"], m.score[m.group == "control"]
t = stats.ttest_ind(mt, mc, equal_var=False)
print(f"B. Mouse scores treated as numbers, Welch t-test: mean {mt.mean():.2f} vs {mc.mean():.2f}, p = {t.pvalue:.2g}")
u = stats.mannwhitneyu(mt, mc, alternative="two-sided", method="asymptotic", use_continuity=True)
print(f"C. Mouse scores ({len(mt)} vs {len(mc)}), Mann-Whitney test: p = {u.pvalue:.2g}")
mf = d.groupby(["group", "mouse"]).score.mean().reset_index()
t = stats.ttest_ind(mf.score[mf.group == "treated"], mf.score[mf.group == "control"], equal_var=False)
print(f"D. Mean of the 9 field scores per mouse, Welch t-test: {mf.score[mf.group == 'treated'].mean():.2f} vs {mf.score[mf.group == 'control'].mean():.2f}, p = {t.pvalue:.2g}")
lev = np.sort(m.score.unique()); y = np.searchsorted(lev, m.score.to_numpy()); x = (m.group == "treated").to_numpy().astype(float); K = len(lev)
def nll(p, use_x=True):                              # proportional odds: logit P(Y <= k) = zeta_k - beta * treated
    zeta = np.cumsum(np.r_[p[0], np.exp(p[1:K - 1])]); eta = p[K - 1] * x if use_x else 0 * x
    cum = np.column_stack([np.zeros(len(y)), 1 / (1 + np.exp(-(zeta[None, :] - eta[:, None]))), np.ones(len(y))])
    return -np.sum(np.log(cum[np.arange(len(y)), y + 1] - cum[np.arange(len(y)), y]))
p0 = np.r_[-1.0, np.zeros(K - 2)]
f1 = optimize.minimize(nll, np.r_[p0, 0.0], method="BFGS", options={"gtol": 1e-9})
f0 = optimize.minimize(lambda p: nll(np.r_[p, 0.0], False), p0, method="BFGS", options={"gtol": 1e-9})
lr = 2 * (f0.fun - f1.fun)
print(f"E. Proportional-odds model of mouse scores: odds ratio for a higher score, treated vs control {np.exp(f1.x[-1]):.2f}, "
      f"likelihood-ratio chi2(1) = {lr:.2f}, p = {stats.chi2.sf(lr, 1):.2g}")