# Tumour growth: vehicle vs drug, 8 mice each, measured every 3 days; mice culled at 1500 mm3.
# Python 3.13, pandas, NumPy, SciPy.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("tumours.csv")                       # columns: mouse, group, day, volume_mm3; one row per measurement
d["logv"] = np.log(d.volume_mm3)
groups = ["vehicle", "drug"]
print("Mice still measured on each day:")
for g in groups:
    c = d[d.group == g].groupby("day").mouse.count()
    print(f"  {g:7s} " + " ".join(f"d{day}:{n}" for day, n in c.items()))

def welch(a, b):                                       # difference b - a with Welch 95% CI and p
    t = stats.ttest_ind(b, a, equal_var=False)
    va, vb = a.var(ddof=1) / len(a), b.var(ddof=1) / len(b)
    df = (va + vb) ** 2 / (va ** 2 / (len(a) - 1) + vb ** 2 / (len(b) - 1))
    est = b.mean() - a.mean(); hw = stats.t.ppf(0.975, df) * np.sqrt(va + vb)
    return est, est - hw, est + hw, t.pvalue

last = d[d.day == 21]
a, b = last.logv[last.group == "vehicle"], last.logv[last.group == "drug"]
est, lo, hi, p = welch(a, b)
print(f"Day 21, surviving mice only (biased): n = {len(a)} vehicle, {len(b)} drug; "
      f"geometric means {np.exp(a.mean()):.0f} and {np.exp(b.mean()):.0f} mm3; T/C = {np.exp(est):.2f}, p = {p:.4f}")

print("Separate test on each day (log volume, Welch):")
ps = []
for day in sorted(d.day.unique()):
    x = d[d.day == day]
    ps.append(welch(x.logv[x.group == "vehicle"], x.logv[x.group == "drug"])[3])
print("  " + " ".join(f"d{day}:{p:.3f}" for day, p in zip(sorted(d.day.unique()), ps)))

slopes = d.groupby(["group", "mouse"]).apply(lambda m: np.polyfit(m.day, m.logv, 1)[0], include_groups=False)
sv, sd_ = slopes["vehicle"], slopes["drug"]
est, lo, hi, p = welch(sv, sd_)
print("Growth rate per mouse (slope of log volume per day, all measurements up to culling):")
print(f"  mean rate: vehicle {sv.mean():.4f}, drug {sd_.mean():.4f} per day; "
      f"doubling time {np.log(2) / sv.mean():.1f} and {np.log(2) / sd_.mean():.1f} days")
print(f"  difference (drug - vehicle) {est:.4f} per day (95% CI {lo:.4f} to {hi:.4f}), p = {p:.2g}; "
      f"rate ratio {sd_.mean() / sv.mean():.2f}")