Chapter 3: Hypothesis Testing and Power Analysis | Inferential and Bayesian Statistics

3.1 The Framework of Hypothesis Testing

📘 Basic Structure of Hypothesis Testing

Null Hypothesis $ H_0 $: The claim that "there is no difference" or "there is no effect"

Alternative Hypothesis $ H_1 $: The claim that "there is a difference" or "there is an effect"

Type I Error: The error of rejecting $ H_0 $ when it is actually true (false positive)

$$ \alpha = P(\text{Type I Error}) = P(\text{reject } H_0 | H_0 \text{ is true}) $$

Type II Error: The error of failing to reject $ H_0 $ when it is actually false (false negative)

$$ \beta = P(\text{Type II Error}) = P(\text{fail to reject } H_0 | H_1 \text{ is true}) $$

Power: The probability of correctly detecting a true effect

$$ \text{Power} = 1 - \beta $$

💻 Code Example 1: z-test (test for a population mean)

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Scenario: Did the new manufacturing process improve the average strength over the conventional one?
# H0: μ = 450 (mean of the conventional process)
# H1: μ > 450 (one-sided test)

mu_0 = 450  # population mean under the null hypothesis
sigma = 30  # population standard deviation (known)
alpha = 0.05  # significance level

# Data from the new process
np.random.seed(42)
n = 25
mu_true = 465  # actual population mean (assumed unknown)
data = np.random.normal(mu_true, sigma, n)

# Computing the test statistic
x_bar = np.mean(data)
se = sigma / np.sqrt(n)
z_stat = (x_bar - mu_0) / se

# Computing the p-value (one-sided test)
p_value = 1 - stats.norm.cdf(z_stat)

# Rejection region
z_critical = stats.norm.ppf(1 - alpha)

