🧩 多重比较

🚨 多重比较问题的本质

🎯 药物试验的统计陷阱

一家制药公司开发了新药，要验证疗效。研究人员设计了实验：

• 对比新药和安慰剂
• 测量10个指标：血压、心率、血糖...
• 每个指标都做显著性检验（α=0.05）

结果：有2个指标显示显著差异！研究人员兴奋地准备发论文...

等等！这里有个陷阱。

假设10个原假设均为真，各检验相互独立，且每次检验的假阳性概率均为5%，那么至少出现一次假阳性的概率是： 1 - (1-0.05)^10 ≈ 40%

这就是多重比较问题——检验次数越多，犯错概率越大！

多重比较的数学基础

Bonferroni不等式： $$P(\text{至少一次拒绝真实原假设}) \leq \sum_{i=1}^k P(\text{拒绝第i个真实原假设})$$

如果每次检验的显著性水平为α，进行k次检验，则： $$\text{FWER} \leq k \times \alpha$$

📊 常用多重比较方法

1. 最小显著差异法（LSD, Least Significant Difference）

$$LSD = t_{\alpha/2, df_{error}} \times SE_{diff}$$

优点：检验力高，容易发现真实差异 缺点：不控制FWER，当组数增加时容易产生假阳性

2. Bonferroni校正

$$\alpha_{adjusted} = \frac{\alpha}{k}$$

优点：严格控制FWER 缺点：过于保守，检验力较低

3. Tukey's HSD (Honestly Significant Difference)

$$HSD = q_{\alpha, k, df_{error}} \times \sqrt{\frac{MS_{error}}{n}}$$

优点：在所有两两比较中控制FWER，较为均衡 适用边界：此处标准HSD公式适用于等样本量；不等样本量可使用Tukey-Kramer扩展，但仍需考虑独立性、正态误差与方差齐性。

4. Scheffe方法

$$\text{临界值} = \sqrt{(k-1) \times F_{\alpha, k-1, df_{error}}}$$

优点：可用于任何对比，不仅限于两两比较 缺点：为所有线性对比提供同时保护，在仅关心两两比较时通常较保守；检验力高低取决于比较范围和数据条件。

💻 Python实现与应用

基础多重比较实现

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.multicomp import pairwise_tukeyhsd, MultiComparison
from statsmodels.stats.libqsturng import psturng, qsturng
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体和绘图样式
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
sns.set_style("whitegrid")

def demonstrate_multiple_comparisons():
    """
    多重比较演示：药物疗效比较
    """
    print("=== 多重比较演示：四种减肥药物疗效比较 ===\n")
    
 # 模拟数据：四种减肥药物的减重效果（kg）
    np.random.seed(42)
    
 # 药物A：轻微效果
    drug_a = np.random.normal(2.5, 1.2, 20)
    
 # 药物B：中等效果
    drug_b = np.random.normal(4.8, 1.5, 18)
    
 # 药物C：显著效果
    drug_c = np.random.normal(7.2, 1.8, 22)
    
 # 安慰剂：无效果
    placebo = np.random.normal(0.8, 1.0, 19)
    
    groups = [placebo, drug_a, drug_b, drug_c]
    group_names = ['安慰剂', '药物A', '药物B', '药物C']
    
    print("各组减重效果统计：")
    for name, group in zip(group_names, groups):
        print(f"{name}: 均值={np.mean(group):.2f}kg, 标准差={np.std(group, ddof=1):.2f}kg, n={len(group)}")
    
 # 首先进行ANOVA检验
    f_stat, p_anova = stats.f_oneway(*groups)
    print(f"\nANOVA结果:")
    print(f"F统计量: {f_stat:.4f}")
    print(f"P值: {p_anova:.4f}")
    
    if p_anova < 0.05:
        print("✅ ANOVA显著，各组均值存在差异，需要进行多重比较")
    else:
        print("❌ ANOVA不显著，各组均值无显著差异")
        return None
    
    return groups, group_names

# 执行演示
demo_groups, demo_names = demonstrate_multiple_comparisons()

LSD最小显著差异法

