07 - Hypothesis Testing: t-test, z-test, p-value

A machine is supposed to fill 500 g packs. You draw 25 from the line, weigh them and find a mean of 497.6 g with a sample standard deviation of 6 g: is the machine out of adjustment, or is it just the normal fluctuation of the sample? With the numbers of the example the statistic is exactly −2.00, but the conclusion depends on which distribution you feed it to: with the Normal the p-value is 0.0455 and you reject the hypothesis that the machine is fine, with Student’s t it is 0.0569 and you do not reject it.

The same data, two opposite verdicts at the 5% level. This article explains why, and takes you from “I have a formula” to “I know what it guarantees”. It is the seventh in the Probability and Statistics series; you need the estimators from the previous article and the Normal from the central limit theorem.

How a hypothesis test reasons

A test starts from two statements about the parameter you do not know. The null hypothesis H₀ is the one to be discarded only in the face of strong evidence (in the example, “the mean of the packs is 500 g”); the alternative H₁ is what you accept if you discard the null (“the mean is different from 500 g”).

Then you choose a test statistic, a function of the data whose distribution, if H₀ is true, is known: in the example it is T = (X̄ − 500)/(s/√n), that is, the gap between the observed mean and the stated one measured in standard errors (here 6/√25 = 1.2 g, hence (497.6 − 500)/1.2 = −2). Finally you fix a rule: you reject H₀ if the statistic falls in a zone so far out that it is unlikely under H₀.

Such a rule can go wrong in two ways. The type I error is rejecting H₀ when it is true (a false alarm); its probability, which you choose before looking at the data, is called the level α (usually 0.05). The type II error is not rejecting H₀ when it is false (a missed opportunity); the probability of avoiding it is the power, and it depends on how large the true effect is, on how much data you have and on α. The definitions of error, level and power are those of the parametric hypothesis testing lecture of the MIT 18.650 course.

The p-value is the smallest level α at which your data would have led to rejecting H₀ (same source): it is equivalent to the probability, computed assuming H₀ is true, of obtaining a statistic at least as extreme as the observed one. In the example, with the Normal, the probability of falling beyond ±2 is 0.0455. The p-value is a random variable, because it depends on the sample. The link with confidence intervals therefore remains: the two-sided test at level 5% rejects H₀ exactly when the stated value (500 g) is outside the 95% confidence interval, as the next program shows.

z or t: it depends on what you know about σ

The z-test applies when the standard deviation σ of the population is known: the statistic Z = (X̄ − μ₀)/(σ/√n) is standard Normal under H₀ if the data are Normal.

When σ is unknown you replace it with the sample standard deviation S, and the statistic becomes the ratio between a Normal and a random quantity; it has heavier tails, and it is a Student’s t with n − 1 degrees of freedom (the number of independent deviations that enter S). For small n the difference is large, for large n it disappears. The program reads the data of the example both ways: the same statistic −2, two different tails.

# The same 25 packs read with a z-test (sigma "known" = 6) and with a t-test (estimated sigma)
import numpy as np
from scipy import stats

n, media, s, mu0 = 25, 497.6, 6.0, 500.0
se = s / np.sqrt(n)                       # standard error = 6/5 = 1.2 g
stat = (media - mu0) / se                 # same numerical statistic in both tests

p_z = 2 * stats.norm.sf(abs(stat))        # two-sided tail of the standard Normal
p_t = 2 * stats.t.sf(abs(stat), df=n - 1) # two-sided tail of the t with 24 degrees of freedom
print(f"statistic = {stat:.3f}")
print(f"z-test: p = {p_z:.4f}  -> at 5% {'rejects' if p_z < 0.05 else 'does not reject'} H0")
print(f"t-test: p = {p_t:.4f}  -> at 5% {'rejects' if p_t < 0.05 else 'does not reject'} H0")

# 95% confidence interval built with the t: it contains 500 iff the t-test does not reject
q = stats.t.ppf(0.975, df=n - 1)
print(f"t(24, 0.975) quantile = {q:.4f} versus z(0.975) = {stats.norm.ppf(0.975):.4f}")
print(f"95% CI (t): [{media - q * se:.2f}, {media + q * se:.2f}]")

