# Comparing a marker's fluorescence between WT and KO: events as n, per-mouse arithmetic mean, per-mouse median on the
# log scale, and the per-mouse median with staining day as a block. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats

e = pd.read_csv("events.csv")                        # one row per gated event: day, mouse, genotype, fluorescence
ko_e = e.genotype == "KO"
t0 = stats.ttest_ind(e.fluorescence[ko_e], e.fluorescence[~ko_e], equal_var=False)
print(f"Events as n ({ko_e.sum()} vs {(~ko_e).sum()} events): Welch's t-test p = {t0.pvalue:.2g}")
m = e.groupby(["day", "mouse", "genotype"]).fluorescence.agg(["mean", "median"]).reset_index()
ko = (m.genotype == "KO").to_numpy()
t1 = stats.ttest_ind(m["mean"][ko], m["mean"][~ko], equal_var=False)
print(f"Per-mouse arithmetic mean (6 vs 6 mice): KO / WT {m['mean'][ko].mean() / m['mean'][~ko].mean():.2f}, Welch's t-test p = {t1.pvalue:.2g}")
lm = np.log(m["median"].to_numpy())
t2 = stats.ttest_ind(lm[ko], lm[~ko], equal_var=False); c2 = t2.confidence_interval()
print(f"Per-mouse median, log scale, ignoring day: KO / WT {np.exp(lm[ko].mean() - lm[~ko].mean()):.2f} "
      f"(95% CI {np.exp(c2.low):.2f} to {np.exp(c2.high):.2f}), p = {t2.pvalue:.2g}")
X = np.column_stack([np.ones(len(m)), ko.astype(float), (m.day == 2).astype(float), (m.day == 3).astype(float)])
b, *_ = np.linalg.lstsq(X, lm, rcond=None); res = lm - X @ b; dfe = len(lm) - X.shape[1]
se = np.sqrt(res @ res / dfe * np.linalg.inv(X.T @ X)[1, 1]); tq = stats.t.ppf(0.975, dfe); p = 2 * stats.t.sf(abs(b[1] / se), dfe)
print(f"Per-mouse median, log scale, staining day as a block: KO / WT {np.exp(b[1]):.2f} "
      f"(95% CI {np.exp(b[1] - tq * se):.2f} to {np.exp(b[1] + tq * se):.2f}), t({dfe}) p = {p:.2g}")
print("Per-day WT median fluorescence (mean of the two WT mice): " + "; ".join(f"day {dd} {m['median'][(~ko) & (m.day == dd)].mean():.0f}" for dd in (1, 2, 3)))