Python统计实战：一题搞定双因子方差分析（交互效应分析）

为了解决特定问题而进行的学习是提高效率的最佳途径。这种方法能够使我们专注于最相关的知识和技能，从而更快地掌握解决问题所需的能力。

（以下练习题来源于《统计学—基于Python》。联系获取完整数据和Python源代码文件。）

练习题

城市道路交通管理部门为研究不同路段和不同时段对行车时间的影响，派一名交通警察分别在3个路段的高峰期与非高峰期亲自驾车进行实验，通过实验共获得30个行车时间（单位：分钟）数据。

（1）分析路段、时段以及路段和时段的交互作用对行车时间的影响（α=0.05）。

（2）计算偏效应量并进行分析。

计算结果和分析

（1）设路段（因子A）对行车时间的效应分别为α1（路段1）、α2（路段2）、α3（路段3）；时段（因子B）对行车时间的效应分别为β1（高峰期）、β2（非高峰期）。

检验路段效应的假设为：H0 = α1 = α2 = α3 = 0（路段对行车时间的影响不显著）；H1：α1，α2，α3至少有一个不等于0（路段对行车时间的影响显著）

检验时段效应的假设为： H0 = β1 = β2 = 0（路段对行车时间的影响不显著）；H1：β1，β2至少有一个不等于0（路段对行车时间的影响显著）

交互效应方差分析结果如下图所示。可以发现，检验时段和路段的P值均小于0.05，表示两个因子对行车时间的影响均显著；而检验交互效应的P值大于0.05，表示交互效应对产量的影响不显著。

# 主效应方差分析
# 将表格的短格式数据转为长格式
import pandas as pd
exercise9_5 = pd.read_csv('exercise9_5.csv')
# exercise9_4.head()

exercise9_5_1 = pd.melt(exercise9_5, id_vars = ['时段'], value_vars = ['路段1', '路段2', '路段3'],
                       var_name = '路段', value_name = '行车时间')

# 输出主效应方差分析表
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

model2i = ols(formula = '行车时间 ~ 路段 + 时段 + 路段:时段', data = exercise9_5_1).fit()  
anova_lm(model2i, typ = 1) # 输出交互效应方差分析表

（2）主效应量与偏效应量如计算结果如下图所示。交互效应的偏效应量为0.000244，表示在排除时段和路段的影响后，交互效应解释了行车时间误差的0.0244%。如此小的偏效应量也佐证了（1）题检验中交互效应不显著的结论。

# 计算交互效应方差分析的效应量
import numpy as np
model2i = ols(formula = '行车时间 ~ 路段 + 时段 + 路段:时段', data = exercise9_5_1).fit()  
df_res = anova_lm(model2i, typ = 1) 
df_res['eta_sq'] = df_res['sum_sq'] / df_res['sum_sq'].sum()  # 计算主效应量
df_res['eta_sq_part'] = df_res['sum_sq'] / (df_res['sum_sq'] + df_res.loc['Residual', 'sum_sq']) # 计算偏效应量
df_res.loc['Residual', 'eta_sq_part'] = np.nan
df_res[['eta_sq', 'eta_sq_part']][:3]  # 输出前3行效应量

绘制主效应和交互效应图

下图中的箱线图反映了两个因子的主效应，可观察路段和时段的不同处理的行车实践差异。折线图反映了两个因子的交互效应。图中的各条折现基本上是平行的，表示两个因子间无明显的交互效应。

# 绘制主效应和交互效应图
from statsmodels.graphics.api import interaction_plot
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Songti SC'] # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False     # 正常显示负号
df = exercise9_5_1

y_col = '行车时间'
x1_col = '时段'
x2_col = '路段'

plt.subplots(2, 2, figsize = (8, 7))
ax1 = plt.subplot(221)
interaction_plot(x = df[x2_col], trace = df[x1_col], response = df[y_col],
                xlabel = y_col, ylabel = x2_col, ax = ax1)
ax1.set_title(f'{y_col} ~ {x2_col}|{x1_col}')

ax2 = plt.subplot(222)
df.boxplot(column = y_col, by = x1_col, grid = False, ax = ax2)
ax2.set_title(f'{y_col} ~ {x1_col}|{x1_col}')

ax3 = plt.subplot(223)
df.boxplot(column = y_col, by = x2_col, grid = False, ax = ax3)
ax3.set_title(f'{y_col} ~ {x2_col}|{x2_col}')

ax4 = plt.subplot(224)
interaction_plot(x = df[x1_col], trace = df[x2_col], response = df[y_col],
                xlabel = y_col, ylabel = x1_col, ax = ax4)
ax4.set_title(f'{y_col} ~ {x1_col}|{x2_col}')

plt.tight_layout()

都读到这里了，不妨关注、点赞一下吧！