# Output:
# statistic = -2.000
# z-test: p = 0.0455  -> at 5% rejects H0
# t-test: p = 0.0569  -> at 5% does not reject H0
# t(24, 0.975) quantile = 2.0639 versus z(0.975) = 1.9600
# 95% CI (t): [495.12, 500.08]

The statistic is identical, −2.000, but the tail of the t with 24 degrees of freedom is 0.0569 against 0.0455 for the Normal, because the 97.5% quantile is 2.0639 and not 1.96. The 95% confidence interval built with the t is [495.12; 500.08] and contains 500: consistent with the t-test not rejecting. In this example σ is not known (6 g is the estimate from the sample), so the correct test is the t-test; the “machine out of adjustment” verdict of the z-test is a false alarm that comes from having used the wrong distribution. This is not an observation about one case: the next program measures how often it happens.

To measure it, I simulate one million samples from a Normal with mean 0 (H₀ true) for each n, and count how often the t-test rejects with its threshold (quantile of the t) and the “z-test with S” with the threshold 1.96. Under H₀ the statistic is exactly a t with n − 1 degrees of freedom, so the exact frequency of the second is 2 · P(t_n−1 > 1.96), which is computed with the survival function of the t.

# False positive rate under H0: z-test with estimated sigma versus t-test, for small samples
import numpy as np
from scipy import stats