print("=== z-test: one-sided test for a population mean ===")
print(f"Null hypothesis H0: μ = {mu_0}")
print(f"Alternative hypothesis H1: μ > {mu_0}")
print(f"Significance level: α = {alpha}")
print(f"\nSample size: {n}")
print(f"Sample mean: {x_bar:.2f}")
print(f"z-statistic: {z_stat:.4f}")
print(f"Critical value: z_{{{1-alpha}}} = {z_critical:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"\nDecision: {'Reject H0' if p_value < alpha else 'Cannot reject H0'}")
print(f"Interpretation: {'The new process significantly improved the strength' if p_value < alpha else 'No significant improvement can be confirmed'}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Distribution and rejection region under the null hypothesis
x = np.linspace(-4, 4, 300)
axes[0].plot(x, stats.norm.pdf(x), 'b-', linewidth=2, label='Distribution under H0')
axes[0].fill_between(x, 0, stats.norm.pdf(x), where=(x >= z_critical),
                      alpha=0.3, color='red', label=f'Rejection region (α={alpha})')
axes[0].axvline(z_stat, color='green', linestyle='--', linewidth=2,
                label=f'Observed z-statistic: {z_stat:.2f}')
axes[0].axvline(z_critical, color='red', linestyle=':', linewidth=2,
                label=f'Critical value: {z_critical:.2f}')
axes[0].set_xlabel('z-value')
axes[0].set_ylabel('Probability density')
axes[0].set_title('Visualization of the hypothesis test (one-sided)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Raw data and the distribution under the null hypothesis
x_data = np.linspace(mu_0 - 4*sigma, mu_0 + 4*sigma, 300)
axes[1].hist(data, bins=12, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Observed data')
axes[1].plot(x_data, stats.norm.pdf(x_data, mu_0, se),
             'r-', linewidth=2, label=f'H0: μ={mu_0} distribution')
axes[1].axvline(x_bar, color='blue', linestyle='--', linewidth=2,
                label=f'Sample mean: {x_bar:.1f}')
axes[1].axvline(mu_0, color='red', linestyle=':', linewidth=2,
                label=f'H0 mean: {mu_0}')
axes[1].set_xlabel('Strength [MPa]')
axes[1].set_ylabel('Density')
axes[1].set_title('Data vs. the null hypothesis')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

📌 The correct interpretation of the p-value
The p-value is "the probability of obtaining the observed data, or data more extreme, given that the null hypothesis is true." It is NOT "the probability that the null hypothesis is correct." p < 0.05 means "a rare phenomenon that occurs with less than 5% probability."

3.2 The t-test (one-sample, two-sample, paired)

💻 Code Example 2: The three patterns of the t-test

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# ===== One-sample t-test =====
print("=== One-sample t-test ===")
# Does the material achieve the target strength of 500 MPa?
mu_target = 500
n1 = 15
data_1sample = np.random.normal(510, 25, n1)

t_stat_1, p_value_1 = stats.ttest_1samp(data_1sample, mu_target)
print(f"H0: μ = {mu_target}")
print(f"Sample mean: {np.mean(data_1sample):.2f}")
print(f"t-statistic: {t_stat_1:.4f}, p-value: {p_value_1:.6f}")
print(f"Decision: {'Reject H0 (significantly different from the target)' if p_value_1 < 0.05 else 'Cannot reject H0'}\n")

# ===== Two-sample t-test (independent) =====
print("=== Two-sample t-test (independent) ===")
# Is there a difference in material strength between two suppliers?
n2a, n2b = 20, 22
data_2sample_A = np.random.normal(480, 20, n2a)
data_2sample_B = np.random.normal(495, 25, n2b)

# Welch's t-test (unequal variances)
t_stat_2, p_value_2 = stats.ttest_ind(data_2sample_A, data_2sample_B, equal_var=False)
print(f"H0: μA = μB")
print(f"Mean A: {np.mean(data_2sample_A):.2f}, Mean B: {np.mean(data_2sample_B):.2f}")
print(f"t-statistic: {t_stat_2:.4f}, p-value: {p_value_2:.6f}")
print(f"Decision: {'Reject H0 (significant difference)' if p_value_2 < 0.05 else 'Cannot reject H0'}\n")

# ===== Paired t-test =====
print("=== Paired t-test ===")
# Did the hardness of the material change before and after heat treatment?
n3 = 12
before = np.random.normal(65, 5, n3)
after = before + np.random.normal(3, 2, n3)  # a mean increase of 3

t_stat_3, p_value_3 = stats.ttest_rel(after, before)
diff = after - before
print(f"H0: μdiff = 0 (no change before and after treatment)")
print(f"Mean difference: {np.mean(diff):.2f}")
print(f"t-statistic: {t_stat_3:.4f}, p-value: {p_value_3:.6f}")
print(f"Decision: {'Reject H0 (significant change)' if p_value_3 < 0.05 else 'Cannot reject H0'}\n")

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# One-sample t-test
axes[0].hist(data_1sample, bins=8, density=True, alpha=0.7,
             color='skyblue', edgecolor='black')
axes[0].axvline(np.mean(data_1sample), color='blue', linestyle='--',
                linewidth=2, label=f'Sample mean: {np.mean(data_1sample):.1f}')
axes[0].axvline(mu_target, color='red', linestyle='--',
                linewidth=2, label=f'Target value: {mu_target}')
axes[0].set_xlabel('Strength [MPa]')
axes[0].set_ylabel('Density')
axes[0].set_title(f'One-sample t-test\np={p_value_1:.4f}')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Two-sample t-test
bp = axes[1].boxplot([data_2sample_A, data_2sample_B], labels=['Supplier A', 'Supplier B'],
                      patch_artist=True)
bp['boxes'][0].set_facecolor('lightblue')
bp['boxes'][1].set_facecolor('lightgreen')
axes[1].set_ylabel('Strength [MPa]')
axes[1].set_title(f'Two-sample t-test\np={p_value_2:.4f}')
axes[1].grid(True, alpha=0.3, axis='y')

# Paired t-test
x_pos = np.arange(n3)
axes[2].plot(x_pos, before, 'bo-', label='Before treatment', alpha=0.6)
axes[2].plot(x_pos, after, 'rs-', label='After treatment', alpha=0.6)
for i in range(n3):
    axes[2].plot([i, i], [before[i], after[i]], 'k-', alpha=0.3)
axes[2].set_xlabel('Sample number')
axes[2].set_ylabel('Hardness [HV]')
axes[2].set_title(f'Paired t-test\np={p_value_3:.4f}')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

3.3 The Chi-Square Test (goodness-of-fit and independence)

💻 Code Example 3: The chi-square test

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# ===== Goodness-of-fit test =====
print("=== Chi-square goodness-of-fit test ===")
# Testing the fairness of a die
observed = np.array([48, 52, 45, 50, 55, 50])  # observed frequencies
expected = np.array([50, 50, 50, 50, 50, 50])  # expected frequencies (1/6 for each face if fair)

chi2_stat, p_value_gof = stats.chisquare(observed, expected)
print(f"H0: the die is fair")
print(f"Observed frequencies: {observed}")
print(f"Expected frequencies: {expected}")
print(f"χ² statistic: {chi2_stat:.4f}")
print(f"p-value: {p_value_gof:.6f}")
print(f"Decision: {'Reject H0 (unfair)' if p_value_gof < 0.05 else 'Cannot reject H0'}\n")

# ===== Test of independence =====
print("=== Chi-square test of independence ===")
# Is there an association between material type and defect occurrence?
# 2×3 contingency table
contingency_table = np.array([
    [15, 25, 10],  # Material A: good, minor defect, major defect
    [20, 18, 12]   # Material B
])

chi2_stat_ind, p_value_ind, dof, expected_ind = stats.chi2_contingency(contingency_table)

print(f"H0: material type and defect severity are independent")
print(f"Observed frequencies:\n{contingency_table}")
print(f"Expected frequencies:\n{expected_ind}")
print(f"χ² statistic: {chi2_stat_ind:.4f}")
print(f"Degrees of freedom: {dof}")
print(f"p-value: {p_value_ind:.6f}")
print(f"Decision: {'Reject H0 (association exists)' if p_value_ind < 0.05 else 'Cannot reject H0'}\n")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Goodness-of-fit test
x = np.arange(1, 7)
width = 0.35
axes[0].bar(x - width/2, observed, width, label='Observed', color='skyblue', edgecolor='black')
axes[0].bar(x + width/2, expected, width, label='Expected', color='orange', edgecolor='black', alpha=0.7)
axes[0].set_xlabel('Die face')
axes[0].set_ylabel('Frequency')
axes[0].set_title(f'Goodness-of-fit test\nχ²={chi2_stat:.2f}, p={p_value_gof:.4f}')
axes[0].set_xticks(x)
axes[0].legend()
axes[0].grid(True, alpha=0.3, axis='y')

# Test of independence
categories = ['Good', 'Minor defect', 'Major defect']
materials = ['Material A', 'Material B']
x = np.arange(len(categories))
width = 0.35

axes[1].bar(x - width/2, contingency_table[0], width, label='Material A',
            color='lightblue', edgecolor='black')
axes[1].bar(x + width/2, contingency_table[1], width, label='Material B',
            color='lightgreen', edgecolor='black')
axes[1].set_xlabel('Quality category')
axes[1].set_ylabel('Frequency')
axes[1].set_title(f'Test of independence\nχ²={chi2_stat_ind:.2f}, p={p_value_ind:.4f}')
axes[1].set_xticks(x)
axes[1].set_xticklabels(categories)
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

3.4 The F-test (test for equality of variances)

💻 Code Example 4: The F-test

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# Are the variances of two processes equal?
n1, n2 = 20, 25
sigma1, sigma2 = 15, 22  # true standard deviations

data1 = np.random.normal(100, sigma1, n1)
data2 = np.random.normal(100, sigma2, n2)

var1 = np.var(data1, ddof=1)
var2 = np.var(data2, ddof=1)

# F-statistic (larger variance / smaller variance)
F_stat = max(var1, var2) / min(var1, var2)
df1 = n1 - 1 if var1 > var2 else n2 - 1
df2 = n2 - 1 if var1 > var2 else n1 - 1

# p-value for the two-sided test
p_value = 2 * min(stats.f.cdf(F_stat, df1, df2),
                  1 - stats.f.cdf(F_stat, df1, df2))

print("=== F-test: test for equality of variances ===")
print(f"H0: σ₁² = σ₂² (equal variances)")
print(f"H1: σ₁² ≠ σ₂²")
print(f"\nData 1: n={n1}, s²={var1:.2f}")
print(f"Data 2: n={n2}, s²={var2:.2f}")
print(f"\nF-statistic: {F_stat:.4f}")
print(f"Degrees of freedom: ({df1}, {df2})")
print(f"p-value: {p_value:.6f}")
print(f"\nDecision: {'Reject H0 (variances differ)' if p_value < 0.05 else 'Cannot reject H0 (equal variances can be assumed)'}")

# Levene's test (more robust)
levene_stat, levene_p = stats.levene(data1, data2)
print(f"\n[Levene's test] (more robust)")
print(f"Statistic: {levene_stat:.4f}, p-value: {levene_p:.6f}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Distributions of the data
axes[0].hist(data1, bins=12, density=True, alpha=0.6,
             color='skyblue', edgecolor='black', label=f'Data 1 (s²={var1:.1f})')
axes[0].hist(data2, bins=12, density=True, alpha=0.6,
             color='lightgreen', edgecolor='black', label=f'Data 2 (s²={var2:.1f})')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Density')
axes[0].set_title('Distributions of the two datasets')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# F-distribution and the observed F-statistic
x = np.linspace(0, 5, 300)
axes[1].plot(x, stats.f.pdf(x, df1, df2), 'b-', linewidth=2,
             label=f'F-distribution ({df1}, {df2})')
axes[1].axvline(F_stat, color='red', linestyle='--', linewidth=2,
                label=f'Observed F-statistic: {F_stat:.2f}')

# Critical value (two-sided 5%)
f_crit_upper = stats.f.ppf(0.975, df1, df2)
axes[1].axvline(f_crit_upper, color='orange', linestyle=':', linewidth=2,
                label=f'Critical value (upper 2.5%): {f_crit_upper:.2f}')
axes[1].fill_between(x, 0, stats.f.pdf(x, df1, df2),
                      where=(x >= f_crit_upper), alpha=0.3, color='red',
                      label='Rejection region')
axes[1].set_xlabel('F-value')
axes[1].set_ylabel('Probability density')
axes[1].set_title(f'Visualization of the F-test\np={p_value:.4f}')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

3.5 Power Analysis and Sample Size Calculation

📘 Power

Power is the probability of correctly detecting a true effect:

$$ \text{Power} = 1 - \beta = P(\text{reject } H_0 | H_1 \text{ is true}) $$

Factors that affect power:

Effect Size: The magnitude of the true difference
Sample size (n): The larger it is, the higher the power
Significance level (α): The larger it is, the higher the power (but Type I error increases)
Data variance: The smaller it is, the higher the power

💻 Code Example 5: Power analysis and sample size calculation

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

# Parameters for the power analysis
mu0 = 100  # mean under the null hypothesis
mu1 = 105  # mean under the alternative hypothesis
sigma = 10  # standard deviation
alpha = 0.05

# Cohen's d (effect size)
effect_size = (mu1 - mu0) / sigma
print("=== Power analysis ===")
print(f"Effect size Cohen's d: {effect_size:.3f}")
print(f"Interpretation: {'small' if effect_size < 0.5 else 'medium' if effect_size < 0.8 else 'large'}")

# Relationship between sample size and power
sample_sizes = np.arange(10, 201, 5)
powers = []

power_analysis = TTestIndPower()
for n in sample_sizes:
    power = power_analysis.power(effect_size, n, alpha, alternative='two-sided')
    powers.append(power)

# Sample size to achieve the target power of 0.8
target_power = 0.8
required_n = power_analysis.solve_power(effect_size, power=target_power,
                                        alpha=alpha, alternative='two-sided')

print(f"\nSample size to achieve the target power {target_power}: {int(np.ceil(required_n))}")

# Relationship between effect size and power
effect_sizes = np.linspace(0.1, 1.5, 50)
powers_by_effect = []
n_fixed = 50

for es in effect_sizes:
    power = power_analysis.power(es, n_fixed, alpha, alternative='two-sided')
    powers_by_effect.append(power)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Distributions of the null and alternative hypotheses
x = np.linspace(80, 120, 300)
se_example = sigma / np.sqrt(50)
axes[0, 0].plot(x, stats.norm.pdf(x, mu0, se_example), 'b-',
                linewidth=2, label=f'H0: μ={mu0}')
axes[0, 0].plot(x, stats.norm.pdf(x, mu1, se_example), 'r-',
                linewidth=2, label=f'H1: μ={mu1}')

# Critical values
z_crit = stats.norm.ppf(1 - alpha/2)
x_crit_upper = mu0 + z_crit * se_example
x_crit_lower = mu0 - z_crit * se_example

axes[0, 0].axvline(x_crit_upper, color='green', linestyle='--',
                   linewidth=2, label='Critical value')
axes[0, 0].axvline(x_crit_lower, color='green', linestyle='--', linewidth=2)

# β region (Type II error)
x_beta = np.linspace(x_crit_lower, x_crit_upper, 100)
axes[0, 0].fill_between(x_beta, 0, stats.norm.pdf(x_beta, mu1, se_example),
                         alpha=0.3, color='orange', label='β (Type II error)')

# Power region
x_power_lower = np.linspace(80, x_crit_lower, 100)
x_power_upper = np.linspace(x_crit_upper, 120, 100)
axes[0, 0].fill_between(x_power_lower, 0, stats.norm.pdf(x_power_lower, mu1, se_example),
                         alpha=0.3, color='green', label='Power')
axes[0, 0].fill_between(x_power_upper, 0, stats.norm.pdf(x_power_upper, mu1, se_example),
                         alpha=0.3, color='green')

axes[0, 0].set_xlabel('Value')
axes[0, 0].set_ylabel('Probability density')
axes[0, 0].set_title(f'Conceptual diagram of power (n=50, Power={powers[8]:.3f})')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Sample size and power
axes[0, 1].plot(sample_sizes, powers, 'b-', linewidth=2)
axes[0, 1].axhline(target_power, color='red', linestyle='--',
                   linewidth=2, label=f'Target power: {target_power}')
axes[0, 1].axvline(required_n, color='green', linestyle='--',
                   linewidth=2, label=f'Required n: {int(np.ceil(required_n))}')
axes[0, 1].set_xlabel('Sample size (per group)')
axes[0, 1].set_ylabel('Power (1-β)')
axes[0, 1].set_title(f'Sample size and power (d={effect_size:.2f})')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Effect size and power
axes[1, 0].plot(effect_sizes, powers_by_effect, 'b-', linewidth=2)
axes[1, 0].axhline(target_power, color='red', linestyle='--',
                   linewidth=2, label=f'Target power: {target_power}')
axes[1, 0].axvline(effect_size, color='orange', linestyle='--',
                   linewidth=2, label=f'Current effect size: {effect_size:.2f}')
axes[1, 0].set_xlabel("Effect size (Cohen's d)")
axes[1, 0].set_ylabel('Power (1-β)')
axes[1, 0].set_title(f'Effect size and power (n={n_fixed})')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Trade-off between α and β
alphas = np.linspace(0.01, 0.2, 50)
betas_at_alphas = []

for a in alphas:
    power = power_analysis.power(effect_size, n_fixed, a, alternative='two-sided')
    betas_at_alphas.append(1 - power)

axes[1, 1].plot(alphas, betas_at_alphas, 'b-', linewidth=2, label='β')
axes[1, 1].plot(alphas, alphas, 'r--', linewidth=2, label='α')
axes[1, 1].axvline(0.05, color='green', linestyle=':', linewidth=2,
                   label='α=0.05 (convention)')
axes[1, 1].set_xlabel('α (Type I error rate)')
axes[1, 1].set_ylabel('Error rate')
axes[1, 1].set_title(f'Trade-off between α and β (n={n_fixed}, d={effect_size:.2f})')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Practical advice
print("\n=== Guidelines for sample size design ===")
print(f"Small effect (d=0.2): n ≈ {int(np.ceil(power_analysis.solve_power(0.2, power=0.8, alpha=0.05)))}")
print(f"Medium effect (d=0.5): n ≈ {int(np.ceil(power_analysis.solve_power(0.5, power=0.8, alpha=0.05)))}")
print(f"Large effect (d=0.8): n ≈ {int(np.ceil(power_analysis.solve_power(0.8, power=0.8, alpha=0.05)))}")

3.6 The Multiple Comparisons Problem and Correction Methods

📘 The Multiple Comparisons Problem

When multiple hypothesis tests are performed simultaneously, the probability of making at least one Type I error (the Family-Wise Error Rate, FWER) increases:

$$ \text{FWER} = 1 - (1-\alpha)^m $$

Here $ m $ is the number of tests.

Bonferroni correction: Set the significance level of each test to $ \alpha/m $

Holm's method: A step-down method (higher power)

FDR control: The Benjamini-Hochberg method (controls the False Discovery Rate)

💻 Code Example 6: Multiple comparison corrections

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.multitest import multipletests

np.random.seed(42)

# Test 10 material properties simultaneously
n_tests = 10
n_samples = 30
alpha = 0.05

# 9 have no difference, only 1 has a true difference
p_values = []
for i in range(n_tests):
    if i == 5:  # only the 6th has a true difference
        data1 = np.random.normal(100, 15, n_samples)
        data2 = np.random.normal(110, 15, n_samples)
    else:  # the others have no difference
        data1 = np.random.normal(100, 15, n_samples)
        data2 = np.random.normal(100, 15, n_samples)

    _, p = stats.ttest_ind(data1, data2)
    p_values.append(p)

p_values = np.array(p_values)

print("=== The multiple comparisons problem ===")
print(f"Number of tests: {n_tests}")
print(f"Significance level of each test: α = {alpha}")
print(f"Theoretical FWER: {1 - (1-alpha)**n_tests:.4f}")
print(f"\nOriginal p-values:\n{p_values}")
print(f"\nSignificant without correction (α={alpha}): {np.sum(p_values < alpha)}")

# Applying the various correction methods
methods = ['bonferroni', 'holm', 'fdr_bh']
method_names = ['Bonferroni', 'Holm', 'Benjamini-Hochberg']

results = {}
for method, name in zip(methods, method_names):
    reject, p_corrected, _, _ = multipletests(p_values, alpha=alpha, method=method)
    results[name] = {'reject': reject, 'p_corrected': p_corrected}
    print(f"\n[{name} method]")
    print(f"  Corrected p-values: {p_corrected}")
    print(f"  Judged significant: {np.sum(reject)} (positions: {np.where(reject)[0]})")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Original p-values
axes[0, 0].bar(range(n_tests), p_values, color='skyblue', edgecolor='black')
axes[0, 0].axhline(alpha, color='red', linestyle='--', linewidth=2,
                   label=f'α = {alpha}')
axes[0, 0].set_xlabel('Test number')
axes[0, 0].set_ylabel('p-value')
axes[0, 0].set_title('Original p-values (no correction)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3, axis='y')

# Bonferroni correction
reject_bonf = results['Bonferroni']['reject']
colors = ['red' if r else 'skyblue' for r in reject_bonf]
axes[0, 1].bar(range(n_tests), p_values, color=colors, edgecolor='black')
axes[0, 1].axhline(alpha/n_tests, color='orange', linestyle='--',
                   linewidth=2, label=f'Bonferroni α = {alpha/n_tests:.4f}')
axes[0, 1].set_xlabel('Test number')
axes[0, 1].set_ylabel('p-value')
axes[0, 1].set_title('Bonferroni correction')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3, axis='y')

# Holm correction
reject_holm = results['Holm']['reject']
colors = ['red' if r else 'skyblue' for r in reject_holm]
axes[1, 0].bar(range(n_tests), p_values, color=colors, edgecolor='black')
# Display the stepwise thresholds of Holm's method
sorted_idx = np.argsort(p_values)
for rank, idx in enumerate(sorted_idx):
    holm_alpha = alpha / (n_tests - rank)
    axes[1, 0].plot(idx, holm_alpha, 'o', color='orange', markersize=8)

axes[1, 0].set_xlabel('Test number')
axes[1, 0].set_ylabel('p-value')
axes[1, 0].set_title('Holm correction (stepwise thresholds)')
axes[1, 0].grid(True, alpha=0.3, axis='y')

# FDR (Benjamini-Hochberg)
reject_fdr = results['Benjamini-Hochberg']['reject']
colors = ['red' if r else 'skyblue' for r in reject_fdr]
sorted_p = np.sort(p_values)
bh_thresholds = [(i+1)/n_tests * alpha for i in range(n_tests)]

axes[1, 1].bar(range(n_tests), sorted_p, color=['red' if sorted_p[i] < bh_thresholds[i] else 'skyblue' for i in range(n_tests)],
               edgecolor='black', alpha=0.7, label='p-value')
axes[1, 1].plot(range(n_tests), bh_thresholds, 'o-', color='orange',
                linewidth=2, markersize=6, label='BH threshold')
axes[1, 1].set_xlabel('Test number (sorted by p-value)')
axes[1, 1].set_ylabel('p-value')
axes[1, 1].set_title('Benjamini-Hochberg (FDR control)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

# Summary
print("\n=== Comparison of correction methods ===")
print(f"No correction: {np.sum(p_values < alpha)} significant")
print(f"Bonferroni: {np.sum(results['Bonferroni']['reject'])} significant (most conservative)")
print(f"Holm: {np.sum(results['Holm']['reject'])} significant")
print(f"Benjamini-Hochberg: {np.sum(results['Benjamini-Hochberg']['reject'])} significant (higher power)")

3.7 Hypothesis Testing in Quality Control

💻 Code Example 7: Comprehensive testing of quality control data

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Quality control scenario: quality data from a production line
np.random.seed(42)

# Conventional process (baseline)
baseline_mean = 500
baseline_std = 15
n_baseline = 50
baseline_data = np.random.normal(baseline_mean, baseline_std, n_baseline)

# Improved process (under evaluation)
improved_mean = 505  # improved mean
improved_std = 12    # reduced variability
n_improved = 60
improved_data = np.random.normal(improved_mean, improved_std, n_improved)

print("=== Statistical testing in quality control ===")

# 1. Test of the mean (t-test)
print("\n[1. Test of mean strength improvement]")
t_stat, p_value_mean = stats.ttest_ind(improved_data, baseline_data, equal_var=False)
print(f"Conventional process mean: {np.mean(baseline_data):.2f} MPa")
print(f"Improved process mean: {np.mean(improved_data):.2f} MPa")
print(f"t-statistic: {t_stat:.4f}, p-value: {p_value_mean:.6f}")
print(f"Decision: {'Significant improvement' if p_value_mean < 0.05 and t_stat > 0 else 'No significant improvement'}")

# 2. Test of variance (F-test)
print("\n[2. Test of variability improvement]")
var_baseline = np.var(baseline_data, ddof=1)
var_improved = np.var(improved_data, ddof=1)
F_stat = var_baseline / var_improved
df1, df2 = n_baseline - 1, n_improved - 1
p_value_var = stats.f.sf(F_stat, df1, df2)  # one-sided test (reduction in variability)

print(f"Conventional process variance: {var_baseline:.2f}")
print(f"Improved process variance: {var_improved:.2f}")
print(f"F-statistic: {F_stat:.4f}, p-value: {p_value_var:.6f}")
print(f"Decision: {'Variability significantly reduced' if p_value_var < 0.05 else 'No significant reduction in variability'}")

# 3. Test of the out-of-spec rate (test of proportions)
print("\n[3. Test of the out-of-spec rate]")
spec_lower = 470  # lower specification limit
spec_upper = 530  # upper specification limit

defect_baseline = np.sum((baseline_data < spec_lower) | (baseline_data > spec_upper))
defect_improved = np.sum((improved_data < spec_lower) | (improved_data > spec_upper))

p_baseline = defect_baseline / n_baseline
p_improved = defect_improved / n_improved

# Test of two proportions
pooled_p = (defect_baseline + defect_improved) / (n_baseline + n_improved)
se_pool = np.sqrt(pooled_p * (1 - pooled_p) * (1/n_baseline + 1/n_improved))
z_prop = (p_baseline - p_improved) / se_pool if se_pool > 0 else 0
p_value_prop = stats.norm.sf(z_prop)  # one-sided test

print(f"Conventional process out-of-spec rate: {p_baseline:.4f} ({defect_baseline}/{n_baseline})")
print(f"Improved process out-of-spec rate: {p_improved:.4f} ({defect_improved}/{n_improved})")
print(f"z-statistic: {z_prop:.4f}, p-value: {p_value_prop:.6f}")
print(f"Decision: {'Out-of-spec rate significantly improved' if p_value_prop < 0.05 else 'No significant improvement'}")

# 4. Comparison of process capability indices
print("\n[4. Process capability indices Cp, Cpk]")
spec_range = spec_upper - spec_lower

Cp_baseline = spec_range / (6 * np.std(baseline_data, ddof=1))
Cpk_baseline_lower = (np.mean(baseline_data) - spec_lower) / (3 * np.std(baseline_data, ddof=1))
Cpk_baseline_upper = (spec_upper - np.mean(baseline_data)) / (3 * np.std(baseline_data, ddof=1))
Cpk_baseline = min(Cpk_baseline_lower, Cpk_baseline_upper)

Cp_improved = spec_range / (6 * np.std(improved_data, ddof=1))
Cpk_improved_lower = (np.mean(improved_data) - spec_lower) / (3 * np.std(improved_data, ddof=1))
Cpk_improved_upper = (spec_upper - np.mean(improved_data)) / (3 * np.std(improved_data, ddof=1))
Cpk_improved = min(Cpk_improved_lower, Cpk_improved_upper)

print(f"Conventional process: Cp={Cp_baseline:.3f}, Cpk={Cpk_baseline:.3f}")
print(f"Improved process: Cp={Cp_improved:.3f}, Cpk={Cpk_improved:.3f}")
print(f"Evaluation: Cpk ≥ 1.33 is excellent, ≥ 1.00 is acceptable")

# Visualization
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3)

# 1. Comparison of the data distributions
ax1 = fig.add_subplot(gs[0, :])
ax1.hist(baseline_data, bins=20, density=True, alpha=0.6,
         color='skyblue', edgecolor='black', label='Conventional process')
ax1.hist(improved_data, bins=20, density=True, alpha=0.6,
         color='lightgreen', edgecolor='black', label='Improved process')
ax1.axvline(spec_lower, color='red', linestyle='--', linewidth=2, label='Lower spec limit')
ax1.axvline(spec_upper, color='red', linestyle='--', linewidth=2, label='Upper spec limit')
ax1.axvline(np.mean(baseline_data), color='blue', linestyle=':', linewidth=2)
ax1.axvline(np.mean(improved_data), color='green', linestyle=':', linewidth=2)
ax1.set_xlabel('Strength [MPa]')
ax1.set_ylabel('Density')
ax1.set_title(f'Effect of process improvement\nMean: {np.mean(baseline_data):.1f} → {np.mean(improved_data):.1f} (p={p_value_mean:.4f}), '
              f'SD: {np.std(baseline_data, ddof=1):.1f} → {np.std(improved_data, ddof=1):.1f}')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Box plot
ax2 = fig.add_subplot(gs[1, 0])
bp = ax2.boxplot([baseline_data, improved_data], labels=['Conventional', 'Improved'],
                  patch_artist=True)
bp['boxes'][0].set_facecolor('lightblue')
bp['boxes'][1].set_facecolor('lightgreen')
ax2.axhline(spec_lower, color='red', linestyle='--', linewidth=1.5, alpha=0.7)
ax2.axhline(spec_upper, color='red', linestyle='--', linewidth=1.5, alpha=0.7)
ax2.set_ylabel('Strength [MPa]')
ax2.set_title('Comparison of variability')
ax2.grid(True, alpha=0.3, axis='y')

# 3. Q-Q plot (normality check)
ax3 = fig.add_subplot(gs[1, 1])
stats.probplot(improved_data, dist="norm", plot=ax3)
ax3.set_title('Normal Q-Q plot of the improved process')
ax3.grid(True, alpha=0.3)

# 4. Comparison of process capability indices
ax4 = fig.add_subplot(gs[2, 0])
metrics = ['Cp', 'Cpk']
baseline_metrics = [Cp_baseline, Cpk_baseline]
improved_metrics = [Cp_improved, Cpk_improved]
x = np.arange(len(metrics))
width = 0.35

ax4.bar(x - width/2, baseline_metrics, width, label='Conventional',
        color='lightblue', edgecolor='black')
ax4.bar(x + width/2, improved_metrics, width, label='Improved',
        color='lightgreen', edgecolor='black')
ax4.axhline(1.33, color='green', linestyle='--', linewidth=2, label='Excellent line')
ax4.axhline(1.00, color='orange', linestyle='--', linewidth=2, label='Acceptable line')
ax4.set_xticks(x)
ax4.set_xticklabels(metrics)
ax4.set_ylabel('Process capability index')
ax4.set_title('Improvement in process capability')
ax4.legend()
ax4.grid(True, alpha=0.3, axis='y')

# 5. Comparison of out-of-spec rates
ax5 = fig.add_subplot(gs[2, 1])
categories = ['Conventional process', 'Improved process']
defect_rates = [p_baseline * 100, p_improved * 100]
bars = ax5.bar(categories, defect_rates, color=['lightblue', 'lightgreen'],
               edgecolor='black')
for bar, rate in zip(bars, defect_rates):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height,
             f'{rate:.2f}%', ha='center', va='bottom', fontsize=12)
ax5.set_ylabel('Out-of-spec rate [%]')
ax5.set_title(f'Improvement in out-of-spec rate (p={p_value_prop:.4f})')
ax5.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print("\n=== Overall evaluation ===")
improvements = []
if p_value_mean < 0.05 and t_stat > 0:
    improvements.append("✓ Mean strength significantly improved")
if p_value_var < 0.05:
    improvements.append("✓ Variability significantly reduced")
if p_value_prop < 0.05:
    improvements.append("✓ Out-of-spec rate significantly improved")
if Cpk_improved > Cpk_baseline:
    improvements.append(f"✓ Process capability index improved ({Cpk_baseline:.2f} → {Cpk_improved:.2f})")

if improvements:
    print("The process improvement is statistically significant:")
    for imp in improvements:
        print(f"  {imp}")
else:
    print("No statistically significant improvement could be confirmed")

📝 Exercises

Explain, with concrete examples, how to choose between a one-sided test and a two-sided test.
Discuss how a p-value of 0.06 should be interpreted.
Calculate the sample size required to achieve a power of 0.8 for effect sizes d=0.3, 0.5, and 0.8.
When performing five hypothesis tests, what is the significance level of each test after Bonferroni correction?

Summary

Hypothesis testing is a framework for statistically evaluating a null hypothesis
Balancing Type I error (α) and Type II error (β) is important
The p-value is the probability of extremeness given that the null hypothesis is true, not the probability of the null hypothesis being true or false
The t-test is used to compare population means, and the F-test to compare variances
The chi-square test is useful for frequency data and tests of independence
Power analysis is an essential step in experimental design
For multiple comparisons, choose an appropriate correction method (Bonferroni, Holm, FDR)
In quality control, evaluation combines multiple statistical metrics