基于 Python 的多重比较分析流程

书接上回，现在来介绍使用 Python 的 scikit_posthocs 库，使用 Tukey_HSD 进行多重比较。为了丰富内容，把基本的一个分析流程也一起过一遍。

导入库与生成数据

# 导入库import numpy as npimport pandas as pdimport seaborn as snsimport matplotlib.pyplotas plt# 正态性检验from scipy.statsimport shapiro# 多重比较Tukey HSDimport scikit_posthocs as sp%matplotlib inline# 设置绘图样式为白底sns.set_theme(style='ticks')# 图中文为宋体，英文为Times NewRoman，解决负号显示问题plt.rcParams['font.sans-serif']=['SimSun','Times New Roman']  # 设置中文字体为宋体plt.rcParams['axes.unicode_minus']=False    # 解决负号显示问题# 生成数据df = pd.DataFrame({'group':['T1','T1','T1','T2','T2','T2','T3','T3','T3','T4','T4','T4'],'value':[84,82,83,83,85,87,87,88,89,90,91,92]})df

输出结果：

正态分布检验

这边采用较为常用的夏皮洛-威尔克检验（Shapiro—Wilk test），我先前的推文也有介绍过。

# W检验# 夏皮洛-威尔克检验（Shapiro—Wilk test），一般又称W检验。from scipy.statsimport shapironormality =shapiro(df.value)normality# 显示结果print("Shapiro-Wilk检验结果：")print(f"统计量: {normality.statistic}")print(f"p值: {normality.pvalue}")# 判断正态性alpha =0.05if normality.pvalue> alpha:print(f"\n结论：因为p值 = {normality.pvalue:.4f} > {alpha}，所以数据服从正态分布，可以进行后续的方差分析和多重比较分析")else:print(f"\n结论：因为p值 = {normality.pvalue:.4f} ≤ {alpha}，所以数据不服从正态分布，建议进行非参数检验或数据转换后再进行分析")

输出结果：

Shapiro-Wilk 检验结果：统计量: 0.9471656341738938p 值: 0.5959878635330442

结论: 因为 p 值 = 0.5960 > 0.05，所以数据服从正态分布，可以进行后续的方差分析和多重比较分析

方差齐性检验

方差齐性推荐使用 Levene 检验

# 方差齐性检验，推荐Levene检验from scipy.statsimport levene# 提取各组的数值数据groups_data =[df[df['group']== group]['value'].valuesfor group in df['group'].unique()]# 执行Levene检验levene_stat, levene_p =levene(*groups_data)print(f"Levene检验统计量: {levene_stat:.4f}")print(f"p值: {levene_p:.4f}")# 判断方差齐性alpha =0.05if levene_p > alpha:print(f"\n结论：因为p值 = {levene_p:.4f} > {alpha}，所以数据满足方差齐性假设")else:print(f"\n结论：因为p值 = {levene_p:.4f} ≤ {alpha}，所以数据不满足方差齐性假设")

输出结果：

Levene 检验统计量: 0.5714p 值: 0.6495

结论：因为 p 值 = 0.6495 > 0.05，所以数据满足方差齐性假设

多重比较

这边输出结果不展示了，可以自行运行查看。

# 进行多重比较,基本结果tukey_result = sp.posthoc_tukey(df, val_col='value', group_col='group')print("\nTukey HSD多重比较结果：")print(tukey_result)# 这边生成的为p值，可自行两两比较，判断是否显著不同# 获取显著性字母标记cld_result = sp.compact_letter_display(tukey_result)print("\n多重比较字母标记法结果：")cld_result# 生成数据表格summary_stats = df.groupby('group')['value'].agg(Mean='mean',     # 均值Std='std'        # 标准差).reset_index()     # 将分组列变回普通列，方便后续合并# 将 summary_stats 按 'Mean' 降序排序，并以此顺序重新生成字母标记summary_stats_sorted = summary_stats.sort_values('Mean', ascending=False)groups_sorted = summary_stats_sorted['group'].tolist()p_values_sorted = tukey_result.loc[groups_sorted, groups_sorted]cld_result_sorted = sp.compact_letter_display(p_values_sorted)# 将字母标记添加到统计摘要表中summary_stats['Sig']= summary_stats['group'].map(cld_result_sorted)# 为确保表格清晰，最终可按均值降序排列（让最重要的结果一目了然）final_table = summary_stats.sort_values('Mean', ascending=False).reset_index(drop=True)# 将小数位数格式化为更美观的格式final_table['Mean']= final_table['Mean'].map('{:.2f}'.format)final_table['Std']= final_table['Std'].map('{:.2f}'.format)final_table = final_table.sort_values("group", ascending=True)# 输出最终结果表print(final_table)# 数据结果变成均值±标准差和显著性标记的形式final_table['Result']= final_table.apply(lambda row: f"{row['Mean']}±{row['Std']} {row['Sig']}", axis=1)print(final_table[['group','Result']])

要注意一点的是，我最后加的 Result 是为了方便放入三线表中，仅供参考。

带误差棒的柱状图绘制