Python statsmodels 0.14 双因素方差分析实战：3步完成可重复/无重复实验检验

沃克森

266人浏览 · 2026-07-03 09:21:15

沃克森 · 2026-07-03 09:21:15 发布

Python statsmodels 0.14 双因素方差分析实战：3步完成可重复/无重复实验检验

在数据分析领域，双因素方差分析（Two-Way ANOVA）是一种强大的统计工具，用于同时研究两个分类变量对连续型因变量的影响。与单因素方差分析相比，双因素分析不仅能评估各因素的独立效应，还能检测因素间的交互作用。本文将使用Python的statsmodels 0.14库，通过实际代码演示如何快速实现这一分析。

对于数据科学家和分析师而言，掌握双因素方差分析的自动化实现至关重要。传统手动计算不仅耗时且容易出错，而现代统计软件如statsmodels则能高效完成从数据准备到结果解读的全流程。我们将重点区分可重复（有交互）和无重复（无交互）两种实验设计，并提供可直接运行的Jupyter Notebook代码示例。

1. 环境准备与数据模拟

在开始分析前，我们需要确保环境配置正确并生成合适的测试数据。statsmodels 0.14对Python 3.8+有最佳支持，建议使用虚拟环境管理依赖。

# 基础环境配置
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt
import seaborn as sns

print(f"statsmodels版本: {sm.__version__}")

1.1 模拟无重复实验数据

无重复双因素实验意味着每个因素组合只有一个观测值，无法评估交互效应。我们模拟一个农业实验场景，研究肥料类型（A/B/C）和灌溉方式（X/Y）对作物产量的影响：

np.random.seed(42)

# 定义因素水平
fertilizers = ['A', 'B', 'C']
irrigation = ['X', 'Y']

# 生成模拟数据
data = []
for fert in fertilizers:
    for irr in irrigation:
        # 基础产量 + 肥料效应 + 灌溉效应 + 随机误差
        yield_base = 50
        fert_effect = {'A': 0, 'B': 5, 'C': 3}[fert]
        irr_effect = {'X': 0, 'Y': 4}[irr]
        noise = np.random.normal(0, 2)
        
        data.append({
            'Fertilizer': fert,
            'Irrigation': irr,
            'Yield': yield_base + fert_effect + irr_effect + noise
        })

df_no_replication = pd.DataFrame(data)

1.2 模拟可重复实验数据

可重复实验每个因素组合有多个观测值，可评估交互作用。我们扩展上述场景，每个组合测量3次：

np.random.seed(42)

data = []
for fert in fertilizers:
    for irr in irrigation:
        for rep in range(3):  # 每个组合3次重复
            # 基础设置
            yield_base = 50
            fert_effect = {'A': 0, 'B': 5, 'C': 3}[fert]
            irr_effect = {'X': 0, 'Y': 4}[irr]
            
            # 交互效应：特定组合有额外影响
            interaction = 0
            if fert == 'B' and irr == 'Y':
                interaction = 3
            elif fert == 'C' and irr == 'X':
                interaction = -2
                
            noise = np.random.normal(0, 2)
            
            data.append({
                'Fertilizer': fert,
                'Irrigation': irr,
                'Yield': yield_base + fert_effect + irr_effect + interaction + noise
            })

df_with_replication = pd.DataFrame(data)

2. 无重复双因素方差分析实现

无重复实验的分析模型较为简单，不考虑交互项。statsmodels提供了简洁的公式接口实现。

2.1 模型构建与拟合

使用ols函数指定模型公式，其中 C() 表示分类变量：

# 无重复实验模型
model_no_interaction = ols('Yield ~ C(Fertilizer) + C(Irrigation)', 
                          data=df_no_replication).fit()

2.2 结果解读与可视化

拟合完成后，使用anova_table查看方差分析结果：

anova_table = sm.stats.anova_lm(model_no_interaction, typ=2)
print(anova_table)

典型输出结果示例：

来源	平方和(SS)	自由度(df)	均方(MS)	F值	P值
Fertilizer	86.72	2	43.36	12.45	0.0079
Irrigation	32.04	1	32.04	9.20	0.0286
Residual	13.94	4	3.48

关键解读点：

Fertilizer效应 ：P=0.0079<0.05，表明肥料类型对产量有显著影响
Irrigation效应 ：P=0.0286<0.05，灌溉方式也有显著影响
无交互项 ：因为实验设计无法评估交互作用

可视化效应大小：

plt.figure(figsize=(10, 4))
sns.barplot(x='Fertilizer', y='Yield', hue='Irrigation', 
            data=df_no_replication, ci=None)
plt.title("无重复实验各因素组合的平均产量")
plt.show()

3. 可重复双因素方差分析实现

可重复实验能评估交互作用，模型需包含交互项。这是更全面的分析方法。

3.1 包含交互项的模型构建

在公式中使用 * 运算符自动包含主效应和交互效应：

model_with_interaction = ols('Yield ~ C(Fertilizer) * C(Irrigation)', 
                            data=df_with_replication).fit()

3.2 交互效应分析与解读

获取完整的方差分析表：

anova_table_interaction = sm.stats.anova_lm(model_with_interaction, typ=2)
print(anova_table_interaction)

典型输出示例：

来源	平方和(SS)	自由度(df)	均方(MS)	F值	P值
Fertilizer	167.56	2	83.78	25.32	<0.0001
Irrigation	124.33	1	124.33	37.58	<0.0001
Fertilizer:Irrigation	58.17	2	29.08	8.79	0.0038
Residual	39.72	12	3.31

关键发现：

显著交互作用 ：交互项P=0.0038<0.01，表明肥料效果依赖于灌溉方式
主效应依然显著 ：但需谨慎解释，因存在交互作用

交互效应可视化：

plt.figure(figsize=(10, 5))
sns.pointplot(x='Fertilizer', y='Yield', hue='Irrigation',
              data=df_with_replication, ci=95, dodge=True)
plt.title("可重复实验中的交互效应图示")
plt.show()

4. 模型诊断与进阶技巧

确保方差分析结果可靠，需进行模型诊断。statsmodels提供多种诊断工具。

4.1 残差分析

残差应满足正态性和方差齐性假设：

# 残差正态性检验
from scipy.stats import shapiro
residuals = model_with_interaction.resid
_, p_value = shapiro(residuals)
print(f"Shapiro-Wilk正态性检验P值: {p_value:.4f}")

# 残差图
plt.figure(figsize=(12, 4))
plt.subplot(121)
sns.histplot(residuals, kde=True)
plt.title("残差分布")

plt.subplot(122)
sns.scatterplot(x=model_with_interaction.fittedvalues, y=residuals)
plt.axhline(y=0, color='r', linestyle='--')
plt.title("拟合值 vs 残差")
plt.show()

4.2 多重比较校正

当主效应显著时，可能需要进行事后检验：

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 对肥料类型进行多重比较
tukey = pairwise_tukeyhsd(endog=df_with_replication['Yield'],
                         groups=df_with_replication['Fertilizer'],
                         alpha=0.05)
print(tukey.summary())

4.3 效应量计算

除P值外，计算效应量（如η²）更全面评估因素重要性：

def calculate_eta_squared(anova_table):
    total_ss = anova_table['sum_sq'].sum()
    return anova_table['sum_sq'] / total_ss

eta_squared = calculate_eta_squared(anova_table_interaction)
print("各效应项的η²值:")
print(eta_squared)