⚙️ Setup¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
# statsmodels = the R/SAS-style stats library (pre-installed on Colab): it computes the
# standard errors, test statistics, intervals, and post-hoc comparisons, so we write less by hand.
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.proportion import proportions_ztest, confint_proportions_2indep, proportion_confint
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.weightstats import DescrStatsW, CompareMeans
from scipy.stats.contingency import association # Cramer's V in one call
IND="#4f46e5"; DEEP="#4338ca"; LIGHT="#818cf8"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"
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/"
rng = np.random.default_rng(75)
alpha=0.05; trials=20000
pvals=np.array([stats.ttest_1samp(rng.normal(0,1,40), 0).pvalue for _ in range(trials)])
print(f"H0 is TRUE in every trial. Rejections at alpha={alpha}: {np.mean(pvals<alpha):.3f} (should be ~0.05)")
fig,ax=plt.subplots(figsize=(7,3.1))
ax.hist(pvals,bins=20,color=LIGHT,alpha=0.85,edgecolor="white")
ax.axhline(trials/20,color=DEEP,ls="--",lw=1.5,label="uniform")
ax.set_title("Under H0, p-values are UNIFORM on [0,1]"); ax.set_xlabel("p-value"); ax.legend()
plt.tight_layout(); plt.show()
H0 is TRUE in every trial. Rejections at alpha=0.05: 0.050 (should be ~0.05)
When H0 is true the p-value is uniform, every value is equally likely, so exactly 5% land below 0.05 purely by chance. That 5% is the Type I error rate you sign up for. Significance is a controlled false-alarm budget, not a guarantee of truth.
def power_sim(true_mean, n, alpha=0.05, B=4000):
return np.mean([stats.ttest_1samp(rng.normal(true_mean,1,n),0).pvalue<alpha for _ in range(B)])
print("Power to detect a true mean of 0.5 (sd 1):")
for n in [10,20,40,80,160]:
print(f" n={n:>3}: power = {power_sim(0.5,n):.2f}")
print("\nPower vs effect size (n=40):")
for eff in [0.2,0.35,0.5,0.8]:
print(f" effect d={eff}: power = {power_sim(eff,40):.2f}")
Power to detect a true mean of 0.5 (sd 1):
n= 10: power = 0.30 n= 20: power = 0.57 n= 40: power = 0.86 n= 80: power = 0.99
n=160: power = 1.00 Power vs effect size (n=40): effect d=0.2: power = 0.24
effect d=0.35: power = 0.58 effect d=0.5: power = 0.86
effect d=0.8: power = 1.00
Bigger samples and bigger effects both raise power. A small study can easily miss a real effect, that is a Type II error, and "we found no significant difference" from an underpowered test means little. Designing for adequate power (often 80%) is as important as controlling alpha.
# Trap 1: multiple comparisons on pure noise
sig=sum(stats.ttest_ind(rng.normal(0,1,50),rng.normal(0,1,50)).pvalue<0.05 for _ in range(20))
print(f"20 tests, all H0 true -> {sig} came out significant (expected ~1) -> false positives multiply")
# Trap 2: significance != importance
big=rng.normal(0.02,1,100000) # a tiny true effect of 0.02
p=stats.ttest_1samp(big,0).pvalue
print(f"\nn=100000, true effect only 0.02: p={p:.2e} (significant!) but effect size d={big.mean()/big.std():.3f} (negligible)")
print("Lesson: always report an EFFECT SIZE and a confidence interval alongside the p-value.")
20 tests, all H0 true -> 2 came out significant (expected ~1) -> false positives multiply n=100000, true effect only 0.02: p=3.19e-10 (significant!) but effect size d=0.020 (negligible) Lesson: always report an EFFECT SIZE and a confidence interval alongside the p-value.
Two of the most common ways p-values mislead: testing many things until one "wins" (fix with corrections like Bonferroni or false-discovery-rate control), and confusing a tiny p with a big effect (fix by reporting effect size and a confidence interval). Statistical significance is not practical importance.
A battery is spec’d at 10.0 hours. An engineer tests 50 units (significance-p-values-and-errors--battery_life.xlsx) and suspects they run short. We test H0: mu = 10 against H1: mu < 10, then ask the harder question: how often would a test this size even catch a real shortfall?
try: bat = pd.read_excel("../../data/significance-p-values-and-errors--battery_life.xlsx", sheet_name="Batteries")
except FileNotFoundError: bat = pd.read_excel(BASE+"significance-p-values-and-errors--battery_life.xlsx", sheet_name="Batteries")
# EXPLORE FIRST: size, missing, and a summary by model
print("shape:", bat.shape, "| missing:", bat.isna().sum().sum())
print(bat.groupby("model").life_hours.agg(["size","mean","std"]).round(3))
x=bat["life_hours"]; mu0=10.0; n=len(x)
res=stats.ttest_1samp(x,mu0); p_one=res.pvalue/2
d=(x.mean()-mu0)/x.std(ddof=1)
print(f"n={n}, mean={x.mean():.3f} h, sd={x.std(ddof=1):.3f}")
print(f"t={res.statistic:.2f}, one-sided p={p_one:.4f} -> {'REJECT' if p_one<0.05 else 'fail to reject'} H0 at 5%")
print(f"effect size Cohen\u2019s d = {d:.2f} (small)")
shape: (50, 5) | missing: 0
size mean std
model
A100 20 9.612 0.759
A200 30 9.950 0.764
n=50, mean=9.815 h, sd=0.773
t=-1.70, one-sided p=0.0481 -> REJECT H0 at 5%
effect size Cohen’s d = -0.24 (small)
# Power via statsmodels (no hand-rolled noncentral t): how often would we catch this shortfall?
from statsmodels.stats.power import TTestPower
analysis = TTestPower()
power = analysis.power(effect_size=d, nobs=n, alpha=0.05, alternative="smaller")
print(f"power at the observed effect (n={n}, alpha=.05, one-sided) = {power:.2f}")
print("how many units for more power at this effect size:")
for nn in [50,100,150,200,300]:
pw = analysis.power(effect_size=d, nobs=nn, alpha=0.05, alternative="smaller")
print(f" n={nn:>3}: power = {pw:.2f}")
print("To reliably catch a shortfall this small, the engineer needs a far larger sample.")
power at the observed effect (n=50, alpha=.05, one-sided) = 0.51 how many units for more power at this effect size: n= 50: power = 0.51 n=100: power = 0.77 n=150: power = 0.90 n=200: power = 0.96 n=300: power = 0.99 To reliably catch a shortfall this small, the engineer needs a far larger sample.
The result is just barely significant (one-sided p ≈ 0.048): we reject the 10-hour claim, but only by a hair. And the power is only about 0.51, this test would miss a real shortfall of this size nearly half the time. A borderline p with low power is fragile evidence; the honest move is to report the effect size (d ≈ −0.24), the confidence interval, and to recommend a larger confirmatory sample rather than declaring victory.