Python案例如何用Statsmodels做方差分析

wen python案例 2026-07-16 223

本文目录导读：

单因素方差分析（One-way ANOVA）
双因素方差分析（Two-way ANOVA）
重复测量ANOVA（Repeated Measures ANOVA）
方差分析前提假设检验
非参数替代方法（Kruskal-Wallis检验）
可视化结果
完整ANOVA分析流程示例
关键注意事项

我来详细介绍如何使用Statsmodels进行方差分析（ANOVA）。

单因素方差分析（One-way ANOVA）

准备数据和基本ANOVA

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
# 创建示例数据：三种不同肥料对植物生长的影响
np.random.seed(42)
# 生成数据
fertilizer_A = np.random.normal(20, 3, 30)  # 肥料A
fertilizer_B = np.random.normal(25, 3, 30)  # 肥料B
fertilizer_C = np.random.normal(22, 3, 30)  # 肥料C
# 创建DataFrame
df = pd.DataFrame({
    'growth': np.concatenate([fertilizer_A, fertilizer_B, fertilizer_C]),
    'fertilizer': ['A']*30 + ['B']*30 + ['C']*30
})
print("数据预览:")
print(df.head())
print(f"\n各组统计量:")
print(df.groupby('fertilizer')['growth'].describe())

执行单因素ANOVA

# 方法1：使用ols模型
model = ols('growth ~ C(fertilizer)', data=df).fit()
anova_results = anova_lm(model)
print("单因素方差分析结果:")
print(anova_results)
# 方法2：直接从statsmodels获取
from statsmodels.stats.oneway import oneway
# 准备数据
groups = [df[df['fertilizer'] == g]['growth'].values for g in ['A', 'B', 'C']]
# 执行One-way ANOVA
f_stat, p_val, means = oneway(groups, use_var='equal')
print(f"\n直接One-way ANOVA结果:")
print(f"F统计量: {f_stat:.4f}")
print(f"p值: {p_val:.4f}")

事后检验（Post-hoc tests）

from statsmodels.stats.multicomp import pairwise_tukeyhsd
# Tukey HSD事后检验
tukey = pairwise_tukeyhsd(df['growth'], df['fertilizer'], alpha=0.05)
print("\nTukey HSD事后检验结果:")
print(tukey)
# 可视化事后检验结果
tukey.plot_simultaneous()'Tukey HSD检验结果')
plt.show()

双因素方差分析（Two-way ANOVA）

# 创建双因素数据：肥料类型和土壤类型对植物生长的影响
np.random.seed(42)
# 生成数据
data = []
fertilizers = ['A', 'B', 'C']
soils = ['clay', 'sandy', 'loamy']
for fert in fertilizers:
    for soil in soils:
        # 不同组合有不同的均值
        if fert == 'A':
            mean = 20 if soil == 'clay' else 22 if soil == 'sandy' else 24
        elif fert == 'B':
            mean = 25 if soil == 'clay' else 23 if soil == 'sandy' else 27
        else:  # C
            mean = 22 if soil == 'clay' else 24 if soil == 'sandy' else 26
        growth = np.random.normal(mean, 2, 10)
        for g in growth:
            data.append({'growth': g, 'fertilizer': fert, 'soil': soil})
df2 = pd.DataFrame(data)
print("双因素数据预览:")
print(df2.head())
print(f"\n数据形状: {df2.shape}")

执行双因素ANOVA

# 执行双因素ANOVA（考虑交互作用）
model2 = ols('growth ~ C(fertilizer) + C(soil) + C(fertilizer):C(soil)', data=df2).fit()
anova_results2 = anova_lm(model2)
print("双因素方差分析结果:")
print(anova_results2)
# 如果交互作用不显著，可以使用加性模型（无交互）
model2_additive = ols('growth ~ C(fertilizer) + C(soil)', data=df2).fit()
anova_results2_additive = anova_lm(model2_additive)
print("\n加性模型（无交互）方差分析结果:")
print(anova_results2_additive)

重复测量ANOVA（Repeated Measures ANOVA）

# 创建重复测量数据
np.random.seed(42)
# 10个被试，在3个时间点的测量
subjects = 10
times = 3
data_rm = []
for subject in range(subjects):
    # 每个被试的基线水平不同
    baseline = np.random.normal(50, 10)
    for time in range(times):
        # 随时间变化
        measurement = baseline + time * 5 + np.random.normal(0, 5)
        data_rm.append({
            'subject': f'Subject_{subject+1}',
            'time': f'Time_{time+1}',
            'measurement': measurement
        })
df_rm = pd.DataFrame(data_rm)
# 使用混合模型进行重复测量ANOVA
from statsmodels.stats.anova import AnovaRM
# 执行重复测量ANOVA
rm_anova = AnovaRM(df_rm, 'measurement', 'subject', within=['time'])
rm_results = rm_anova.fit()
print("重复测量方差分析结果:")
print(rm_results)

方差分析前提假设检验

# 1. 正态性检验（Shapiro-Wilk test）
print("正态性检验 (Shapiro-Wilk):")
for fert in ['A', 'B', 'C']:
    stat, p_value = stats.shapiro(df[df['fertilizer'] == fert]['growth'])
    print(f"{fert}: 统计量={stat:.4f}, p值={p_value:.4f}")