def lsd_analysis(groups, group_names, alpha=0.05):
    """
    LSD最小显著差异法分析
    """
    print("\n=== LSD最小显著差异法 ===")
    
 # 计算ANOVA相关统计量
    all_data = np.concatenate(groups)
    k = len(groups)  # 组数
    n_total = len(all_data)
    
 # 计算组内均方误差(MSE)
    group_means = [np.mean(group) for group in groups]
    grand_mean = np.mean(all_data)
    
 # 组内平方和
    ss_within = sum(sum((x - group_means[i])**2 for x in groups[i]) for i in range(k))
    df_within = n_total - k
    mse = ss_within / df_within
    
    print(f"组内均方误差(MSE): {mse:.4f}")
    print(f"自由度: {df_within}")
    
 # t临界值
    t_critical = stats.t.ppf(1 - alpha/2, df_within)
    print(f"t临界值(α={alpha}): {t_critical:.4f}")
    
    print("\n=== LSD两两比较结果 ===")
    print(f"{'比较组':<20} {'均值差':<10} {'标准误':<10} {'t值':<10} {'p值':<10} {'显著性':<10}")
    print("-" * 80)
    
 # 进行所有两两比较
    comparisons = []
    for i in range(k):
        for j in range(i+1, k):
 # 计算均值差
            mean_diff = group_means[i] - group_means[j]
            
 # 计算标准误
            n_i, n_j = len(groups[i]), len(groups[j])
            se_diff = np.sqrt(mse * (1/n_i + 1/n_j))
            
 # 计算t统计量
            t_stat = abs(mean_diff) / se_diff
            
 # 计算p值
            p_value = 2 * (1 - stats.t.cdf(t_stat, df_within))
            
 # 判断显著性
            is_significant = p_value < alpha
            significance = "是" if is_significant else "否"
            
            comparison_name = f"{group_names[i]} vs {group_names[j]}"
            
            print(f"{comparison_name:<20} {mean_diff:<10.3f} {se_diff:<10.3f} {t_stat:<10.3f} {p_value:<10.4f} {significance:<10}")
            
            comparisons.append({
                'groups': (group_names[i], group_names[j]),
                'mean_diff': mean_diff,
                'se_diff': se_diff,
                't_stat': t_stat,
                'p_value': p_value,
                'significant': is_significant
            })
    
    print(f"\n⚠️  注意：LSD方法未校正多重比较误差，实际显著性水平可能高于{alpha}")
    
    return comparisons

# 执行LSD分析
lsd_results = lsd_analysis(demo_groups, demo_names)

Bonferroni校正

def bonferroni_analysis(groups, group_names, alpha=0.05):
    """
    Bonferroni校正多重比较
    """
    print("\n=== Bonferroni校正方法 ===")
    
    k = len(groups)
    n_comparisons = k * (k - 1) // 2  # 两两比较次数
    alpha_adjusted = alpha / n_comparisons
    
    print(f"总比较次数: {n_comparisons}")
    print(f"校正后α水平: {alpha_adjusted:.6f}")
    
    print("\n=== Bonferroni校正比较结果 ===")
    print(f"{'比较组':<20} {'均值差':<10} {'p值':<12} {'校正p值':<12} {'显著性':<10}")
    print("-" * 85)
    
    bonferroni_results = []
    significant_pairs = []
    
    for i in range(k):
        for j in range(i+1, k):
 # 进行t检验
            t_stat, p_value = stats.ttest_ind(groups[i], groups[j])
            
 # Bonferroni校正
            p_adjusted = min(p_value * n_comparisons, 1.0)  # 确保不超过1
            
 # 判断显著性
            is_significant = p_adjusted < alpha
            significance = "是" if is_significant else "否"
            
            comparison_name = f"{group_names[i]} vs {group_names[j]}"
            mean_diff = np.mean(groups[i]) - np.mean(groups[j])
            
            print(f"{comparison_name:<20} {mean_diff:<10.3f} {p_value:<12.6f} {p_adjusted:<12.6f} {significance:<10}")
            
            if is_significant:
                significant_pairs.append((group_names[i], group_names[j]))
            
            bonferroni_results.append({
                'groups': (group_names[i], group_names[j]),
                'mean_diff': mean_diff,
                'p_value': p_value,
                'p_adjusted': p_adjusted,
                'significant': is_significant
            })
    
    print(f"\n✅ Bonferroni校正后显著差异的组对:")
    if significant_pairs:
        for pair in significant_pairs:
            print(f"   • {pair[0]} vs {pair[1]}")
    else:
        print("   • 无显著差异的组对")
    
    return bonferroni_results

# 执行Bonferroni分析
bonferroni_results = bonferroni_analysis(demo_groups, demo_names)

Tukey's HSD检验

