# Bacterial burden (CFU per spleen), treated vs control: raw CFU, log10 CFU with values below the limit of detection
# set to LOD/2, a rank test, and a censored (Tobit) model. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats, optimize
d = pd.read_csv("cfu.csv")                           # one row per mouse: group, mouse, dilution (10^-d plated), colonies
VOL_ML, HOMOG_ML = 0.1, 1.0                          # 100 uL plated; organ homogenised in 1 mL
d["cfu"] = d.colonies * 10.0 ** d.dilution / VOL_ML * HOMOG_ML
d["lod"] = 10.0 ** d.dilution / VOL_ML * HOMOG_ML  # one colony on the plate read: 10 CFU per organ on the undiluted plate
d["below"] = d.colonies == 0; LOD = d.lod[d.below].max() if d.below.any() else d.lod.min()
d["log_half"] = np.log10(np.where(d.below, d.lod / 2, d.cfu))
for g, x in d.groupby("group"):
    det = x.cfu[~x.below]
    print(f"{g}: {len(x)} mice, {x.below.sum()} below the limit of detection ({LOD:g} CFU); detected mice range {det.min():g} to {det.max():g} CFU")
c, t_ = d.group == "control", d.group == "treated"
t = stats.ttest_ind(d.cfu[t_], d.cfu[c], equal_var=False)
print(f"A. Raw CFU, Welch t-test: mean treated {d.cfu[t_].mean():.0f}, control {d.cfu[c].mean():.0f}, p = {t.pvalue:.2g}")
t = stats.ttest_ind(d.log_half[t_], d.log_half[c], equal_var=False); ci = t.confidence_interval()
print(f"B. log10 CFU, below-LOD values set to LOD/2, Welch t-test: difference {d.log_half[t_].mean() - d.log_half[c].mean():.2f} log10 (95% CI {ci.low:.2f} to {ci.high:.2f}), p = {t.pvalue:.2g}")
u = stats.mannwhitneyu(d.cfu[t_], d.cfu[c], alternative="two-sided", method="asymptotic", use_continuity=True)
print(f"C. Mann-Whitney test on CFU (below-LOD values tied at the bottom): p = {u.pvalue:.2g}")
y = np.log10(np.where(d.below, d.lod, d.cfu)); cens = d.below.to_numpy(); x = t_.to_numpy().astype(float)
def nll(p, use_x=True):                              # left-censored normal (Tobit) likelihood
    mu = p[0] + (p[1] * x if use_x else 0); s = np.exp(p[-1])
    return -(np.sum(stats.norm.logpdf(y[~cens], mu[~cens] if use_x else mu, s)) + np.sum(stats.norm.logcdf((y[cens] - (mu[cens] if use_x else mu)) / s)))
f1 = optimize.minimize(nll, [y.mean(), 0, np.log(y.std())], method="BFGS", options={"gtol": 1e-10})
f0 = optimize.minimize(lambda p: nll(p, False), [y.mean(), np.log(y.std())], method="BFGS", options={"gtol": 1e-10})
lr = 2 * (f0.fun - f1.fun)
h = 1e-4; k = len(f1.x); H = np.zeros((k, k))                # observed information by central differences
for i in range(k):
    for j in range(k):
        e_i, e_j = np.eye(k)[i] * h, np.eye(k)[j] * h
        H[i, j] = (nll(f1.x + e_i + e_j) - nll(f1.x + e_i - e_j) - nll(f1.x - e_i + e_j) + nll(f1.x - e_i - e_j)) / (4 * h * h)
se = np.sqrt(np.linalg.inv(H)[1, 1])
print(f"D. Censored (Tobit) model of log10 CFU, LOD treated as censoring: difference {f1.x[1]:.2f} log10 (95% Wald CI {f1.x[1] - 1.96 * se:.2f} to {f1.x[1] + 1.96 * se:.2f}), "
      f"likelihood-ratio chi2(1) = {lr:.2f}, p = {stats.chi2.sf(lr, 1):.2g}")