Chapter 36: Hypothesis Testing

The Logic of Hypothesis Testing

Every hypothesis test begins with two complementary claims about the world. The null hypothesis $H_0$ represents the status quo or the skeptical position — typically a statement of no effect, no difference, or no association. The alternative hypothesis $H_1$ (or $H_a$) captures the claim we actually want to investigate. The asymmetry is intentional and important: we never prove $H_0$; we can only accumulate evidence that makes it implausible.

The classical setup is:

for a two-sided test, or $H_1: \mu > \mu_0$ / $H_1: \mu < \mu_0$ for one-sided alternatives. The choice of one- vs. two-sided should be made before looking at the data based on the scientific question — not retroactively based on which direction happens to look significant.

Test Statistics and Null Distributions

A test statistic is a single number computed from the sample that summarizes the evidence against $H_0$. Its key property is that its sampling distribution is known (or can be approximated) when $H_0$ is true. For a one-sample t-test:

where $\bar{X}$ is the sample mean, $s$ is the sample standard deviation, and $n$ is the sample size. Under $H_0$, this statistic follows a $t$-distribution with $n-1$ degrees of freedom — exactly when the data are normal, approximately via the Central Limit Theorem for large $n$.

The P-Value: Definition and Correct Interpretation

The p-value is the probability of observing a test statistic at least as extreme as the one computed, assuming $H_0$ is true:

A small p-value means the observed data would be rare under $H_0$ — it is evidence against the null. A large p-value means the data are consistent with the null, not that the null is true.

What a P-Value Is Not

Three misconceptions dominate practice:

• The p-value is not the probability that <!--MATHBLOCK20--> is true. That would require a prior probability on <!--MATHBLOCK21-->, which is the domain of Bayesian inference.
• A p-value below 0.05 does not mean the effect is practically important. With <!--MATHBLOCK22-->, a 0.001-unit difference will be wildly significant.
• A non-significant result does not confirm <!--MATHBLOCK23-->. It means the data were insufficient to detect an effect — which could exist but be small or the study underpowered.

The 0.05 threshold, enshrined by Fisher in the 1920s, is a convention, not a law of nature. Many journals now recommend reporting the exact p-value alongside effect sizes and confidence intervals.

Confidence Intervals as Companions to P-Values

A 95% confidence interval (CI) and a two-sided test at $\alpha = 0.05$ are dual: you reject $H_0: \mu = \mu_0$ at level $\alpha$ if and only if $\mu_0$ falls outside the $(1-\alpha)$ CI. Reporting the CI alongside the p-value communicates the direction and magnitude of the effect, not just its significance.

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
data = rng.normal(loc=5.3, scale=2.0, size=40)  # true mean slightly above 5

# One-sample t-test: H0: mu = 5
t_stat, p_val = stats.ttest_1samp(data, popmean=5.0)
ci = stats.t.interval(0.95, df=len(data)-1,
                       loc=np.mean(data),
                       scale=stats.sem(data))

print(f"t = {t_stat:.3f}, p = {p_val:.4f}")
print(f"95% CI: ({ci[0]:.3f}, {ci[1]:.3f})")
print(f"Sample mean: {np.mean(data):.3f}")

The output reveals whether $\mu_0 = 5$ is plausible given the data, and the CI communicates the range of mean values consistent with the data.

Assumptions and Their Violations

The t-test assumes:

