⚙️ Setup & data¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try: df = pd.read_csv("../../data/biomedical_clinical_trials.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"biomedical_clinical_trials.csv")
print("loaded:", df.shape)
df.head()
loaded: (1000, 5)
| patient_id | treatment_group | baseline_systolic_bp | post_treatment_systolic_bp | bp_change_mmHg | |
|---|---|---|---|---|---|
| 0 | PT_1000 | Placebo | 140.5 | 131.3 | -9.2 |
| 1 | PT_1001 | Placebo | 140.4 | 142.7 | 2.3 |
| 2 | PT_1002 | Placebo | 119.4 | 117.3 | -2.1 |
| 3 | PT_1003 | Placebo | 129.6 | 134.6 | 5.0 |
| 4 | PT_1004 | Placebo | 131.4 | 133.3 | 1.9 |
active = df[df.treatment_group=="Active_Compound_X"]["bp_change_mmHg"]
placebo = df[df.treatment_group=="Placebo"]["bp_change_mmHg"]
print(f"treatment groups: {df.treatment_group.value_counts().to_dict()}")
tstat, pval = stats.ttest_1samp(active, popmean=0)
print(f"\nActive arm (n={len(active)}): mean bp change = {active.mean():.3f} mmHg, sd = {active.std():.3f}")
print(f"t = (xbar - 0)/(s/sqrt(n)) = {tstat:.2f}, df = {len(active)-1}, p-value = {pval:.2e}")
print("-> reject H0: BP dropped significantly in the treatment arm" if pval<0.05 else "-> not significant")
treatment groups: {'Placebo': 515, 'Active_Compound_X': 485}
Active arm (n=485): mean bp change = -2.580 mmHg, sd = 4.010
t = (xbar - 0)/(s/sqrt(n)) = -14.17, df = 484, p-value = 2.31e-38
-> reject H0: BP dropped significantly in the treatment arm
Blood pressure fell by about 2.6 mmHg in the treatment arm, and the one-sample t-test is overwhelmingly significant (p far below 0.05). The t-statistic, (x̄ − 0)/(s/√n), uses the sample standard deviation s, which is exactly why the t-distribution is the right reference. Case closed, the drug works? Not so fast.
tstat2, pval2 = stats.ttest_ind(active, placebo, equal_var=False)
print(f"Active arm: mean bp change = {active.mean():.3f} mmHg")
print(f"Placebo arm: mean bp change = {placebo.mean():.3f} mmHg")
print(f"difference (Active - Placebo) = {active.mean()-placebo.mean():+.3f} mmHg")
print(f"two-sample t = {tstat2:.3f}, p-value = {pval2:.3f}")
na, npl = len(active), len(placebo)
sp = np.sqrt(((na-1)*active.std(ddof=1)**2 + (npl-1)*placebo.std(ddof=1)**2)/(na+npl-2)) # pooled SD
print(f"Cohen d = {(active.mean()-placebo.mean())/sp:+.2f} (effect size; |0.2| small, |0.5| medium, |0.8| large)")
print("-> Active significantly beats Placebo" if pval2<0.05 else "-> NO significant difference: the drug adds nothing beyond placebo")
Active arm: mean bp change = -2.580 mmHg Placebo arm: mean bp change = -2.737 mmHg difference (Active - Placebo) = +0.157 mmHg two-sample t = 0.640, p-value = 0.522 Cohen d = +0.04 (effect size; |0.2| small, |0.5| medium, |0.8| large) -> NO significant difference: the drug adds nothing beyond placebo
fig,ax=plt.subplots(figsize=(6.5,3.2))
ax.boxplot([placebo,active],tick_labels=["Placebo","Active_Compound_X"],patch_artist=True,
boxprops=dict(facecolor="#fde68a"),medianprops=dict(color=PINK,linewidth=2))
ax.axhline(0,color=INK,ls=":",lw=1); ax.set_ylabel("BP change (mmHg)"); ax.set_title("BP change: Placebo vs Active (nearly identical)")
plt.tight_layout(); plt.show()
Here is the twist: the placebo arm dropped just as much (−2.74 vs −2.58 mmHg), and the two-sample t-test gives p ≈ 0.52, no significant difference. The apparent improvement in Beat 1 was real but not caused by the drug; blood pressure fell in both groups. This is the entire reason trials are placebo-controlled: testing against zero confuses a real-but-spurious drop with genuine efficacy, and only the comparison reveals the truth.
for n in [5,10,15,30,100]:
tcrit=stats.t.ppf(0.975,n-1); print(f"n={n:3d} (df={n-1:3d}): 95% t critical = {tcrit:.3f} vs normal z = 1.960")
n= 5 (df= 4): 95% t critical = 2.776 vs normal z = 1.960 n= 10 (df= 9): 95% t critical = 2.262 vs normal z = 1.960 n= 15 (df= 14): 95% t critical = 2.145 vs normal z = 1.960 n= 30 (df= 29): 95% t critical = 2.045 vs normal z = 1.960 n=100 (df= 99): 95% t critical = 1.984 vs normal z = 1.960
At n = 5 the 95% critical value is 2.78, far above the normal's 1.96, the t's way of demanding more evidence when σ is estimated from few points. As n grows the gap closes, and by n = 100 the t is essentially the normal. Using z on a tiny trial would overstate significance; the t keeps small-sample conclusions honest.