def tukey_hsd_analysis(groups, group_names, alpha=0.05):
    """
    Tukey's HSD检验
    """
    print("\n=== Tukey's HSD检验 ===")
    
 # 准备数据
    all_data = []
    all_labels = []
    for i, group in enumerate(groups):
        all_data.extend(group)
        all_labels.extend([group_names[i]] * len(group))
    
 # 使用statsmodels进行Tukey检验
    mc = MultiComparison(all_data, all_labels)
    tukey_result = mc.tukeyhsd(alpha=alpha)
    
    print("Tukey HSD结果摘要:")
    print(tukey_result)
    
 # 手动计算验证（适用于等样本量情况）
    print(f"\n=== 手动计算验证 ===")
    
 # 检查样本量是否相等
    sample_sizes = [len(group) for group in groups]
    if len(set(sample_sizes)) == 1:
        print("✅ 各组样本量相等，可使用标准Tukey方法")
        n = sample_sizes[0]
        
 # 计算MSE
        k = len(groups)
        all_data_array = np.concatenate(groups)
        group_means = [np.mean(group) for group in groups]
        
        ss_within = sum(sum((x - group_means[i])**2 for x in groups[i]) for i in range(k))
        df_within = len(all_data_array) - k
        mse = ss_within / df_within
        
 # 计算HSD临界值
 # 使用Studentized Range分布
        q_critical = qsturng(1-alpha, k, df_within)
        hsd = q_critical * np.sqrt(mse / n)
        
        print(f"MSE: {mse:.4f}")
        print(f"q临界值: {q_critical:.4f}")
        print(f"HSD临界值: {hsd:.4f}")
        
        print(f"\n{'比较组':<20} {'均值差':<10} {'|均值差|':<12} {'显著性':<10}")
        print("-" * 65)
        
        for i in range(k):
            for j in range(i+1, k):
                mean_diff = group_means[i] - group_means[j]
                abs_mean_diff = abs(mean_diff)
                is_significant = abs_mean_diff > hsd
                significance = "是" if is_significant else "否"
                
                comparison_name = f"{group_names[i]} vs {group_names[j]}"
                print(f"{comparison_name:<20} {mean_diff:<10.3f} {abs_mean_diff:<12.3f} {significance:<10}")
    
    else:
        print("⚠️  各组样本量不等，使用Tukey-Kramer方法（由statsmodels自动处理）")
    
    return tukey_result

# 执行Tukey分析
tukey_results = tukey_hsd_analysis(demo_groups, demo_names)

多种方法结果比较

def compare_multiple_methods(groups, group_names, alpha=0.05):
    """
    比较不同多重比较方法的结果
    """
    print("\n=== 多种多重比较方法结果对比 ===")
    
    k = len(groups)
    comparisons = []
    
 # 生成所有组对
    for i in range(k):
        for j in range(i+1, k):
            comparisons.append((group_names[i], group_names[j], i, j))
    
 # 存储结果
    results_df = pd.DataFrame()
    comparison_names = [f"{comp[0]} vs {comp[1]}" for comp in comparisons]
    results_df['比较组'] = comparison_names
    
 # 1. 未校正t检验（参考基准）
    uncorrected_p = []
    for comp in comparisons:
        _, p_val = stats.ttest_ind(groups[comp[2]], groups[comp[3]])
        uncorrected_p.append(p_val)
    results_df['未校正p值'] = uncorrected_p
    results_df['未校正显著'] = [p < alpha for p in uncorrected_p]
    
 # 2. Bonferroni校正
    n_comparisons = len(comparisons)
    bonferroni_p = [min(p * n_comparisons, 1.0) for p in uncorrected_p]
    results_df['Bonferroni_p'] = bonferroni_p
    results_df['Bonferroni显著'] = [p < alpha for p in bonferroni_p]
    
 # 3. Holm-Bonferroni校正（逐步方法）
    sorted_indices = sorted(range(len(uncorrected_p)), key=lambda x: uncorrected_p[x])
    holm_significant = [False] * len(uncorrected_p)
    
    for idx, sort_idx in enumerate(sorted_indices):
        correction_factor = n_comparisons - idx
        corrected_p = uncorrected_p[sort_idx] * correction_factor
        
        if corrected_p < alpha:
            holm_significant[sort_idx] = True
        else:
            break  # 一旦有不显著，后续都不显著
    
    results_df['Holm显著'] = holm_significant
    
 # 4. FDR (Benjamini-Hochberg)
    sorted_p = sorted(uncorrected_p)
    fdr_significant = [False] * len(uncorrected_p)
    
    for i in range(len(sorted_p) - 1, -1, -1):
        original_idx = sorted_indices[i]
        critical_value = (i + 1) / len(sorted_p) * alpha
        
        if sorted_p[i] <= critical_value:
 # 这个和所有更小的p值都显著
            for j in range(i + 1):
                orig_idx = sorted_indices[j]
                fdr_significant[orig_idx] = True
            break
    
    results_df['FDR显著'] = fdr_significant
    
    print("方法比较结果:")
    print(results_df.to_string(index=False))
    
 # 统计各方法检出的显著比较数
    print(f"\n=== 各方法检出显著差异数量 ===")
    print(f"未校正t检验: {sum(results_df['未校正显著'])}/{len(comparisons)}")
    print(f"Bonferroni校正: {sum(results_df['Bonferroni显著'])}/{len(comparisons)}")
    print(f"Holm-Bonferroni: {sum(results_df['Holm显著'])}/{len(comparisons)}")
    print(f"FDR校正: {sum(results_df['FDR显著'])}/{len(comparisons)}")
    
    print(f"\n=== 方法特点总结 ===")
    print("• 未校正：检验力最高，但FWER不受控制")
    print("• Bonferroni：最保守，严格控制FWER") 
    print("• Holm-Bonferroni：比Bonferroni略宽松，仍控制FWER")
    print("• FDR：控制假发现率，在多重检验中平衡发现能力和错误控制")
    
    return results_df

