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

Python statsmodels 0.14 双因素方差分析实战从数据导入到结果解读在数据分析的实际工作中我们常常需要同时考察两个因素对结果变量的影响。比如在农业试验中研究不同品种和不同施肥量对作物产量的影响在市场营销中分析不同广告渠道和促销力度对销售额的作用。这时候双因素方差分析Two-Way ANOVA就成为了一个强有力的工具。本文将带你使用Python的statsmodels 0.14库通过一个完整的案例一步步实现双因素方差分析的全流程。不同于教科书上的理论推导我们聚焦于实际应用场景中的代码实现和结果解读让你能够快速将这一统计方法应用到自己的项目中。1. 环境准备与数据加载在开始分析之前我们需要确保环境配置正确并准备好待分析的数据集。statsmodels是Python中一个强大的统计建模库特别适合进行各种方差分析。首先安装必要的库如果尚未安装pip install statsmodels pandas numpy假设我们有一个关于植物生长的实验数据集研究两种因素对植物高度的影响光照条件低、中、高肥料类型A、Bimport 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), datadf).fit() anova_table sm.stats.anova_lm(model, typ2) print(anova_table)2.2 结果解读运行上述代码后你将得到一个标准的方差分析表来源平方和(SS)自由度(df)均方(MS)F值P值C(light)56.11228.0628.060.000016C(fertilizer)72.25172.2572.250.000001残差8.00140.57解读要点P值若小于显著性水平通常0.05则认为该因素对结果有显著影响F值反映因素效应与随机误差的比值越大表示影响越显著自由度light有3个水平自由度为2fertilizer有2个水平自由度为12.3 可视化分析结果为了更直观地理解因素效应我们可以绘制均值图import seaborn as sns import matplotlib.pyplot as plt plt.figure(figsize(10, 5)) sns.pointplot(xlight, yheight, huefertilizer, datadf, ci95, dodgeTrue) 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), datadf_interaction).fit() anova_table_interaction sm.stats.anova_lm(model_interaction, typ2) print(anova_table_interaction)3.3 交互作用分析得到的方差分析表将包含交互作用项来源平方和(SS)自由度(df)均方(MS)F值P值C(light)112.22256.1156.110.0001C(fertilizer)144.501144.50144.500.0001C(light):C(fertilizer)16.0028.008.000.0012残差36.00301.20交互作用显著P0.0012意味着光照和肥料对植物高度的影响不是独立的。我们需要进一步分析这种交互模式。3.4 交互作用可视化使用交互作用图来直观展示plt.figure(figsize(10, 6)) sns.pointplot(xlight, yheight, huefertilizer, datadf_interaction, ci95, dodgeTrue) 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, lines) 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(Cooks距离检测离群值) plt.xlabel(观测序号) plt.ylabel(Cooks距离) plt.show()5. 事后检验与多重比较当主效应显著时我们需要进一步了解哪些组别之间存在显著差异。5.1 Tukey HSD检验对所有组别进行两两比较from statsmodels.stats.multicomp import pairwise_tukeyhsd tukey pairwise_tukeyhsd(endogdf_interaction[height], groupsdf_interaction[light] _ df_interaction[fertilizer], alpha0.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), datasubset).fit() print(f\n光照水平: {light_level}) print(sm.stats.anova_lm(model, typ2))6. 实际应用中的注意事项在实际项目中应用双因素方差分析时有几个关键点需要特别注意实验设计阶段确保样本量足够特别是考虑交互作用时随机化实验顺序避免混淆变量影响平衡设计各组合样本量相同能提高检验效能数据分析阶段总是先检查交互作用再解释主效应当交互作用显著时主效应的解释需要谨慎考虑使用效应量如η²补充P值信息结果报告阶段报告完整的方差分析表包括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), datadf_unbalanced).fit() anova_unbalanced sm.stats.anova_lm(model_unbalanced, typ3) print(anova_unbalanced)7.2 混合效应模型当数据存在嵌套结构或重复测量时可以考虑混合效应模型import statsmodels.formula.api as smf # 混合效应模型示例 mixed_model smf.mixedlm(height ~ C(light)*C(fertilizer), datadf, groupsdf[block]) result mixed_model.fit() print(result.summary())7.3 协方差分析(ANCOVA)当需要考虑连续型协变量时可以使用协方差分析# 包含协变量的模型 ancova_model ols(height ~ C(light) C(fertilizer) initial_height, datadf).fit() print(sm.stats.anova_lm(ancova_model, typ2))