R, BLOCCO, z = 1_000_000, 100_000, stats.norm.ppf(0.975)
rng = np.random.default_rng(99)
print("  n  t-test (sim.)  z with S (sim.)  z with S (exact)")
for n in (5, 10, 30, 100):
    crit_t = stats.t.ppf(0.975, df=n - 1)              # correct threshold: quantile of the t
    rif_t = rif_z = 0
    for _ in range(R // BLOCCO):
        X = rng.normal(0.0, 2.0, size=(BLOCCO, n))     # H0 true: mean 0 (sigma = 2, but the test does not use it)
        stat = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
        rif_t += (np.abs(stat) > crit_t).sum()
        rif_z += (np.abs(stat) > z).sum()              # wrong threshold: quantile of the Normal
    esatto = 2 * stats.t.sf(z, df=n - 1)               # under H0 the statistic is a t(n-1)
    print(f"{n:3d}  {rif_t / R:12.4f}  {rif_z / R:15.4f}  {esatto:17.4f}")
print(f"standard error of each simulated value, with R = {R}: about {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
#   n  t-test (sim.)  z with S (sim.)  z with S (exact)
#   5        0.0503           0.1222             0.1216
#  10        0.0501           0.0816             0.0816
#  30        0.0503           0.0599             0.0597
# 100        0.0501           0.0529             0.0528
# standard error of each simulated value, with R = 1000000: about 0.0002

The t-test rejects 5.0% of the time for every n, as promised. The “z-test with S” rejects 12.2% of the time with 5 observations (exact 12.16%), 8.2% with 10, 6.0% with 30 and 5.3% with 100: the simulation coincides with the exact calculation within about 0.0006. About a hundred observations are needed to make the error of using the Normal in place of the t negligible, and this is why the rule “with more than 30 data points you can use the Normal” is an acceptable shortcut only if a true level of about 6% in place of 5% is good enough for you.

Two samples: Student or Welch

To compare the means of two independent groups, for example the response time of two versions of a page, there are two t-tests. Student’s assumes that the two groups have the same variance, estimates it with a common value (pooled) and uses n₁ + n₂ − 2 degrees of freedom. Welch’s estimates the two variances separately, uses √(s₁²/n₁ + s₂²/n₂) as the standard error and the approximate Welch-Satterthwaite degrees of freedom. In SciPy scipy.stats.ttest_ind has equal_var=True as its default value, that is, it runs Student’s: the documentation of the function says so. To get Welch’s you must write equal_var=False.

The program writes the two tests in vectorized form, checks that the p-values coincide with those of SciPy on one sample, then simulates 200,000 pairs of samples with the same mean (H₀ true) but different variances (standard deviations 3 and 1), in three arrangements: small group with the larger variance, small group with the smaller variance, equal sizes. For this comparison there is no simple closed formula, so the check is the simulation alone, with its standard error.

# Two samples, equal means but different variances: Student's t (equal variances) versus Welch
import numpy as np
from scipy import stats

R = 200_000
rng = np.random.default_rng(8)

def p_values(a, b):
    """Two-sided p-values of Student (common variance) and Welch; each row is a sample."""
    n1, n2 = a.shape[1], b.shape[1]
    m1, m2 = a.mean(axis=1), b.mean(axis=1)
    v1, v2 = a.var(axis=1, ddof=1), b.var(axis=1, ddof=1)
    # Student: estimated common variance (pooled), n1 + n2 - 2 degrees of freedom
    sp2 = ((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2)
    t_st = (m1 - m2) / np.sqrt(sp2 * (1 / n1 + 1 / n2))
    p_st = 2 * stats.t.sf(np.abs(t_st), df=n1 + n2 - 2)
    # Welch: separate standard errors, Welch-Satterthwaite degrees of freedom
    e1, e2 = v1 / n1, v2 / n2
    t_w = (m1 - m2) / np.sqrt(e1 + e2)
    df_w = (e1 + e2) ** 2 / (e1 ** 2 / (n1 - 1) + e2 ** 2 / (n2 - 1))
    p_w = 2 * stats.t.sf(np.abs(t_w), df=df_w)
    return p_st, p_w

# check: the formulas agree with scipy on one sample
a0, b0 = rng.normal(0, 3, (1, 10)), rng.normal(0, 1, (1, 40))
mio = p_values(a0, b0)
sci = (stats.ttest_ind(a0[0], b0[0], equal_var=True).pvalue, stats.ttest_ind(a0[0], b0[0], equal_var=False).pvalue)
print(f"p Student: mine {mio[0][0]:.6f}, scipy {sci[0]:.6f}; p Welch: mine {mio[1][0]:.6f}, scipy {sci[1]:.6f}")
print("case                                      Student   Welch")
casi = {"small sample has the larger variance": (10, 3.0, 40, 1.0),
        "small sample has the smaller variance": (10, 1.0, 40, 3.0),
        "equal sizes, different variances": (25, 3.0, 25, 1.0)}
for nome, (n1, s1, n2, s2) in casi.items():
    p_st, p_w = p_values(rng.normal(0, s1, (R, n1)), rng.normal(0, s2, (R, n2)))   # same mean: H0 true
    print(f"{nome:42s} {(p_st < 0.05).mean():6.4f}  {(p_w < 0.05).mean():6.4f}")
print(f"standard error with R = {R}: about {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# p Student: mine 0.040557, scipy 0.040557; p Welch: mine 0.268662, scipy 0.268662
# case                                      Student   Welch
# small sample has the larger variance       0.2568  0.0501
# small sample has the smaller variance      0.0014  0.0496
# equal sizes, different variances           0.0528  0.0493
# standard error with R = 200000: about 0.0005

The first two lines are the important verdicts. If the small group (10 observations) is also the more dispersed one, Student rejects a true H₀ 25.7% of the time instead of 5%, because the estimated common variance is dominated by the large group and underestimates the uncertainty of the small one. In the opposite case it almost never rejects (0.14%), with a loss of power that cannot be seen from the level.

Welch stays around 5% in all three cases (0.0493-0.0501), and even with equal sizes Student holds (0.0528). Practical conclusion: Welch costs almost nothing when the variances are equal and protects when they are not, so with SciPy it pays to write equal_var=False out of habit.

How much data you need: power

Power answers a question to ask before collecting the data: if the true difference is the one I care about, with how much data will I discover it with high probability? The effect is measured in units of standard deviation, δ = (μ₁ − μ₀)/σ. Under the alternative the statistic of the t-test follows a noncentral t with noncentrality parameter δ√n; scipy.stats.nct computes its probabilities. The formula found in textbooks, n ≥ ((z_1−α/2 + z_0.80)/δ)², uses the Normal and therefore underestimates n slightly. The program compares the formula with the exact power and confirms it with a simulation.

# Power of the one-sample t-test (two-sided, alpha = 0.05) and sample size
import numpy as np
from scipy import stats

alfa, delta = 0.05, 0.5            # delta = (mu1 - mu0)/sigma, "medium" effect according to Cohen
def potenza(n, d):
    """Exact: under H1 the statistic is a noncentral t with parameter d*sqrt(n)."""
    crit = stats.t.ppf(1 - alfa / 2, df=n - 1)
    nc = d * np.sqrt(n)
    return stats.nct.sf(crit, n - 1, nc) + stats.nct.cdf(-crit, n - 1, nc)

z = stats.norm.ppf
n_formula = ((z(1 - alfa / 2) + z(0.80)) / delta) ** 2     # formula with the Normal
print(f"Normal formula: n >= {n_formula:.2f} -> {int(np.ceil(n_formula))}")
for n in (32, 33, 34):
    print(f"n = {n}: exact power of the t-test = {potenza(n, delta):.4f}")
n_t = next(n for n in range(2, 1000) if potenza(n, delta) >= 0.80)
print(f"first n with power >= 0.80 (t-test): {n_t}")

# Check by simulation with n = 34
R, n = 200_000, 34
rng = np.random.default_rng(2026)
X = rng.normal(delta, 1.0, size=(R, n))
t = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
sim = (np.abs(t) > stats.t.ppf(1 - alfa / 2, df=n - 1)).mean()
print(f"n = {n}: simulated power = {sim:.4f} +/- {np.sqrt(sim * (1 - sim) / R):.4f}")
for d in (0.4, 0.2):
    primo = next(n for n in range(2, 2000) if potenza(n, d) >= 0.80)
    print(f"delta = {d}: first n with power >= 0.80 = {primo} (Normal formula: "
          f"{int(np.ceil(((z(1 - alfa / 2) + z(0.80)) / d) ** 2))})")
print(f"machine: with n = 25 and delta = 0.4 the exact power is {potenza(25, 0.4):.4f}")

# Output:
# Normal formula: n >= 31.40 -> 32
# n = 32: exact power of the t-test = 0.7823
# n = 33: exact power of the t-test = 0.7954
# n = 34: exact power of the t-test = 0.8078
# first n with power >= 0.80 (t-test): 34
# n = 34: simulated power = 0.8086 +/- 0.0009
# delta = 0.4: first n with power >= 0.80 = 52 (Normal formula: 50)
# delta = 0.2: first n with power >= 0.80 = 199 (Normal formula: 197)
# machine: with n = 25 and delta = 0.4 the exact power is 0.4840

For a medium effect (δ = 0.5) the formula says n ≥ 32, but with 32 observations the exact power of the t-test is 0.7823, with 33 it is 0.7954 and only with 34 does it reach 0.8078. The simulation with 34 observations gives 0.8086 with a standard error of 0.0009.

For the machine exercise, if you want to notice a gap of 2.4 g with σ = 6 g (δ = 0.4) you need 52 packs, not 25: with 25 the exact power is 0.4840, that is, less than one time in two, and a “non-significant” result like the one in the example does not prove that the machine is fine. With δ = 0.2 you need 199 observations: halving the effect to be detected requires about four times the data.

What a p-value says, and what it does not

The statement of the American Statistical Association of 7 March 2016 lists six principles on the use of the p-value, reported in the ASA press release; the full statement, by Wasserstein and Lazar, is in The American Statistician, 70(2), with DOI 10.1080/00031305.2016.1154108. Two principles concern the most widespread misreadings: a p-value indicates how incompatible the data are with a specified model, but it does not measure the probability that the hypothesis under study is true, nor the probability that the data were produced by chance alone; and it does not measure the size of an effect or its importance. A third says that conclusions should not be based only on whether a threshold is crossed.

The first program measures three facts about the p-value, with data that contain no effect at all. First, that under H₀ it is uniform on [0, 1], that is, that “p < 0.05” happens in 5% of cases, and also “p < 0.50” in 50%. Second, that 20 independent tests give a probability 1 − 0.95²⁰ = 0.6415 of at least one p < 0.05, while the Bonferroni threshold 0.05/20 brings that probability back under 5%. Third, that looking at the data every 10 observations up to 100 and stopping at the first p < 0.05 inflates the false positives; here there is no closed formula, the simulation is what holds.

# Under H0 the p-value is uniform on [0, 1]; what happens with many tests or peeking at the data along the way
import numpy as np
from scipy import stats

def pvalue_t(X):
    """Two-sided p-value of the one-sample t-test (H0: mean 0) on the rows of X."""
    n = X.shape[-1]
    t = np.sqrt(n) * X.mean(axis=-1) / X.std(axis=-1, ddof=1)
    return 2 * stats.t.sf(np.abs(t), df=n - 1)

rng = np.random.default_rng(20)

# 1) uniformity of the p-value when H0 is true
R, n = 100_000, 30
p = pvalue_t(rng.normal(0, 1, size=(R, n)))
ks = stats.kstest(p, "uniform")
print(f"1) p < 0.05: {(p < 0.05).mean():.4f}; p < 0.50: {(p < 0.50).mean():.4f}; KS D = {ks.statistic:.4f}")

# 2) twenty independent tests, no true effect: probability of at least one p < threshold
R, m = 20_000, 20
P = pvalue_t(rng.normal(0, 1, size=(R, m, n)))            # shape (R, m)
for nome, soglia in (("threshold 0.05   ", 0.05), ("Bonferroni 0.05/20", 0.05 / m)):
    sim = (P < soglia).any(axis=1).mean()
    print(f"2) {nome}: at least one false positive, simulated {sim:.4f}, exact {1 - (1 - soglia) ** m:.4f}")

# 3) looking at the data every 10 observations, up to 100, and stopping at the first p < 0.05
R = 100_000
X = rng.normal(0, 1, size=(R, 100))
S1, S2 = X.cumsum(axis=1), (X ** 2).cumsum(axis=1)
fermati = np.zeros(R, dtype=bool)
for k in range(10, 101, 10):
    media = S1[:, k - 1] / k
    var = (S2[:, k - 1] - k * media ** 2) / (k - 1)
    pk = 2 * stats.t.sf(np.abs(np.sqrt(k) * media / np.sqrt(var)), df=k - 1)
    fermati |= pk < 0.05
print(f"3) 10 interim checks: false positives {fermati.mean():.4f} (nominal 0.05); only at the end: "
      f"{(pk < 0.05).mean():.4f}")

# Output:
# 1) p < 0.05: 0.0491; p < 0.50: 0.4995; KS D = 0.0015
# 2) threshold 0.05   : at least one false positive, simulated 0.6459, exact 0.6415
# 2) Bonferroni 0.05/20: at least one false positive, simulated 0.0500, exact 0.0488
# 3) 10 interim checks: false positives 0.1987 (nominal 0.05); only at the end: 0.0504

The three results are consistent with the theory. The p-value has frequency 0.0491 below 0.05 and 0.4995 below 0.50, and the Kolmogorov-Smirnov statistic against the uniform is 0.0015, so there is no observable departure. With twenty tests the probability of at least one false positive is 0.6459 simulated against 0.6415 exact (standard error 0.0034 with 20,000 simulations), and 0.0500 against 0.0488 with Bonferroni.

With ten interim checks the false positives rise to 0.1987, almost four times the stated level, while the single test at the end stays at 5.04%. A test at level 5% is valid only if the plan (how many tests, when to stop) is decided before looking at the data.

The most damaging misunderstanding remains: thinking that “p < 0.05” means “H₀ has a 5% probability of being true”. The probability that H₀ is true given a significant result cannot be computed from the p-value: it also depends on the fraction of true hypotheses you are putting to the test and on the power.

The program builds a world in which a fraction π₀ of the hypotheses are null and the others have an effect δ = 0.5, with 34 observations per test (power 0.8078), and counts among the significant results how many are true nulls. The exact formula is π₀α/(π₀α + (1 − π₀)·power), and it is obtained with Bayes’ rule.

# Among the "significant" results, how many are true null hypotheses? It depends on how many true hypotheses are put to the test.
import numpy as np
from scipy import stats

alfa, n, delta = 0.05, 34, 0.5
crit = stats.t.ppf(1 - alfa / 2, df=n - 1)
nc = delta * np.sqrt(n)
potenza = stats.nct.sf(crit, n - 1, nc) + stats.nct.cdf(-crit, n - 1, nc)
print(f"alpha = {alfa}, power (n = {n}, effect {delta}) = {potenza:.4f}")

rng = np.random.default_rng(5)
R, BLOCCO = 2_000_000, 100_000
print("share of true null hypotheses   exact   simulated  (among results with p < 0.05)")
for pi0 in (0.5, 0.9, 0.99):
    significativi = nulle_fra_significativi = 0
    for _ in range(R // BLOCCO):
        vera_nulla = rng.random(BLOCCO) < pi0              # H0 true with probability pi0
        media = np.where(vera_nulla, 0.0, delta)           # effect 0 or delta
        X = rng.normal(media[:, None], 1.0, size=(BLOCCO, n))
        t = np.sqrt(n) * X.mean(axis=1) / X.std(axis=1, ddof=1)
        signif = np.abs(t) > crit
        significativi += signif.sum()
        nulle_fra_significativi += (vera_nulla & signif).sum()
    esatto = pi0 * alfa / (pi0 * alfa + (1 - pi0) * potenza)
    print(f"{pi0:29.2f}  {esatto:7.4f}  {nulle_fra_significativi / significativi:10.4f}")

# Output:
# alpha = 0.05, power (n = 34, effect 0.5) = 0.8078
# share of true null hypotheses   exact   simulated  (among results with p < 0.05)
#                          0.50   0.0583      0.0577
#                          0.90   0.3578      0.3589
#                          0.99   0.8597      0.8592

If half of the hypotheses you test are true, 5.8% of the significant results are false alarms. If one in ten is true, they are 35.8%, and if one in a hundred, 86%: the simulation with two million hypotheses coincides with the formula within 0.11 percentage points. The p-value is the same (below 0.05) in all three worlds, but the probability that the result is a false alarm changes fifteen-fold. This is why a single p-value is not enough: you need the estimated effect with its interval, the prior plausibility of the hypothesis, and an independent confirmation.

Five mistakes that are really made

• Using the Normal in place of the t with small samples and estimated σ: the true level rises from 5% to 12% with five observations.
• Running Student with ttest_ind without noticing that it is the default, when the groups have different variances and sizes: the true level reaches 25.7%.
• Reading a non-significant result as proof that there is no difference: with 25 observations and an effect of 0.4 standard deviations the power is 0.4840.
• Stopping data collection as soon as p falls below 0.05 or running many tests and reporting the best one: the nominal level no longer holds.
• Reading the p-value as the probability that the null hypothesis is true: without the fraction of true hypotheses and without the power, that probability cannot be computed.

The five mistakes all come from the same swap: treating the guarantee of a test, which concerns the procedure repeated over many samples and under precise conditions, as a guarantee on the single result. For each one you saw above a number that shows by how much the guarantee breaks, and each one can be checked with a short program.

You can now set up a test (hypotheses, statistic, rule, level), choose between z, t, Student and Welch according to what you know about the data, compute with scipy.stats.nct how much data you need for a given power, and say what a p-value asserts and what it does not. In the article on ANOVA you will extend the comparison to more than two groups, where the multiple testing problem seen here becomes the central theme, and in the one on the chi-square you will test categorical data.

The same series also includes Simple and Multiple Linear Regression, Nonparametric Statistics: Wilcoxon and Mann-Whitney and Exploratory Data Analysis (EDA) with Python.