# 执行方法比较
method_comparison = compare_multiple_methods(demo_groups, demo_names)

🏢 实际应用场景

1. 医学研究：药物剂量效应分析

def medical_dose_response_analysis():
    """
    医学研究案例：不同剂量药物的疗效比较
    """
    print("=== 医学案例：降血压药物剂量效应研究 ===")
    print("研究目标：比较四种不同剂量对收缩压降低效果的差异\n")
    
 # 模拟数据：收缩压降低量（mmHg）
    np.random.seed(123)
    
 # 安慰剂组
    placebo = np.random.normal(2, 4, 25)
    placebo = np.clip(placebo, -5, 15)  # 限制在合理范围
    
 # 低剂量组
    low_dose = np.random.normal(8, 5, 24)
    low_dose = np.clip(low_dose, 0, 25)
    
 # 中剂量组
    medium_dose = np.random.normal(15, 6, 26)
    medium_dose = np.clip(medium_dose, 5, 35)
    
 # 高剂量组
    high_dose = np.random.normal(22, 7, 23)
    high_dose = np.clip(high_dose, 10, 45)
    
    dose_groups = [placebo, low_dose, medium_dose, high_dose]
    dose_names = ['安慰剂', '低剂量', '中剂量', '高剂量']
    
    print("各组降压效果统计：")
    for name, group in zip(dose_names, dose_groups):
        print(f"{name}: 平均降压{np.mean(group):.1f}mmHg, 标准差{np.std(group, ddof=1):.1f}mmHg, n={len(group)}")
    
 # ANOVA检验
    f_stat, p_anova = stats.f_oneway(*dose_groups)
    print(f"\nANOVA结果: F={f_stat:.3f}, p={p_anova:.6f}")
    
    if p_anova < 0.05:
        print("✅ 存在显著差异，进行多重比较")
        
 # 使用Tukey检验
        all_data = []
        all_labels = []
        for i, group in enumerate(dose_groups):
            all_data.extend(group)
            all_labels.extend([dose_names[i]] * len(group))
        
        mc = MultiComparison(all_data, all_labels)
        tukey_result = mc.tukeyhsd(alpha=0.05)
        
        print("\nTukey HSD结果:")
        print(tukey_result)
        
 # 临床意义解释
        print("\n=== 临床意义解释 ===")
        group_means = [np.mean(group) for group in dose_groups]
        
        print("💊 剂量-效应关系分析:")
        for i, (name, mean) in enumerate(zip(dose_names, group_means)):
            if i == 0:
                print(f"• {name}: 基线效应 {mean:.1f}mmHg")
            else:
                net_effect = mean - group_means[0]
                print(f"• {name}: 净效应 {net_effect:.1f}mmHg")
                
 # 教学假设：以净降压>10mmHg、5-10mmHg和<5mmHg演示分层输出；这些阈值不是本文已建立的临床判断标准
                if net_effect > 10:
                    print(f"  示例阈值分层：净降压>10mmHg，不代表已证实临床意义")
                elif net_effect > 5:
                    print(f"  示例阈值分层：净降压5-10mmHg，不代表已证实临床意义")
                else:
                    print(f"  示例阈值分层：净降压<5mmHg，不代表临床意义判断")
        
 # 最佳剂量推荐
        best_dose_idx = np.argmax(group_means)
        print(f"\n💡 模拟数据中平均降压最大的组别: {dose_names[best_dose_idx]}（平均降压{group_means[best_dose_idx]:.1f}mmHg）；未评估安全性，不构成剂量推荐")
        
 # 成本效益分析
        costs = [0, 100, 180, 250]  # 每月药物成本（元）
        print("\n💰 成本效益分析:")
        for i, (name, cost) in enumerate(zip(dose_names, costs)):
            if i > 0:
                net_effect = group_means[i] - group_means[0]
                if net_effect > 0:
                    cost_per_mmhg = cost / net_effect
                    print(f"• {name}: {cost}元/月, 成本效益比 {cost_per_mmhg:.1f}元/mmHg")
        
    return dose_groups, dose_names

