# Values below the limit of detection (LOD = 5 pg/mL): detected = 0 means "< 5" (recorded as 5).
# Python 3.13, NumPy, SciPy (left-censored likelihood written out, to check against R's survreg).
import numpy as np, pandas as pd
from scipy import stats, optimize

LOD = 5.0
d = pd.read_csv("cytokine.csv")
for k in ["control", "treated"]:
    x = d[d.group == k]; ld = np.log(x.conc_pg_ml[x.detected == 1])
    print(f"{k}: n = {len(x)}, below the LOD = {(x.detected == 0).sum()}, log-SD of detected values = {ld.std(ddof=1):.3f}")

sub = np.where(d.detected == 1, d.conc_pg_ml, LOD / 2)          # LOD/2 substitution
lc, lt = np.log(sub[d.group == "control"]), np.log(sub[d.group == "treated"])
t = stats.ttest_ind(lt, lc, equal_var=False)
print(f"LOD/2 substitution, Welch's t-test on log values: ratio of geometric means {np.exp(lt.mean() - lc.mean()):.2f}, "
      f"log-SD control {lc.std(ddof=1):.3f}, p = {t.pvalue:.4f}")
u = stats.mannwhitneyu(d.conc_pg_ml[d.group == "treated"], d.conc_pg_ml[d.group == "control"], method="asymptotic", use_continuity=True)
print(f"Mann-Whitney test (values below the LOD tied at the bottom): p = {u.pvalue:.4f}")

y = np.log(d.conc_pg_ml.to_numpy()); ev = d.detected.to_numpy(); xg = (d.group == "treated").to_numpy().astype(float)
def nll(q):                                           # log-normal model; a value below the LOD contributes P(X < LOD)
    z = (y - q[0] - q[1] * xg) / np.exp(q[2])
    return -(ev * (stats.norm.logpdf(z) - q[2]) + (1 - ev) * stats.norm.logcdf(z)).sum()
fit = optimize.minimize(nll, [np.log(8), 0, np.log(0.8)], method="BFGS")
E = np.eye(3) * 1e-4; p0 = fit.x
H = np.array([[(nll(p0 + a + b) - nll(p0 + a - b) - nll(p0 - a + b) + nll(p0 - a - b)) / 4e-8 for b in E] for a in E])
se = np.sqrt(np.diag(np.linalg.inv(H)))
b1, s1 = p0[1], se[1]
print(f"Left-censored log-normal model (maximum likelihood): ratio of geometric means {np.exp(b1):.2f} "
      f"(95% CI {np.exp(b1 - 1.96 * s1):.2f} to {np.exp(b1 + 1.96 * s1):.2f}), log-SD {np.exp(p0[2]):.3f}, p = {2 * stats.norm.sf(abs(b1 / s1)):.4f}")