import numpy as np, pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
d = pd.read_csv("pressure.csv")
d["group"] = pd.Categorical(d.group, ["vehicle", "drug"])
d["change"] = d.week4 - d.baseline
d["pct"] = 100 * d.change / d.baseline
fmt = lambda p: "< 0.001" if p < 0.001 else f"= {p:.3f}"
print(f"Baseline means: vehicle {d.baseline[d.group == 'vehicle'].mean():.1f}, drug {d.baseline[d.group == 'drug'].mean():.1f}; "
      f"baseline-week 4 correlation r = {np.corrcoef(d.baseline, d.week4)[0, 1]:.2f}")
for col, label in [("week4", "Week 4 only"), ("change", "Change from baseline"), ("pct", "Percent change")]:
    a, b = d.loc[d.group == "drug", col], d.loc[d.group == "vehicle", col]
    r = stats.ttest_ind(a, b, equal_var=False); ci = r.confidence_interval()
    print(f"{label}: drug - vehicle = {a.mean() - b.mean():.1f} (95% CI {ci.low:.1f} to {ci.high:.1f}), Welch t = {r.statistic:.2f}, df = {r.df:.1f}, p {fmt(r.pvalue)}")
m = smf.ols("week4 ~ group + baseline", d).fit()
lo, hi = m.conf_int().loc["group[T.drug]"]
print(f"ANCOVA (week 4 adjusted for baseline): drug - vehicle = {m.params['group[T.drug]']:.1f} (95% CI {lo:.1f} to {hi:.1f}), "
      f"t = {m.tvalues['group[T.drug]']:.2f}, df = {m.df_resid:.0f}, p {fmt(m.pvalues['group[T.drug]'])}")