# 执行医学案例分析
medical_groups, medical_names = medical_dose_response_analysis()

2. 市场研究：广告策略效果比较

def marketing_strategy_analysis():
    """
    市场研究案例：不同广告策略的转化率比较
    """
    print("\n=== 市场案例：广告策略效果研究 ===")
    print("研究目标：比较五种广告策略对产品转化率的影响\n")
    
 # 模拟数据：转化率（%）
    np.random.seed(456)
    
    strategies = {
        '传统电视广告': np.random.normal(2.3, 0.8, 30),
        '社交媒体广告': np.random.normal(4.1, 1.2, 32),  
        '搜索引擎广告': np.random.normal(5.8, 1.5, 28),
        '影响者营销': np.random.normal(6.9, 1.8, 25),
        '内容营销': np.random.normal(3.7, 1.1, 27)
    }
    
 # 确保转化率为正值
    for key in strategies:
        strategies[key] = np.clip(strategies[key], 0.1, None)
    
    strategy_groups = list(strategies.values())
    strategy_names = list(strategies.keys())
    
    print("各策略转化率统计：")
    for name, group in zip(strategy_names, strategy_groups):
        print(f"{name}: 平均{np.mean(group):.2f}%, 标准差{np.std(group, ddof=1):.2f}%, 样本量{len(group)}")
    
 # ANOVA检验
    f_stat, p_anova = stats.f_oneway(*strategy_groups)
    print(f"\nANOVA结果: F={f_stat:.3f}, p={p_anova:.6f}")
    
    if p_anova < 0.05:
        print("✅ 策略间存在显著差异")
        
 # 多重比较 - 使用Welch t检验配合Bonferroni校正，不是Games-Howell检验
        print("\n=== 多重比较结果 ===")
        
 # 由于样本量不等且可能方差不齐，使用Welch t检验配合Bonferroni校正
        from itertools import combinations
        
        n_comparisons = len(list(combinations(range(len(strategy_groups)), 2)))
        alpha_bonf = 0.05 / n_comparisons
        
        print(f"Bonferroni校正α水平: {alpha_bonf:.4f}")
        print(f"{'策略对比':<30} {'均值差':<10} {'p值':<12} {'显著性':<10}")
        print("-" * 65)
        
        significant_pairs = []
        all_results = []
        
        for i, j in combinations(range(len(strategy_groups)), 2):
 # Welch t检验（不假设方差相等）
            t_stat, p_val = stats.ttest_ind(strategy_groups[i], strategy_groups[j], equal_var=False)
            mean_diff = np.mean(strategy_groups[i]) - np.mean(strategy_groups[j])
            
            is_significant = p_val < alpha_bonf
            significance = "是" if is_significant else "否"
            
            comparison_name = f"{strategy_names[i]} vs {strategy_names[j]}"
            print(f"{comparison_name:<30} {mean_diff:<10.3f} {p_val:<12.6f} {significance:<10}")
            
            if is_significant:
                significant_pairs.append((strategy_names[i], strategy_names[j], mean_diff))
            
            all_results.append({
                'comparison': comparison_name,
                'mean_diff': mean_diff,
                'p_value': p_val,
                'significant': is_significant
            })
        
 # 策略排序和推荐
        print(f"\n=== 策略效果排序 ===")
        strategy_means = [(name, np.mean(group)) for name, group in zip(strategy_names, strategy_groups)]
        strategy_means.sort(key=lambda x: x[1], reverse=True)
        
        print("🏆 转化率排序:")
        for i, (name, mean) in enumerate(strategy_means, 1):
            print(f"{i}. {name}: {mean:.2f}%")
        
 # ROI分析
        print(f"\n💰 基于假设成本的成本效益演示（非ROI）:")
        costs = [50000, 30000, 40000, 80000, 25000]  # 每月广告成本
        cost_dict = dict(zip(strategy_names, costs))
        
        print(f"{'策略':<15} {'转化率':<10} {'月成本':<10} {'成本效益':<15}")
        print("-" * 55)
        
        for name, mean in strategy_means:
            cost = cost_dict[name]
            cost_effectiveness = cost / mean  # 每1%转化率的成本
            print(f"{name:<15} {mean:<10.2f} {cost:<10} {cost_effectiveness:<15.0f}")
        
 # 推荐最佳组合
        print(f"\n💡 策略建议:")
        best_performance = strategy_means[0]
        best_cost_effective = min(strategy_means, key=lambda x: cost_dict[x[0]]/x[1])
        
        print(f"• 最佳效果: {best_performance[0]} (转化率{best_performance[1]:.2f}%)")
        print(f"• 最佳性价比: {best_cost_effective[0]} (每%转化成本{cost_dict[best_cost_effective[0]]/best_cost_effective[1]:.0f}元)")
        
    return strategy_groups, strategy_names

