import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("doses.csv")
groups = [d.response[d.level == k].to_numpy() for k in range(4)]
means = np.array([g.mean() for g in groups]); n = np.array([len(g) for g in groups])
print("Group means (vehicle, 1, 3, 10 mg/kg): " + ", ".join(f"{m:.2f}" for m in means))
# Each dose against vehicle, Dunnett's adjustment (pooled variance)
dn = stats.dunnett(*groups[1:], control=groups[0], rng=np.random.default_rng(1))
print("Dunnett, each dose vs vehicle: p = " + ", ".join(f"{p:.2f}" for p in dn.pvalue))
# Linear trend contrast across the four ordered groups (weights -3, -1, 1, 3), within-group variance
c = np.array([-3, -1, 1, 3])
df = int(n.sum() - len(n)); mse = sum(((g - g.mean()) ** 2).sum() for g in groups) / df
est = (c * means).sum(); se = np.sqrt(mse * (c ** 2 / n).sum())
t = est / se; p = 2 * stats.t.sf(abs(t), df)
print(f"Linear trend contrast: estimate = {est:.2f}, t = {t:.2f}, df = {df}, p = {p:.4f}")
# Change per step of the ordering (scores 0, 1, 2, 3): contrast divided by sum(c * score), with its 95% CI
scores = np.array([0, 1, 2, 3]); k = (c * scores).sum()
tq = stats.t.ppf(0.975, df)
print(f"Change per dose step (scores 0-3): {est / k:.2f} (95% CI {(est - tq * se) / k:.2f} to {(est + tq * se) / k:.2f})")