大家好，我是小寒
今天给大家分享统计学中的一个关键概念，方差分析
方差分析（Analysis of Variance，ANOVA）是统计学中最重要的假设检验方法之一，由英国统计学家 Ronald A. Fisher 于 20 世纪提出。
它的核心目的只有一个：判断多个总体均值是否存在显著差异。
例如：
如果只有两组数据，我们通常使用 t 检验；而当组数大于等于 3 时，就需要使用方差分析。
核心原理：变异的分解
方差分析的核心思想是将数据的总变异拆分为两部分：组间变异和组内变异。
- 总变异：所有观测值与其总体均值之间的偏离程度，由因变量自身的波动以及实验因素、随机误差共同导致。
- 组间变异：各个处理组均值与整体均值之间的差异。这种差异既包含了随机误差，也包含了由于实施了不同的处理（如不同的药量、不同的工艺）所带来的处理效应。
- 组内变异：同一组内部个体之间的差异。因为同一组处理条件相同，这种差异纯粹是由随机测量误差、个体生理差异等不可控因素引起的。
如果组间变异显著大于组内变异，说明“处理因素”起了关键作用，各个组的均值不全相等；反之，如果组间差异和组内随机噪声差不多，就无法拒绝“各组均值相等”的零假设。
数学公式
下面我们以单因素方差分析为例进行说明。
假设有 个独立的组（处理水平），第 组有 个观测值，总样本量为 。
1.原假设与备择假设
- **：各组总体均值 不全相等（至少存在一对 使得 ）
2.平方和分解
总平方和
衡量所有数据点偏离总体均值的程度
自由度 。
组间平方和
衡量各组样本均值偏离总体均值的程度（组间变异）
自由度 。
组内平方和
衡量各组内部观测值偏离本组均值的程度（随机误差引起的变异）：
自由度 。
恒等式
自由度同样满足关系：。
3.均方
因为平方和的大小受自由度（样本量或组数）影响，直接比较平方和不公平，需要除以各自的自由度得到均方（均方差）。
4. 统计量
构造统计量 ：
在 为真的条件下，统计量 服从自由度为 的 分布：
5.决策规则
对给定的显著性水平 （如 0.05 或 0.01），计算得到 值或查找临界值 ：
- 若 （或 ），则拒绝零假设 ，认为不同处理水平之间的总体均值存在显著差异。
- 若 （或 ），则无法拒绝零假设 ，没有足够证据表明变量间存在显著效应。
方差分析的类型
单因素方差分析
用于检验一个分类自变量对一个连续因变量是否有显著影响。
例如，比较三种不同教学方法对学生期末考试成绩的影响是否有显著差异。
双因素/多因素方差分析
用于同时检验两个或多个分类自变量各自的主效应，以及它们之间的交互作用对一个连续因变量的影响。
例如，研究 “不同药物剂量（低、中、高）” 与 “患者性别（男、女）” 对降低血压（连续数值）的共同与独立影响。
重复测量方差分析
用于分析同一批被试对象在不同时间点或不同实验条件下重复测量得到的连续数据，能有效控制个体差异导致的噪音。
例如，追踪同一组患者在接受某种新药治疗前、治疗 1 周后、治疗 1 个月后以及治疗 3 个月后的血糖变化情况。
案例分享
某农业研究所想知道三种不同配方的肥料（Fertilizer A, Fertilizer B, Fertilizer C）对某种水稻产量的影响是否存在显著差异。
研究所选择了 块土质、面积、日照等条件基本相同的小型试验田。
将这 30 块田随机分为 3 组，每组 块。
在收获季节，测量并记录每块试验田的水稻产量（单位：千克/公顷）。
import numpy as npimport pandas as pdimport scipy.stats as statsimport statsmodels.api as smfrom statsmodels.formula.api import olsfrom statsmodels.stats.multicomp import pairwise_tukeyhsdimport matplotlib.pyplot as pltimport seaborn as sns# 1. 设定随机种子以确保结果可复现，并生成模拟实验数据np.random.seed(42)yield_A = np.round(np.random.normal(loc=14.2, scale=1.2, size=10), 2)yield_B = np.round(np.random.normal(loc=18.5, scale=1.1, size=10), 2)yield_C = np.round(np.random.normal(loc=16.1, scale=1.3, size=10), 2)df = pd.DataFrame({'Fertilizer': ['Fertilizer_A'] * 10 + ['Fertilizer_B'] * 10 + ['Fertilizer_C'] * 10,'Yield': np.concatenate([yield_A, yield_B, yield_C])})# 2. 描述性统计print("=== 各组描述性统计量 ===")summary = df.groupby('Fertilizer')['Yield'].agg(['count', 'mean', 'std', 'min', 'max'])print(summary)# 3. 前提假设检验# (1) 正态性检验 (Shapiro-Wilk)print("\n=== 正态性检验 (Shapiro-Wilk Test) ===")for fert, group in df.groupby('Fertilizer'):stat, p_val = stats.shapiro(group['Yield'])print(f"{fert}: W = {stat:.4f}, p-value = {p_val:.4f}")# (2) 方差齐性检验 (Levene's Test)print("\n=== 方差齐性检验 (Levene's Test) ===")lev_stat, lev_p = stats.levene(yield_A, yield_B, yield_C)print(f"Levene Statistic = {lev_stat:.4f}, p-value = {lev_p:.4f}")# 4. 单因素方差分析 (One-Way ANOVA)model = ols('Yield ~ C(Fertilizer)', data=df).fit()anova_table = sm.stats.anova_lm(model, typ=1)# 计算效应量 Eta-squared (η²)ss_between = anova_table.loc['C(Fertilizer)', 'sum_sq']ss_total = anova_table['sum_sq'].sum()eta_sq = ss_between / ss_totalprint("\n=== 方差分析表 (ANOVA Table) ===")print(anova_table)print(f"\n效应量 Effect Size (Eta-squared η²) = {eta_sq:.4f}")# 5. 事后检验 (Tukey's HSD Post-Hoc Test)tukey = pairwise_tukeyhsd(endog=df['Yield'], groups=df['Fertilizer'], alpha=0.05)print("\n=== 事后检验 (Tukey HSD Test) ===")print(tukey)# 6. 图表绘制（绘制 2x2 综合分析仪表盘）plt.style.use('seaborn-v0_8-whitegrid'if'seaborn-v0_8-whitegrid'in plt.style.available else'default')fig, axes = plt.subplots(2, 2, figsize=(14, 10))# 图 1：各组产量数据分布（箱线图 + 抖动散点图）sns.boxplot(x='Fertilizer', y='Yield', data=df, ax=axes[0, 0], palette='Set2', width=0.4)sns.stripplot(x='Fertilizer', y='Yield', data=df, ax=axes[0, 0], color='black', alpha=0.6, jitter=0.2, size=7)axes[0, 0].set_title('1. Tomato Crop Yield Distribution by Fertilizer', fontsize=12, fontweight='bold')axes[0, 0].set_xlabel('Fertilizer Type', fontsize=11)axes[0, 0].set_ylabel('Yield (kg / plant)', fontsize=11)axes[0, 0].set_xticklabels(['Fertilizer A\n(Standard NPK)', 'Fertilizer B\n(Bio-Fertilizer)', 'Fertilizer C\n(Organic Compound)'])# 图 2：残差 Q-Q 图（正态性诊断）residuals = model.residstats.probplot(residuals, dist="norm", plot=axes[0, 1])axes[0, 1].set_title('2. Residual Q-Q Plot (Normality Check)', fontsize=12, fontweight='bold')axes[0, 1].get_lines()[0].set_markerfacecolor('crimson')axes[0, 1].get_lines()[0].set_markeredgecolor('crimson')axes[0, 1].get_lines()[1].set_color('black')# 图 3：理论 F 分布密度曲线与拒绝域df_num, df_den = 2, 33f_stat = anova_table.loc['C(Fertilizer)', 'F']f_critical = stats.f.ppf(1 - 0.05, df_num, df_den)x_f = np.linspace(0, 32, 500)y_f = stats.f.pdf(x_f, df_num, df_den)axes[1, 0].plot(x_f, y_f, 'b-', lw=2, label=f'F-distribution ({df_num}, {df_den})')x_shade = np.linspace(f_critical, 32, 200)axes[1, 0].fill_between(x_shade, 0, stats.f.pdf(x_shade, df_num, df_den), color='red', alpha=0.3, label=f'Rejection Region (alpha=0.05, F_crit={f_critical:.2f})')axes[1, 0].axvline(f_stat, color='darkgreen', linestyle='--', linewidth=2, label=f'Observed F = {f_stat:.2f}')axes[1, 0].set_title('3. Theoretical F-Distribution & Critical Region', fontsize=12, fontweight='bold')axes[1, 0].set_xlabel('F value', fontsize=11)axes[1, 0].set_ylabel('Probability Density', fontsize=11)axes[1, 0].legend(loc='upper right')# 图 4：Tukey HSD 事后检验均值差 95% 置信区间tukey_res = tukey.summary()tukey_data = pd.DataFrame(tukey_res.data[1:], columns=tukey_res.data[0])comparison_labels = [f"{row['group1']} vs {row['group2']}"for _, row in tukey_data.iterrows()]means_diff = tukey_data['meandiff'].astype(float)lower_ci = tukey_data['lower'].astype(float)upper_ci = tukey_data['upper'].astype(float)errors = [means_diff - lower_ci, upper_ci - means_diff]y_pos = np.arange(len(comparison_labels))axes[1, 1].errorbar(means_diff, y_pos, xerr=errors, fmt='o', color='purple', ecolor='purple', elinewidth=2, capsize=6, capthick=2)axes[1, 1].axvline(0, color='gray', linestyle='--')axes[1, 1].set_yticks(y_pos)axes[1, 1].set_yticklabels(comparison_labels)axes[1, 1].set_title('4. Tukey HSD Pairwise Differences (95% CI)', fontsize=12, fontweight='bold')axes[1, 1].set_xlabel('Mean Difference in Yield (kg / plant)', fontsize=11)plt.tight_layout()plt.show()