# 执行市场案例分析
marketing_groups, marketing_names = marketing_strategy_analysis()

3. 教育研究：教学方法效果评估

def education_method_analysis():
    """
    教育研究案例：不同教学方法对学习效果的影响
    """
    print("\n=== 教育案例：编程教学方法效果研究 ===")
    print("研究目标：比较四种编程教学方法对学生成绩的影响\n")
    
 # 模拟数据：期末考试成绩（0-100分）
    np.random.seed(789)
    
    methods = {
        '传统讲授': np.random.normal(72, 12, 35),
        '项目驱动': np.random.normal(81, 10, 33),
        '同伴编程': np.random.normal(78, 11, 36),
        '翻转课堂': np.random.normal(85, 9, 34)
    }
    
 # 确保成绩在合理范围内
    for key in methods:
        methods[key] = np.clip(methods[key], 0, 100)
    
    method_groups = list(methods.values())
    method_names = list(methods.keys())
    
    print("各教学方法成绩统计：")
    for name, group in zip(method_names, method_groups):
        print(f"{name}: 平均{np.mean(group):.1f}分, 标准差{np.std(group, ddof=1):.1f}分, 学生数{len(group)}")
    
 # 效应量计算（eta squared）
    all_scores = np.concatenate(method_groups)
    grand_mean = np.mean(all_scores)
    
    ss_total = np.sum((all_scores - grand_mean)**2)
    ss_between = sum(len(group) * (np.mean(group) - grand_mean)**2 for group in method_groups)
    eta_squared = ss_between / ss_total
    
    print(f"\n效应量(η²): {eta_squared:.3f}")
    if eta_squared > 0.14:
        effect_size = "大"
    elif eta_squared > 0.06:
        effect_size = "中等"  
    elif eta_squared > 0.01:
        effect_size = "小"
    else:
        effect_size = "极小"
    print(f"效应量大小: {effect_size}")
    
 # ANOVA检验
    f_stat, p_anova = stats.f_oneway(*method_groups)
    print(f"\nANOVA结果: F={f_stat:.3f}, p={p_anova:.6f}")
    
    if p_anova < 0.05:
        print("✅ 教学方法间存在显著差异")
        
 # Tukey HSD检验
        all_data = []
        all_labels = []
        for i, group in enumerate(method_groups):
            all_data.extend(group)
            all_labels.extend([method_names[i]] * len(group))
        
        mc = MultiComparison(all_data, all_labels)
        tukey_result = mc.tukeyhsd(alpha=0.05)
        
        print("\nTukey HSD多重比较结果:")
        
 # 提取显著的配对
        tukey_df = pd.DataFrame(data=tukey_result._results_table.data[1:], 
                               columns=tukey_result._results_table.data[0])
        
        significant_comparisons = tukey_df[tukey_df['reject'] == True]
        
        if len(significant_comparisons) > 0:
            print("\n显著差异的教学方法配对:")
            for idx, row in significant_comparisons.iterrows():
                group1, group2 = row['group1'], row['group2']
                mean_diff = row['meandiff']
                p_adj = row['p-adj']
                print(f"• {group1} vs {group2}: 分数差异 {mean_diff:.1f}分 (p={p_adj:.4f})")
        
 # 教学效果分析
        print(f"\n=== 教学效果分析 ===")
        method_means = [(name, np.mean(group)) for name, group in zip(method_names, method_groups)]
        method_means.sort(key=lambda x: x[1], reverse=True)
        
        print("🎓 教学方法效果排序:")
        baseline_score = min(mean for _, mean in method_means)
        
        for i, (name, mean) in enumerate(method_means, 1):
            improvement = mean - baseline_score
            print(f"{i}. {name}: {mean:.1f}分 (相对提升 {improvement:.1f}分)")
        
 # 实施建议
        print(f"\n📚 教学实施建议:")
        best_method = method_means[0]
        print(f"• 模拟数据中平均分最高的是'{best_method[0]}'教学法（平均分{best_method[1]:.1f}分），不构成真实教学效果结论")
        
 # 计算需要的样本量来检测差异
        from scipy.stats import power
        print(f"\n📊 统计功效分析:")
        print(f"• 当前研究的效应量为{effect_size}({eta_squared:.3f})")
        print(f"• 每组30名学生仅可作为示例规划假设，不能确保统计功效；应根据目标效应、变异、显著性水平、多重比较方案和目标功效计算样本量")
        
    return method_groups, method_names

