Python statsmodels 0.14 双因素方差分析实战：从数据导入到结果解读

在数据分析的实际工作中，我们常常需要同时考察两个因素对结果变量的影响。比如在农业试验中，研究不同品种和不同施肥量对作物产量的影响；在市场营销中，分析不同广告渠道和促销力度对销售额的作用。这时候，双因素方差分析（Two-Way ANOVA）就成为了一个强有力的工具。

本文将带你使用Python的statsmodels 0.14库，通过一个完整的案例，一步步实现双因素方差分析的全流程。不同于教科书上的理论推导，我们聚焦于实际应用场景中的代码实现和结果解读，让你能够快速将这一统计方法应用到自己的项目中。

1. 环境准备与数据加载

在开始分析之前，我们需要确保环境配置正确，并准备好待分析的数据集。statsmodels是Python中一个强大的统计建模库，特别适合进行各种方差分析。

首先安装必要的库（如果尚未安装）：

pip install statsmodels pandas numpy

假设我们有一个关于植物生长的实验数据集，研究两种因素对植物高度的影响：

• 光照条件（低、中、高）
• 肥料类型（A、B）

import pandas as pd import statsmodels.api as sm from statsmodels.formula.api import ols # 创建示例数据集 data = { 'height': [12, 14, 16, 13, 15, 18, 14, 17, 20, 15, 18, 21, 16, 19, 23, 17, 20, 24], 'light': ['low', 'low', 'low', 'medium', 'medium', 'medium', 'high', 'high', 'high'] * 2, 'fertilizer': ['A']*9 + ['B']*9 } df = pd.DataFrame(data)

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

当两个因素之间没有交互作用，或者我们只关心主效应时，可以使用无重复双因素方差分析。这种情况下，每个因素组合只需要一个观测值。

2.1 模型构建与拟合

使用statsmodels的OLS（普通最小二乘法）接口来构建方差分析模型：

# 无重复双因素方差分析模型 model = ols('height ~ C(light) + C(fertilizer)', data=df).fit() anova_table = sm.stats.anova_lm(model, typ=2) print(anova_table)

2.2 结果解读

运行上述代码后，你将得到一个标准的方差分析表：

解读要点：

• P值：若小于显著性水平（通常0.05），则认为该因素对结果有显著影响
• F值：反映因素效应与随机误差的比值，越大表示影响越显著
• 自由度：light有3个水平，自由度为2；fertilizer有2个水平，自由度为1

2.3 可视化分析结果

为了更直观地理解因素效应，我们可以绘制均值图：

import seaborn as sns import matplotlib.pyplot as plt plt.figure(figsize=(10, 5)) sns.pointplot(x='light', y='height', hue='fertilizer', data=df, ci=95, dodge=True) plt.title('植物高度在不同光照和肥料条件下的均值') plt.show()

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

当我们需要考虑两个因素之间的交互作用时，就要使用可重复双因素方差分析。这意味着每个因素组合需要有多个观测值（重复实验）。

3.1 数据准备

扩展我们的数据集，加入交互作用：

# 扩展数据集，包含重复测量 data_interaction = { 'height': [12, 14, 11, 16, 15, 17, 13, 15, 14, 18, 17, 19, 14, 16, 15, 20, 19, 21, 15, 17, 16, 21, 20, 22, 16, 18, 17, 23, 22, 24, 17, 19, 18, 24, 23, 25], 'light': ['low']*6 + ['medium']*6 + ['high']*6 + ['low']*6 + ['medium']*6 + ['high']*6, 'fertilizer': ['A']*18 + ['B']*18, 'replicate': [1,2,3]*12 } df_interaction = pd.DataFrame(data_interaction)

3.2 包含交互项的模型

在模型公式中加入交互项（使用*符号）：

# 可重复双因素方差分析模型（含交互项） model_interaction = ols('height ~ C(light) * C(fertilizer)', data=df_interaction).fit() anova_table_interaction = sm.stats.anova_lm(model_interaction, typ=2) print(anova_table_interaction)

3.3 交互作用分析

得到的方差分析表将包含交互作用项：

