统计假设检验

概述假设检验提供了一个框架，通过检验观察到的差异是否具有统计学意义或是否由偶然因素引起，从而做出数据驱动的决策。

测试框架

零假设 (H0): 不存在效应或差异
备择假设 (H1): 存在效应或差异
显著性水平 (α): 拒绝H0的阈值（通常为0.05）
P值: 如果H0为真，观察到数据的概率

常见测试

T检验: 比较两组之间的均值
方差分析 (ANOVA): 比较多个组的均值
卡方检验: 测试分类变量的独立性
曼-惠特尼U检验: T检验的非参数替代方法
克鲁斯卡尔-瓦利斯检验: ANOVA的非参数替代方法

使用Python实现 import pandas as pd import numpy as np from scipy import stats import matplotlib.pyplot as plt

示例数据

group_a = np.random.normal(100, 15, 50) # 均值=100, SD=15 group_b = np.random.normal(105, 15, 50) # 均值=105, SD=15

测试1: 独立样本T检验

t_stat, p_value = stats.ttest_ind(group_a, group_b) print(f"T检验: t={t_stat:.4f}, p值={p_value:.4f}") if p_value < 0.05: print(“拒绝零假设: 组间存在显著差异”) else: print(“未能拒绝零假设: 没有显著差异”)

测试2: 配对T检验（同一受试者，两种条件）

before = np.array([85, 90, 88, 92, 87, 89, 91, 86, 88, 90]) after = np.array([92, 95, 91, 98, 94, 96, 99, 93, 95, 97])

t_stat, p_value = stats.ttest_rel(before, after) print(f" 配对T检验: t={t_stat:.4f}, p值={p_value:.4f}")

测试3: 单因素方差分析（多组）

group1 = np.random.normal(100, 10, 30) group2 = np.random.normal(105, 10, 30) group3 = np.random.normal(102, 10, 30)

f_stat, p_value = stats.f_oneway(group1, group2, group3) print(f" 方差分析: F={f_stat:.4f}, p值={p_value:.4f}")

测试4: 卡方检验（分类变量）

创建列联表

contingency = np.array([ [50, 30], # 控制组: 成功, 失败 [45, 35] # 治疗组: 成功, 失败 ])

chi2, p_value, dof, expected = stats.chi2_contingency(contingency) print(f" 卡方检验: χ²={chi2:.4f}, p值={p_value:.4f}")

测试5: 曼-惠特尼U检验（非参数）

u_stat, p_value = stats.mannwhitneyu(group_a, group_b) print(f" 曼-惠特尼U检验: U={u_stat:.4f}, p值={p_value:.4f}")

可视化

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

分布比较

axes[0, 0].hist(group_a, alpha=0.5, label=‘组A’, bins=20) axes[0, 0].hist(group_b, alpha=0.5, label=‘组B’, bins=20) axes[0, 0].set_title(‘组分布’) axes[0, 0].legend()

Q-Q图正态性检验

stats.probplot(group_a, dist=“norm”, plot=axes[0, 1]) axes[0, 1].set_title(‘Q-Q图（组A）’)

前后比较

axes[1, 0].plot(before, ‘o-’, label=‘前’, alpha=0.7) axes[1, 0].plot(after, ‘s-’, label=‘后’, alpha=0.7) axes[1, 0].set_title(‘配对比较’) axes[1, 0].legend()

效应量（科恩d）

cohens_d = (np.mean(group_a) - np.mean(group_b)) / np.sqrt( ((len(group_a)-1)*np.var(group_a, ddof=1) + (len(group_b)-1)*np.var(group_b, ddof=1)) / (len(group_a) + len(group_b) - 2) ) axes[1, 1].text(0.5, 0.5, f"科恩的d = {cohens_d:.4f}", ha=‘center’, va=‘center’, fontsize=14) axes[1, 1].axis(‘off’)

plt.tight_layout() plt.show()

正态性检验（Shapiro-Wilk）

stat, p = stats.shapiro(group_a) print(f" Shapiro-Wilk正态性检验: W={stat:.4f}, p值={p:.4f}")

效应量计算

def calculate_effect_size(group1, group2): n1, n2 = len(group1), len(group2) var1, var2 = np.var(group1, ddof=1), np.var(group2, ddof=1) pooled_std = np.sqrt(((n1-1)*var1 + (n2-1)*var2) / (n1+n2-2)) cohens_d = (np.mean(group1) - np.mean(group2)) / pooled_std return cohens_d

effect_size = calculate_effect_size(group_a, group_b) print(f"效应量（科恩的d）: {effect_size:.4f}")

置信区间

from scipy.stats import t as t_dist def calculate_ci(data, confidence=0.95): n = len(data) mean = np.mean(data) se = np.std(data, ddof=1) / np.sqrt(n) margin = t_dist.ppf((1 + confidence) / 2, n - 1) * se return mean - margin, mean + margin

ci = calculate_ci(group_a) print(f"95% CI for Group A: ({ci[0]:.2f}, {ci[1]:.2f})")

其他测试和可视化

测试6: Levene’s test for equal variances

stat_levene, p_levene = stats.levene(group_a, group_b) print(f" Levene’s Test for Equal Variance:“) print(f"Statistic: {stat_levene:.4f}, P-value: {p_levene:.4f}”)

测试7: Welch’s t-test (doesn’t assume equal variance)

t_stat_welch, p_welch = stats.ttest_ind(group_a, group_b, equal_var=False) print(f" Welch’s t-test (unequal variance):“) print(f"t-stat: {t_stat_welch:.4f}, p-value: {p_welch:.4f}”)

功效分析

from scipy.stats import nct def calculate_power(effect_size, sample_size, alpha=0.05): t_critical = stats.t.ppf(1 - alpha/2, 2sample_size - 2) ncp = effect_size * np.sqrt(sample_size / 2) power = 1 - stats.nct.cdf(t_critical, 2sample_size - 2, ncp) return power

power = calculate_power(abs(effect_size), len(group_a)) print(f" Statistical Power: {power:.2%}")

Bootstrap置信区间

def bootstrap_ci(data, n_bootstrap=10000, ci=95): bootstrap_means = [] for _ in range(n_bootstrap): sample = np.random.choice(data, size=len(data), replace=True) bootstrap_means.append(np.mean(sample)) lower = np.percentile(bootstrap_means, (100-ci)/2) upper = np.percentile(bootstrap_means, ci + (100-ci)/2) return lower, upper

boot_ci = bootstrap_ci(group_a) print(f" Bootstrap 95% CI for Group A: ({boot_ci[0]:.2f}, {boot_ci[1]:.2f})")

多重测试校正（Bonferroni）

num_tests = 4 bonferroni_alpha = 0.05 / num_tests print(f" Bonferroni Corrected Alpha: {bonferroni_alpha:.4f}“) print(f"Use this threshold for {num_tests} tests”)

测试8: Kruskal-Wallis test (non-parametric ANOVA)

h_stat, p_kw = stats.kruskal(group1, group2, group3) print(f" Kruskal-Wallis Test (non-parametric ANOVA):“) print(f"H-statistic: {h_stat:.4f}, p-value: {p_kw:.4f}”)

效应量ANOVA

f_stat, p_anova = stats.f_oneway(group1, group2, group3)

计算eta平方

grand_mean = np.mean([group1, group2, group3]) ss_between = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in [group1, group2, group3]) ss_total = sum((x - grand_mean)**2 for g in [group1, group2, group3] for x in g) eta_squared = ss_between / ss_total print(f" Effect Size (Eta-squared): {eta_squared:.4f}")