# Single-dose cell viability as percent of vehicle control: wells or independent experiments as n, percent scale or
# log ratio. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("viability.csv")                     # one row per well: experiment, condition (blank/vehicle/treated), well, absorbance
d["signal"] = d.absorbance - d.experiment.map(d[d.condition == "blank"].groupby("experiment").absorbance.mean())   # blank-subtracted
w = d[d.condition != "blank"].copy()
veh = w[w.condition == "vehicle"].groupby("experiment").signal.mean(); trt = w[w.condition == "treated"].groupby("experiment").signal.mean()
w["pct"] = 100 * w.signal / w.experiment.map(veh)    # percent of the same experiment's mean vehicle signal
for e in veh.index:
    print(f"{e}: vehicle signal {veh[e]:.3f}, treated {trt[e]:.3f}, treated {100 * trt[e] / veh[e]:.1f}% of vehicle")
tw, vw = w.pct[w.condition == "treated"], w.pct[w.condition == "vehicle"]
a = stats.ttest_1samp(tw, 100)
print(f"A. Wells as n ({len(tw)} treated wells), one-sample t-test against 100%: mean {tw.mean():.1f}%, p = {a.pvalue:.2g}")
b = stats.ttest_ind(tw, vw, equal_var=False)
print(f"B. Wells as n ({len(tw)} vs {len(vw)}), Welch t-test on percent of vehicle: difference {tw.mean() - vw.mean():.1f} points, p = {b.pvalue:.2g}")
pe = 100 * trt / veh; c = stats.ttest_1samp(pe, 100); ci = c.confidence_interval()
print(f"C. Experiments as n ({len(pe)}), one-sample t-test of percent of vehicle against 100%: mean {pe.mean():.1f}% (95% CI {ci.low:.1f} to {ci.high:.1f}), p = {c.pvalue:.2g}")
assert (trt > 0).all() and (veh > 0).all(), "a blank-subtracted mean is 0 or below: check the plate before taking logs"
lr = np.log(trt / veh); t = stats.ttest_1samp(lr, 0); ci = t.confidence_interval()
print(f"D. Experiments as n ({len(lr)}), one-sample t-test of log(treated / vehicle): ratio {np.exp(lr.mean()):.3f} (95% CI {np.exp(ci.low):.3f} to {np.exp(ci.high):.3f}), p = {t.pvalue:.2g}")
p = stats.ttest_rel(trt, veh)
print(f"E. Experiments as n ({len(trt)}), paired t-test on blank-subtracted signal: difference {(trt - veh).mean():.3f}, p = {p.pvalue:.2g}")