# 2. 方差齐性检验（Levene's test）
from scipy import stats as scipy_stats
groups_data = [df[df['fertilizer'] == g]['growth'].values for g in ['A', 'B', 'C']]
levene_stat, levene_p = scipy_stats.levene(*groups_data)
print(f"\n方差齐性检验 (Levene):")
print(f"统计量={levene_stat:.4f}, p值={levene_p:.4f}")
# 3. 残差分析
model = ols('growth ~ C(fertilizer)', data=df).fit()
residuals = model.resid
fitted = model.fittedvalues
# Q-Q图检验正态性
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Q-Q图
stats.probplot(residuals, dist="norm", plot=axes[0])
axes[0].set_title('Q-Q图 - 残差正态性检验')
# 残差vs拟合值图
axes[1].scatter(fitted, residuals, alpha=0.7)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_xlabel('拟合值')
axes[1].set_ylabel('残差')
axes[1].set_title('残差 vs 拟合值')
plt.tight_layout()
plt.show()
# 4. 对非正态数据进行变换
from scipy import stats as scipy_stats
# Box-Cox变换
df['growth_boxcox'], lambda_val = scipy_stats.boxcox(df['growth'])
print(f"\nBox-Cox变换的lambda值: {lambda_val:.4f}")
# 变换后的正态性检验
for fert in ['A', 'B', 'C']:
    stat, p_value = scipy_stats.shapiro(df[df['fertilizer'] == fert]['growth_boxcox'])
    print(f"变换后 {fert}: 统计量={stat:.4f}, p值={p_value:.4f}")

非参数替代方法（Kruskal-Wallis检验）

from scipy.stats import kruskal
# 当方差分析假设不满足时的非参数替代
h_stat, p_value = kruskal(*groups_data)
print("Kruskal-Wallis非参数检验结果:")
print(f"H统计量: {h_stat:.4f}")
print(f"p值: {p_value:.4f}")
# 事后比较（Dunn检验）
from scikit_posthocs import posthoc_dunn
# 准备数据
df_dunn = df.copy()
# 执行Dunn检验
dunn_results = posthoc_dunn(df_dunn, val_col='growth', group_col='fertilizer', p_adjust='bonferroni')
print("\nDunn事后检验结果:")
print(dunn_results)

可视化结果

# 创建综合可视化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
# 1. 箱线图
sns.boxplot(x='fertilizer', y='growth', data=df, ax=axes[0, 0])
axes[0, 0].set_title('各肥料组生长情况分布')
# 2. 小提琴图
sns.violinplot(x='fertilizer', y='growth', data=df, ax=axes[0, 1])
axes[0, 1].set_title('分布密度图')
# 3. 点图+误差线
sns.pointplot(x='fertilizer', y='growth', data=df, 
              ci=95, capsize=0.2, ax=axes[0, 2])
axes[0, 2].set_title('均值与95%置信区间')
# 4. 双因素交互作用图
sns.pointplot(x='fertilizer', y='growth', hue='soil', data=df2, 
              ci=95, capsize=0.2, ax=axes[1, 0])
axes[1, 0].set_title('双因素交互作用图')
# 5. 频率分布直方图
for i, fert in enumerate(['A', 'B', 'C']):
    data_filtered = df[df['fertilizer'] == fert]['growth']
    axes[1, 1].hist(data_filtered, alpha=0.5, label=f'肥料{fert}', bins=10)
axes[1, 1].set_title('分布直方图')
axes[1, 1].legend()
# 6. 热力图显示组间差异
means_matrix = df2.pivot_table(values='growth', index='fertilizer', columns='soil', aggfunc='mean')
sns.heatmap(means_matrix, annot=True, fmt='.1f', cmap='YlOrRd', ax=axes[1, 2])
axes[1, 2].set_title('各组均值热力图')
plt.tight_layout()
plt.show()

完整ANOVA分析流程示例

def complete_anova_analysis(data, dependent_var, group_var, alpha=0.05):
    """
    完整的方差分析函数
    Parameters:
    -----------
    data : DataFrame
    dependent_var : str, 因变量
    group_var : str, 分组变量
    alpha : float, 显著性水平
    """
    results = {}
    # 1. 描述性统计
    results['descriptive'] = data.groupby(group_var)[dependent_var].describe()
    # 2. 执行ANOVA
    formula = f'{dependent_var} ~ C({group_var})'
    model = ols(formula, data=data).fit()
    anova_table = anova_lm(model)
    results['anova'] = anova_table
    # 3. 假设检验
    # 正态性
    shapiro_results = {}
    for group in data[group_var].unique():
        stat, p = scipy_stats.shapiro(data[data[group_var] == group][dependent_var])
        shapiro_results[group] = {'statistic': stat, 'p_value': p}
    results['normality'] = shapiro_results
    # 方差齐性
    groups_list = [data[data[group_var] == g][dependent_var].values 
                   for g in data[group_var].unique()]
    levene_stat, levene_p = scipy_stats.levene(*groups_list)
    results['homogeneity'] = {'statistic': levene_stat, 'p_value': levene_p}
    # 4. 如果需要，执行事后检验
    if anova_table['Pr(>F)'][0] < alpha:
        tukey = pairwise_tukeyhsd(data[dependent_var], data[group_var], alpha=alpha)
        results['posthoc'] = tukey
    return results
# 使用示例
results = complete_anova_analysis(df, 'growth', 'fertilizer')
print("描述性统计:")
print(results['descriptive'])
print("\n方差分析结果:")
print(results['anova'])
print("\n正态性检验:")
for group, stats in results['normality'].items():
    print(f"{group}: p值={stats['p_value']:.4f}")
print(f"\n方差齐性检验: p值={results['homogeneity']['p_value']:.4f}")