交互作用显著（P=0.0012）意味着光照和肥料对植物高度的影响不是独立的。我们需要进一步分析这种交互模式。

3.4 交互作用可视化

使用交互作用图来直观展示：

plt.figure(figsize=(10, 6)) sns.pointplot(x='light', y='height', hue='fertilizer', data=df_interaction, ci=95, dodge=True) plt.title('光照和肥料对植物高度的交互作用') plt.show()

4. 模型诊断与验证

在进行方差分析后，我们需要验证模型假设是否满足，确保分析结果的可靠性。

4.1 正态性检验

检查残差是否服从正态分布：

from scipy import stats import numpy as np residuals = model_interaction.resid _, p_value = stats.normaltest(residuals) print(f'正态性检验P值: {p_value:.4f}') # Q-Q图 sm.qqplot(residuals, line='s') plt.title('残差Q-Q图') plt.show()

4.2 方差齐性检验

检查各组方差是否相等（Levene检验）：

from statsmodels.stats.diagnostic import het_breuschpagan _, p_value, _, _ = het_breuschpagan(residuals, model_interaction.model.exog) print(f'方差齐性检验P值: {p_value:.4f}')

4.3 离群值检测

识别可能的异常观测值：

influence = model_interaction.get_influence() cooks_d = influence.cooks_distance[0] plt.stem(cooks_d) plt.title("Cook's距离检测离群值") plt.xlabel('观测序号') plt.ylabel("Cook's距离") plt.show()

5. 事后检验与多重比较

当主效应显著时，我们需要进一步了解哪些组别之间存在显著差异。

5.1 Tukey HSD检验

对所有组别进行两两比较：

from statsmodels.stats.multicomp import pairwise_tukeyhsd tukey = pairwise_tukeyhsd(endog=df_interaction['height'], groups=df_interaction['light'] + '_' + df_interaction['fertilizer'], alpha=0.05) print(tukey.summary())

5.2 简单效应分析

当交互作用显著时，可以分析一个因素在不同水平下的简单效应：

# 在不同光照水平下分析肥料效应 for light_level in ['low', 'medium', 'high']: subset = df_interaction[df_interaction['light'] == light_level] model = ols('height ~ C(fertilizer)', data=subset).fit() print(f'\n光照水平: {light_level}') print(sm.stats.anova_lm(model, typ=2))

6. 实际应用中的注意事项

在实际项目中应用双因素方差分析时，有几个关键点需要特别注意：

1. 实验设计阶段：

   • 确保样本量足够，特别是考虑交互作用时
   • 随机化实验顺序，避免混淆变量影响
   • 平衡设计（各组合样本量相同）能提高检验效能

2. 数据分析阶段：

   • 总是先检查交互作用，再解释主效应
   • 当交互作用显著时，主效应的解释需要谨慎
   • 考虑使用效应量（如η²）补充P值信息

3. 结果报告阶段：

   • 报告完整的方差分析表
   • 包括F值、自由度、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("效应量(η²):\n", eta_squared)

7. 扩展应用与进阶技巧

掌握了基本方法后，我们可以进一步探索一些高级应用场景：

7.1 处理非平衡设计

当各组样本量不等时，需要使用Type III平方和：

# 非平衡数据的方差分析 model_unbalanced = ols('height ~ C(light) * C(fertilizer)', data=df_unbalanced).fit() anova_unbalanced = sm.stats.anova_lm(model_unbalanced, typ=3) print(anova_unbalanced)

7.2 混合效应模型

当数据存在嵌套结构或重复测量时，可以考虑混合效应模型：

import statsmodels.formula.api as smf # 混合效应模型示例 mixed_model = smf.mixedlm('height ~ C(light)*C(fertilizer)', data=df, groups=df['block']) result = mixed_model.fit() print(result.summary())

7.3 协方差分析(ANCOVA)

当需要考虑连续型协变量时，可以使用协方差分析：

# 包含协变量的模型 ancova_model = ols('height ~ C(light) + C(fertilizer) + initial_height', data=df).fit() print(sm.stats.anova_lm(ancova_model, typ=2))