# 执行教育案例分析
education_groups, education_names = education_method_analysis()

📊 多重比较可视化

def visualize_multiple_comparisons(groups, group_names, title="多重比较结果可视化"):
    """
    多重比较结果综合可视化
    """
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    fig.suptitle(f'{title}', fontsize=16, fontweight='bold')
    
 # 准备数据
    all_data = []
    all_labels = []
    for i, group in enumerate(groups):
        all_data.extend(group)
        all_labels.extend([group_names[i]] * len(group))
    
    df_viz = pd.DataFrame({'value': all_data, 'group': all_labels})
    
 # 1. 箱线图对比
    sns.boxplot(data=df_viz, x='group', y='value', ax=axes[0, 0])
    axes[0, 0].set_title('各组数据分布对比')
    axes[0, 0].tick_params(axis='x', rotation=45)
    axes[0, 0].grid(True, alpha=0.3)
    
 # 2. 均值置信区间图
    group_stats = df_viz.groupby('group')['value'].agg(['mean', 'std', 'count']).reindex(group_names).reset_index()
    group_stats['se'] = group_stats['std'] / np.sqrt(group_stats['count'])
    group_stats['ci'] = stats.t.ppf(0.975, group_stats['count'] - 1) * group_stats['se']  # 正态独立样本下各组均值的逐组95%区间，未经多重比较校正
    
    x_pos = range(len(group_names))
    axes[0, 1].errorbar(x_pos, group_stats['mean'], yerr=group_stats['ci'], 
                       fmt='o', capsize=5, capthick=2, markersize=8)
    axes[0, 1].set_xticks(x_pos)
    axes[0, 1].set_xticklabels(group_names, rotation=45)
    axes[0, 1].set_title('均值及逐组95%置信区间（未校正多重比较）')
    axes[0, 1].grid(True, alpha=0.3)
    
 # 3. Tukey HSD结果可视化
    mc = MultiComparison(all_data, all_labels)
    tukey_result = mc.tukeyhsd(alpha=0.05)
    
 # 创建显著性矩阵
    n_groups = len(group_names)
    sig_matrix = np.zeros((n_groups, n_groups))
    
    tukey_df = pd.DataFrame(data=tukey_result._results_table.data[1:], 
                           columns=tukey_result._results_table.data[0])
    
    for idx, row in tukey_df.iterrows():
        group1_idx = group_names.index(row['group1'])
        group2_idx = group_names.index(row['group2'])
        if row['reject']:
            sig_matrix[group1_idx, group2_idx] = 1
            sig_matrix[group2_idx, group1_idx] = 1
    
 # 绘制显著性热图
    sns.heatmap(sig_matrix, annot=True, cmap='RdYlBu_r', 
                xticklabels=group_names, yticklabels=group_names,
                ax=axes[0, 2], cbar_kws={'label': '显著差异'})
    axes[0, 2].set_title('Tukey HSD显著性矩阵')
    
 # 4. 效应量可视化
    group_means = [np.mean(group) for group in groups]
    grand_mean = np.mean(all_data)
    effect_sizes = [(mean - grand_mean) for mean in group_means]
    
    colors = ['red' if es < 0 else 'blue' for es in effect_sizes]
    bars = axes[1, 0].bar(group_names, effect_sizes, color=colors, alpha=0.7)
    axes[1, 0].axhline(y=0, color='black', linestyle='-', alpha=0.8)
    axes[1, 0].set_title('各组相对于总均值的效应')
    axes[1, 0].tick_params(axis='x', rotation=45)
    
 # 添加数值标签
    for bar, es in zip(bars, effect_sizes):
        height = bar.get_height()
        axes[1, 0].text(bar.get_x() + bar.get_width()/2., height + (0.1 if height >= 0 else -0.3),
                        f'{es:.2f}', ha='center', va='bottom' if height >= 0 else 'top')
    
 # 5. p值分布图
    p_values = []
    comparison_labels = []
    
    for i in range(len(groups)):
        for j in range(i+1, len(groups)):
            _, p_val = stats.ttest_ind(groups[i], groups[j])
            p_values.append(p_val)
            comparison_labels.append(f"{group_names[i][:3]}-{group_names[j][:3]}")
    
    axes[1, 1].scatter(range(len(p_values)), p_values, s=50, alpha=0.7)
    axes[1, 1].axhline(y=0.05, color='red', linestyle='--', label='α=0.05')
    axes[1, 1].axhline(y=0.05/len(p_values), color='orange', linestyle='--', 
                      label=f'Bonferroni α={0.05/len(p_values):.3f}')
    axes[1, 1].set_xticks(range(len(p_values)))
    axes[1, 1].set_xticklabels(comparison_labels, rotation=45)
    axes[1, 1].set_ylabel('p值')
    axes[1, 1].set_title('各组对比p值分布')
    axes[1, 1].legend()
    axes[1, 1].set_yscale('log')
    
 # 6. 统计摘要
    axes[1, 2].axis('off')
    
 # ANOVA结果
    f_stat, p_anova = stats.f_oneway(*groups)
    
 # 效应量
    ss_total = np.sum((all_data - grand_mean)**2)
    ss_between = sum(len(group) * (np.mean(group) - grand_mean)**2 for group in groups)
    eta_squared = ss_between / ss_total
    
 # 显著比较计数
    n_significant_tukey = sum(tukey_df['reject'])
    n_total_comparisons = len(tukey_df)
    
    summary_text = f"""
    多重比较分析摘要:
    
    ANOVA结果:
    • F统计量: {f_stat:.3f}
    • p值: {p_anova:.6f}
    • 效应量(η²): {eta_squared:.3f}
    
    Tukey HSD结果:
    • 总比较数: {n_total_comparisons}
    • 显著差异数: {n_significant_tukey}
    • 显著比例: {n_significant_tukey/n_total_comparisons*100:.1f}%
    
    组间差异:
    • 最大均值: {max(group_means):.2f}
    • 最小均值: {min(group_means):.2f}
    • 极差: {max(group_means)-min(group_means):.2f}
    
    描述与方法:
    • 样本均值最高组别: {group_names[np.argmax(group_means)]}（不代表最优选择）
    • 当前比较方法: Tukey HSD（适用性需结合独立性、正态误差与方差齐性判断）
    • Tukey同时推断的名义置信水平: 95%
    """
    
    axes[1, 2].text(0.1, 0.9, summary_text, fontsize=10, verticalalignment='top',
                     bbox=dict(boxstyle="round,pad=0.3", facecolor="lightgray", alpha=0.8))
    
    plt.tight_layout()
    plt.show()

# 对演示数据进行可视化
if demo_groups is not None:
    visualize_multiple_comparisons(demo_groups, demo_names, "药物疗效多重比较分析")

⚠️ 注意事项与最佳实践

📚 扩展学习资源

理论深化

• 📖 《Multiple Comparisons and Multiple Tests》 - Hochberg & Tamhane
• 📖 《统计学中的多重比较》 - 深入理解各种校正方法
• 🎓 在线课程：Coursera《实验设计与分析》

实践工具

• 🐍 Python库：SciPy、Statsmodels、Pingouin、Scikit-posthocs
• 📊 专业软件：R (multcomp包)、SPSS、SAS
• 🔧 在线工具：GraphPad Prism、JASP

进阶主题

• 方差分析（ANOVA） - 组均值多重比较的常见背景，并非所有多重比较的必要前提
• 实验设计 - 如何设计避免多重比较问题
• 假发现率控制 - FDR方法详解

#统计学 #多重比较 #事后检验 #ANOVA #假设检验 #实验设计 #Python统计分析