1. Independence of observations.
2. Approximate normality of the sampling distribution of <!--MATHBLOCK30--> (not necessarily of the raw data; CLT helps for <!--MATHBLOCK31-->).
3. For two-sample tests: homogeneity of variance (relaxed by Welch's t-test, which is now the default in most software).

When these assumptions are doubtful, Section 4 covers nonparametric alternatives. Independence is the hardest to repair statistically — it requires redesigning the data collection.

The Two Ways to Be Wrong

Every hypothesis test can make one of two errors, and reducing one tends to inflate the other:

• Type I error (false positive): Reject <!--MATHBLOCK3--> when it is actually true. Its probability is <!--MATHBLOCK4-->, the significance level — the threshold you set before the test.
• Type II error (false negative): Fail to reject <!--MATHBLOCK5--> when it is actually false. Its probability is <!--MATHBLOCK6-->.

The relationship is governed by the effect size, sample size, variability, and $\alpha$. You cannot simultaneously make both $\alpha$ and $\beta$ arbitrarily small without increasing $n$.

Statistical Power

Power is $1 - \beta$: the probability of correctly detecting an effect that truly exists. Power depends on four interrelated quantities:

1. Effect size <!--MATHBLOCK12--> — how large the true departure from <!--MATHBLOCK13--> is (often standardized as Cohen's <!--MATHBLOCK14-->).
2. Sample size <!--MATHBLOCK15-->.
3. Significance level <!--MATHBLOCK16-->.
4. Population variability <!--MATHBLOCK17-->.

For a one-sample z-test against $H_0: \mu = \mu_0$ with known $\sigma$, the power of a two-sided test at level $\alpha$ is:

where $\Phi$ is the standard normal CDF. For practical purposes (moderate $n$), the second term is negligible.

Why Power Matters

An underpowered study (say, power = 0.30) wastes resources: when it does find significance, the effect estimate is likely to be inflated (the "winner's curse"), and when it does not, the null result is uninformative. Conventional minimum power is 0.80, meaning at most a 20% chance of missing a real effect. Many professional guidelines (e.g., FDA, APA) require a prospective power analysis.

Computing Sample Size

Rearranging the power formula gives the required $n$ for a desired power $1-\beta$:

For $\alpha = 0.05$, power = 0.80, $z_{\alpha/2} = 1.96$, $z_\beta = 0.84$:

Substituting Cohen's $d = \delta/\sigma$ gives $n \approx (2.80/d)^2$. For a medium effect ($d = 0.5$): $n \approx 31$ per group.

from statsmodels.stats.power import TTestIndPower
import numpy as np
import matplotlib.pyplot as plt

analysis = TTestIndPower()

# Required sample size for various effect sizes
for d in [0.2, 0.5, 0.8]:
    n = analysis.solve_power(effect_size=d, power=0.80,
                              alpha=0.05, alternative='two-sided')
    print(f"Cohen's d={d:.1f} -> n={n:.0f} per group")

# Power curve: power vs sample size for d=0.5
n_range = np.arange(10, 200, 5)
powers = analysis.power(effect_size=0.5, nobs1=n_range,
                         alpha=0.05, alternative='two-sided')

plt.figure(figsize=(7, 4))
plt.plot(n_range, powers, 'b-o', markersize=3)
plt.axhline(0.80, color='red', linestyle='--', label='Power = 0.80')
plt.xlabel('Sample size per group'); plt.ylabel('Power')
plt.title("Power curve: two-sample t-test, d=0.5, α=0.05")
plt.legend(); plt.tight_layout(); plt.show()

Practical Guidelines for Power Analysis

• Estimate effect size from prior literature, pilot data, or the minimum practically meaningful difference — not from the current sample (that inflates estimates).
• Report power analyses in papers as a justification for <!--MATHBLOCK32-->, not as post-hoc rationalization.
• Post-hoc power (computed after a non-significant result using the observed effect size) is widely criticized as circular and uninformative — prefer confidence intervals instead.
• One-sided tests have higher power but require scientific justification for the directional prediction made before data collection.

Common Pitfalls

• Using <!--MATHBLOCK33--> for thousands of tests without correction inflates the familywise error rate dramatically (see Section 5).
• "Peeking" at data and stopping early when <!--MATHBLOCK34--> invalidates the Type I error rate unless a sequential testing procedure is used.
• Increasing <!--MATHBLOCK35--> indefinitely: with enough subjects, even trivially small effects become significant — always pair p-values with effect sizes.

Choosing the Right Parametric Test

Parametric tests assume a specific distributional form (usually normality) and achieve higher power than nonparametric alternatives when assumptions hold. The choice of test depends on:

• Number of groups: one, two, or more than two.
• Whether samples are independent or paired/repeated.
• Data type: continuous (t-tests, ANOVA) vs. categorical counts (chi-square).

T-Tests

One-Sample T-Test

Tests whether the population mean equals a specified value $\mu_0$. Statistic: $t = (\bar{X} - \mu_0)/(s/\sqrt{n})$, df = $n-1$.

Two-Sample (Independent) T-Test

Compares means of two independent groups. Welch's version (unequal variances) is almost always preferred over Student's version:

with degrees of freedom approximated by the Welch-Satterthwaite equation. In scipy.stats.ttest<em>ind, pass equal</em>var=False (the default in R's t.test).

Paired T-Test

When observations are matched (before/after, same subject, matched pairs), compute $d_i = X_{i1} - X_{i2}$ and run a one-sample t-test on $d_i$. Pairing removes between-subject variability and increases power substantially.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
# Two independent groups
group_a = rng.normal(50, 10, 35)
group_b = rng.normal(55, 12, 35)
t, p = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"Welch's t-test: t={t:.3f}, p={p:.4f}")

# Paired: same subjects before/after treatment
before = rng.normal(70, 8, 25)
after = before + rng.normal(3, 4, 25)  # treatment adds ~3 units
t_p, p_p = stats.ttest_rel(before, after)
print(f"Paired t-test: t={t_p:.3f}, p={p_p:.4f}")

# Cohen's d for paired test
d = np.mean(before - after) / np.std(before - after, ddof=1)
print(f"Paired Cohen's d = {d:.3f}")

Chi-Square Tests

Goodness-of-Fit

Tests whether observed categorical frequencies match expected proportions under $H_0$. For $k$ categories:

Distributed as $\chi^2_{k-1}$ under $H_0$ when all $E_i \geq 5$ (rule of thumb).

Test of Independence

For a contingency table with $r$ rows and $c$ columns, tests whether the row and column variables are independent. Degrees of freedom: $(r-1)(c-1)$. Note: association does not imply causation, and chi-square gives no information about the direction or magnitude of an association — report Cramér's V as the effect size.

from scipy.stats import chi2_contingency, chisquare
import numpy as np

# Goodness-of-fit: is a die fair?
observed = np.array([18, 22, 19, 17, 24, 20])
expected = np.full(6, np.sum(observed) / 6)
chi2, p = chisquare(observed, f_exp=expected)
print(f"Die fairness: chi2={chi2:.3f}, p={p:.4f}")

# Independence test: gender vs. preference
table = np.array([[45, 55],   # male: product A, B
                  [60, 40]])  # female: product A, B
chi2, p, dof, expected_counts = chi2_contingency(table)
cramers_v = np.sqrt(chi2 / (table.sum() * (min(table.shape) - 1)))
print(f"Independence: chi2={chi2:.3f}, p={p:.4f}, df={dof}")
print(f"Cramér's V = {cramers_v:.3f}")

One-Way ANOVA

ANOVA (Analysis of Variance) tests whether the means of three or more independent groups are all equal:

The F-statistic compares between-group variance to within-group variance:

A significant F only tells you some means differ — follow up with post-hoc tests (Tukey HSD, Bonferroni correction) to locate which pairs differ.

Two-Way ANOVA

Extends to two categorical factors and their interaction. An interaction means the effect of one factor depends on the level of the other — always inspect the interaction term before interpreting main effects.

import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
import pandas as pd

rng = np.random.default_rng(1)
# Three fertilizer groups
groups = [rng.normal(30, 5, 20),
          rng.normal(35, 5, 20),
          rng.normal(33, 5, 20)]

F, p = stats.f_oneway(*groups)
print(f"One-way ANOVA: F={F:.3f}, p={p:.4f}")

# Post-hoc Tukey HSD via statsmodels
df = pd.DataFrame({'yield': np.concatenate(groups),
                   'group': np.repeat(['A','B','C'], 20)})
model = ols('yield ~ C(group)', data=df).fit()
aov = sm.stats.anova_lm(model, typ=1)
print(aov)

from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(df['yield'], df['group'], alpha=0.05)
print(tukey)

Assumptions and Diagnostics

ANOVA assumes independence, normality within groups, and homogeneity of variances (Levene's test or Bartlett's test). In practice, ANOVA is fairly robust to moderate violations of normality when group sizes are equal. Unequal group sizes combined with unequal variances, however, can seriously inflate Type I error.

When Parametric Assumptions Fail

Parametric tests extract power by exploiting assumed distributional structure. When that structure is absent or uncertain — with small samples, heavy-tailed distributions, ordinal data, or clear outliers — nonparametric tests provide valid inference without assuming normality. The tradeoff: they generally have lower power when the parametric assumptions actually hold (typically ~5-15% efficiency loss), but they can have higher power when those assumptions are badly violated.

Nonparametric tests typically work on ranks of the data rather than the raw values, making them resistant to outliers and appropriate for ordinal measurements.

Mann-Whitney U (Wilcoxon Rank-Sum) Test

The Mann-Whitney U test is the nonparametric alternative to the two-sample t-test. Under $H_0$, it tests whether the two samples come from the same distribution (equivalently, whether one distribution is stochastically greater than the other). The test statistic $U$ counts the number of times an observation from group 1 exceeds an observation from group 2 across all pairs.

For large samples, $U$ is approximately normal:

The effect size is the rank-biserial correlation $r = 1 - 2U/(n_1 n_2)$, ranging from -1 to +1.

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
# Skewed data: log-normal distributions
group1 = rng.lognormal(mean=0.0, sigma=1.0, size=30)
group2 = rng.lognormal(mean=0.5, sigma=1.2, size=30)

# Mann-Whitney U
u_stat, p_mw = stats.mannwhitneyu(group1, group2, alternative='two-sided')
print(f"Mann-Whitney U={u_stat:.1f}, p={p_mw:.4f}")

# Compare with t-test (inappropriate here)
t, p_t = stats.ttest_ind(group1, group2, equal_var=False)
print(f"Welch's t={t:.3f}, p={p_t:.4f}  (less reliable for skewed data)")

# Rank-biserial correlation (effect size)
r_rb = 1 - 2*u_stat / (len(group1)*len(group2))
print(f"Rank-biserial r = {r_rb:.3f}")

Wilcoxon Signed-Rank Test

The nonparametric counterpart to the paired t-test. It ranks the absolute differences $|d_i| = |X_{i1} - X_{i2}|$, then computes the sum of ranks for positive and negative differences separately. Under $H_0$ (symmetric distribution of $d_i$ around 0), the positive and negative rank sums should be equal.

before = np.array([85, 90, 78, 92, 88, 76, 95, 83, 79, 91])
after  = np.array([88, 93, 79, 94, 92, 80, 96, 87, 82, 95])

# Wilcoxon signed-rank test
w_stat, p_wsr = stats.wilcoxon(before, after, alternative='two-sided')
print(f"Wilcoxon signed-rank: W={w_stat:.1f}, p={p_wsr:.4f}")

Kruskal-Wallis Test

The nonparametric analogue of one-way ANOVA. All $N$ observations are ranked together; the test statistic measures whether the mean ranks differ across $k$ groups:

Under $H_0$, $H \sim \chi^2_{k-1}$ (approximately). A significant result requires post-hoc pairwise comparisons (e.g., Dunn's test with appropriate correction).

from scipy.stats import kruskal

rng = np.random.default_rng(5)
g1 = rng.exponential(scale=2.0, size=25)
g2 = rng.exponential(scale=3.0, size=25)
g3 = rng.exponential(scale=2.5, size=25)

H, p_kw = kruskal(g1, g2, g3)
print(f"Kruskal-Wallis: H={H:.3f}, p={p_kw:.4f}")

# Post-hoc: pairwise Mann-Whitney with Bonferroni correction
pairs = [(g1,g2,'G1 vs G2'), (g1,g3,'G1 vs G3'), (g2,g3,'G2 vs G3')]
for a, b, label in pairs:
    _, p = stats.mannwhitneyu(a, b, alternative='two-sided')
    print(f"  {label}: p={p:.4f} (Bonferroni adj: p={min(p*3,1):.4f})")

Kolmogorov-Smirnov and Anderson-Darling Tests

Sometimes the question is not about central tendency but about the entire distribution. The two-sample Kolmogorov-Smirnov test tests whether two samples come from the same distribution by comparing their empirical CDFs:

The Anderson-Darling test places more weight on the tails and is generally more powerful for detecting tail differences.

# KS test for distributional equality
g_normal = rng.normal(0, 1, 100)
g_heavy  = rng.standard_t(df=3, size=100)  # heavier tails

ks_stat, p_ks = stats.ks_2samp(g_normal, g_heavy)
print(f"KS test: D={ks_stat:.3f}, p={p_ks:.4f}")

Practical Decision Guide

• Two independent groups, non-normal or small sample: Mann-Whitney U.
• Paired/matched samples, non-normal: Wilcoxon signed-rank.
• Three or more groups, non-normal: Kruskal-Wallis + Dunn post-hoc.
• Testing full distributional equality: KS or Anderson-Darling.
• Ordinal data of any kind: prefer nonparametric tests by default.

Nonparametric does not mean assumption-free. Mann-Whitney U, for instance, assumes the two distributions are identical under $H_0$ (same shape, just potentially shifted). If the shapes differ, the test detects distributional differences more broadly, which may or may not be what you want.

The Multiple Comparisons Problem

Every individual hypothesis test has a Type I error rate $\alpha$. But when you run $m$ independent tests simultaneously, each at level $\alpha = 0.05$, the probability of making at least one false rejection grows rapidly:

For $m = 20$: $1 - 0.95^{20} \approx 0.64$. For $m = 100$: $\approx 0.99$. This is the familywise error rate (FWER) problem. In genomics (testing 20,000 genes), untreated multiple comparisons make virtually every study a false-discovery machine.

Bonferroni Correction

The simplest FWER control: test each hypothesis at $\alpha^* = \alpha / m$. For $m = 20$ tests and $\alpha = 0.05$, use the threshold $0.05/20 = 0.0025$. Equivalently, multiply each p-value by $m$ and use the original $\alpha$ threshold (capping adjusted p-values at 1).

Bonferroni is exact (not just approximate) when tests are independent, and conservative (guards against more than necessary) when they are positively correlated. It is appropriate when any false positive is catastrophic — but its conservatism costs power and is often excessive in exploratory analysis.

Holm-Bonferroni (Step-Down) Correction

Holm's method is uniformly more powerful than Bonferroni while still controlling the FWER strongly. Algorithm:

1. Sort p-values: <!--MATHBLOCK15-->.
2. For the <!--MATHBLOCK16-->-th smallest, compare against <!--MATHBLOCK17-->.
3. Reject <!--MATHBLOCK18--> only if all <!--MATHBLOCK19--> with <!--MATHBLOCK20--> were also rejected.

Because the threshold increases as we move to larger p-values, Holm rejects more hypotheses than Bonferroni.

Benjamini-Hochberg (BH) False Discovery Rate

When many hypotheses are tested in exploratory research (gene expression, feature selection), controlling FWER is often too stringent. The False Discovery Rate (FDR) controls the expected proportion of false discoveries among all rejections:

where $V$ is the number of false positives and $R$ is the total number of rejections (with the convention $0/0 = 0$).

The Benjamini-Hochberg procedure at level $q$:

1. Sort p-values: <!--MATHBLOCK25-->.
2. Find the largest <!--MATHBLOCK26--> such that <!--MATHBLOCK27-->.
3. Reject all <!--MATHBLOCK28--> for <!--MATHBLOCK29-->.

BH controls FDR at level $q$ exactly when tests are independent, and at $q \cdot m_0/m$ when they are positively correlated (Benjamini-Yekutieli for arbitrary dependence). In practice, BH is the standard correction in genomics, neuroimaging, and any setting where hundreds of tests are expected to include many true positives.

import numpy as np
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng(42)

# Simulate 100 tests: 90 nulls, 10 true effects
m = 100
null_p  = rng.uniform(0, 1, 90)          # p-values under H0: uniform
effect_p = rng.beta(0.5, 10, 10)         # p-values with true effect: small
p_values = np.concatenate([null_p, effect_p])
true_null = np.array([True]*90 + [False]*10)

methods = ['bonferroni', 'holm', 'fdr_bh']
for method in methods:
    reject, p_adj, _, _ = multipletests(p_values, alpha=0.05, method=method)
    tp = np.sum(reject & ~true_null)   # true positives
    fp = np.sum(reject & true_null)    # false positives
    print(f"{method:12s}: reject={reject.sum():3d}, "
          f"TP={tp}, FP={fp}, "
          f"FDR={fp/max(reject.sum(),1):.2f}")

When to Use Each Method

• Bonferroni: Confirmatory studies, regulatory contexts, small number of tests (<!--MATHBLOCK32-->), any false positive is costly.
• Holm-Bonferroni: Same contexts as Bonferroni but uniformly preferred — strictly more powerful with the same FWER guarantee.
• Benjamini-Hochberg: Exploratory studies, large <!--MATHBLOCK33--> (<!--MATHBLOCK34-->), some false positives are acceptable, downstream validation is planned.
• Benjamini-Yekutieli: When test statistics are negatively correlated or dependence structure is unknown and you need FDR control without assuming independence.

Practical Considerations

Defining the family of hypotheses is the most subjective step. Should you correct across all tests in a paper, or only within each experiment? There is no universal answer, but the correction scope should be stated explicitly. Common guidelines:

• Correct within a single coherent set of related tests (e.g., all post-hoc pairwise comparisons from one ANOVA).
• Do not combine tests from logically unrelated experiments into one family.
• Always report the uncorrected p-values alongside corrected ones so readers can evaluate the evidence directly.

import matplotlib.pyplot as plt

# Visualize BH procedure
idx = np.argsort(p_values)
sorted_p = p_values[idx]
m_tests = len(sorted_p)
bh_thresholds = (np.arange(1, m_tests+1) / m_tests) * 0.05

plt.figure(figsize=(8, 5))
plt.scatter(range(1, m_tests+1), sorted_p,
            c=['red' if not true_null[i] else 'blue' for i in idx],
            alpha=0.6, s=20, label='p-values (red=true effect)')
plt.plot(range(1, m_tests+1), bh_thresholds, 'k-', lw=2,
         label='BH threshold line')
plt.xlabel('Rank'); plt.ylabel('p-value')
plt.title('Benjamini-Hochberg procedure visualization')
plt.legend(); plt.tight_layout(); plt.show()