【Python・统计学】（两）多因素方差分析（原理及代码）_python多因素方差分析

前言

自学笔记，分享给对统计学原理不太清楚但需要在论文中用到的小伙伴，欢迎大佬们补充或绕道。ps：本文不涉及公式讲解（文科生小白友好体质）～（部分定义等来源于知乎）

本文重点：多因素方差分析

【1.多因素方差分析简单原理】

【2.多因素方差分析步骤】

【3.多因素方差分析代码（配对/独立+事后检验）】

1.多因素方差分析简单原理

   多因素方差分析，用于研究一个因变量是否受到多个自变量（也称为因素）的影响，它检验多个因素取值水平的不同组合之间，因变量的均值之间是否存在显著的差异。多因素方差分析既可以分析单个因素的作用（主效应），也可以分析因素之间的交互作用（交互效应），还可以进行协方差分析，以及各个因素变量与协变量的交互作用。

*主效应(Main Effect):

• 主效应指的是每个因素单独对观测指标(如作物产量)的影响。
• 在这个例子中,化肥类型的主效应就是不同化肥类型之间作物产量的差异,而土壤类型的主效应就是不同土壤类型之间作物产量的差异。
• 主效应反映了因素的总体效应,不考虑其他因素的影响。

*交互效应(Interaction Effect):

• 交互效应指的是两个或多个因素共同作用对观测指标产生的影响,这种影响不能简单地用各因素的主效应相加来解释。
• 在这个例子中,如果某种化肥在特定的土壤类型上特别有效,而在其他土壤类型上效果一般,这就说明化肥类型和土壤类型存在交互效应。
• 交互效应反映了因素之间的相互影响和依赖关系。

也就是说，

如果只有主效应显著,说明各因素独立地影响观测指标。

如果交互效应显著,说明因素之间存在相互影响,需要进一步分析交互效应的性质和方向。

2.多因素方差分析步骤

之前学的时候写的笔记（虽然是日语）

①交互作用の検定
②交互作用OK: 単純主効果の検定→３水準以上（多重比較）
　　　　　交互作用No:　主効果の検定→３水準以上（多重比較）

①首先在做多因素方差分析之前要看数据是否满足前提条件

  检验方差是否满足方差齐性，以及数据的正态性（详情参考之前的t检验和单因素方差分析）

②检验交互效应（因素之间是否存在显著交互）

• 如果交互效应显著,说明因素之间存在相互作用,需要进一步分析单纯主效应。
• 如果交互效应不显著,说明因素之间不存在显著的相互作用,可以直接分析主效应。

③(如果交互效应显著)分析单纯主效应

• 在每个因素的每个水平下,检验另一个因素的效应是否显著。
• 如果单纯主效应显著,并且因素有三个或更多水平,进行多重比较以确定哪些水平之间存在显著差异。

④(如果交互效应不显著)分析主效应:

• 检验每个因素的主效应是否显著。
• 如果主效应显著,并且因素有三个或更多水平,进行多重比较以确定哪些水平之间存在显著差异。

3.多因素方差分析以及事后比较（独立样本代码）

①读取数据，进行方差齐性检验

使用的是一组教育学数据进行两因素方差分析

有实验组(Group1)和对照组(Group2)进行两种不同的方法教学，其中再把她们分成10人，20人和30人的班（Size），集中教学一段时间后进行测试（成绩为Score）

检验group（教学方法）和班级规模（size）的两因素，是否对最后的学生测试成绩有影响

import pandas as pd
from scipy.stats import levene
from statsmodels.stats.anova import AnovaRM
 
# 从 CSV 文件读取数据
df = pd.read_csv('/content/taionasi.csv')
balanced_df = pd.DataFrame(columns=df.columns)
 
# 按照 group 进行抽样，确保每个组的样本大小相同
for group in df['group'].unique():
    group_data = df[df['group'] == group]
 
    # 对每个组进行抽样，使得样本大小相同
    min_size = min(group_data['size'])
    sampled_data = group_data.groupby('size', group_keys=False).apply(lambda x: x.sample(min_size, random_state=42))
 
    # 将抽样后的数据添加到平衡的 DataFrame 中
    balanced_df = balanced_df.append(sampled_data, ignore_index=True)
 
# 检查平衡后的数据结构
print(balanced_df.head())
 
# 对平衡后的数据进行 Levene's Test
grouped_data_balanced = [balanced_df[balanced_df['group'] == group]['score'].values for group in balanced_df['group'].unique()]
statistic, p_value = levene(*grouped_data_balanced)
# 打印 Levene's Test 的结果
print(f"Levene's Test Statistic: {statistic}")
print(f"P-value: {p_value}")
 
# 根据 P-value 判断方差齐性
alpha = 0.05  # 有意水準
if p_value < alpha:
    print('帰無仮説を棄却：データのグループ間に等分散性がない可能性が高いです。')
else:
    print('帰無仮説は棄却されません：データのグループ間に等分散性があると考えられます。')

*注意：检验方差齐性的时候为了确保每组数据的个数一样，进行了数据抽样。如果不进行数据抽样虽然也不会报错但是数据结果截然相反。 

结果：抽样数据显示满足方差齐性，可以进行多因素方差分析 

  group size score
0     1   10    62
1     1   10    70
2     1   10    83
3     1   10    61
4     1   10    79
Levene's Test Statistic: 1.2774939867588024
P-value: 0.26301796308207803
帰無仮説は棄却されません：データのグループ間に等分散性があると考えられます。

②用图表检查是因素之间是否存在交互效应 

#交互作用プロット
from statsmodels.graphics.factorplots import interaction_plot
from matplotlib import pyplot as plt # import matplotlib.pyplot as plt グラフを描画するライブラリ
 
fig = interaction_plot(df['size'],  #x軸のデータ
                       df['group'],  #層別するデータ
                       df['score'],  #y軸のデータ
                       colors=['red', 'blue'],  #グラフの色を設定
                       ms=15)  #グラフの点の大きさの設定
plt.show()

 数据存在交点，证明存在交互效应（如果没有交互效应，两条线是平行的）

 ③进行多因素方差分析

• 导入所需的库:statsmodels
• 使用ols函数拟合线性模型,模型公式为'score ~ group + size + group:size',表示以score为因变量,group和size为自变量,group:size表示group和size的交互效应
• 使用anova_lm函数执行方差分析,type=2表示使用II型平方和
• 将结果表的列名替换为日语
• 打印方差分析结果表,并注明"交互作用ありの結果"

#二元配置分散分析（交互作用あり）
# 統計モデルを推定するライブラリ
import statsmodels.formula.api as smf
import statsmodels.api as sm
 
anova_mod = smf.ols('score ~ group + size + group:size', data=df).fit()
taionasi=sm.stats.anova_lm(anova_mod, type=2)
taionasi.columns = ["自由度","平方和","平均平方和","F値","p値"] #列名を日本語に差し替え
print("交互作用ありの結果："+"\n"+str(taionasi))

交互作用ありの結果：
              自由度           平方和        平均平方和         F値            p値
group         1.0     26.450000    26.450000   0.261825  6.095178e-01
size          2.0  11238.033333  5619.016667  55.621907  2.108493e-19
group:size    2.0   4057.500000  2028.750000  20.082329  1.421430e-08
Residual    174.0  17577.766667   101.021648        NaN           NaN

 ④事后检验

• 导入所需的库:statsmodels、itertools、pandas
• 将group和size列的数据类型转换为字符串
• 初始化一个空的数据框combined_results,用于存储多重比较结果
• 定义函数run_tukey_comparison,用于执行Tukey's HSD多重比较
  • 对每个group的不同size水平进行两两比较
  • 使用pairwise_tukeyhsd函数进行Tukey's HSD多重比较
  • 提取多重比较结果的重要信息(均值差、置信区间、是否拒绝原假设),并存储到comparison_results数据框中
  • 将comparison_results添加到combined_results数据框中
• 分别对group1和group2执行run_tukey_comparison函数,进行多重比较
• 调整combined_results数据框的列顺序
• 打印合并的多重比较结果

from statsmodels.stats.multicomp import pairwise_tukeyhsd
from itertools import combinations
import pandas as pd
from IPython.display import display, HTML
 
# df 是你的数据框
# 调整数据类型
df['group'] = df['group'].astype(str)
df['size'] = df['size'].astype(str)
 
# 初始化一个空的数据框，用于存储多重比较结果
combined_results = pd.DataFrame(columns=['Group', 'Comparison', 'Mean Difference', 'Lower CI', 'Upper CI', 'Reject Null'])
 
# 定义函数执行多重比较
def run_tukey_comparison(group_label, group_data, size_values):
    group_combinations = combinations(group_data, 2)
    for size1, size2 in group_combinations:
        m_comp = pairwise_tukeyhsd(endog=pd.concat([size1['score'], size2['score']]),
                                   groups=pd.concat([pd.Series([f'{group_label}_size_{size_values[size1.iloc[0]["size"]]}'] * len(size1)),
                                                    pd.Series([f'{group_label}_size_{size_values[size2.iloc[0]["size"]]}'] * len(size2))]),
                                   alpha=0.05)
        # 提取多重比较结果的重要信息
        comparison_results = pd.DataFrame({
            'Group': [f'{group_label}'],
            'Comparison': [f"{size1.iloc[0]['size']} vs {size2.iloc[0]['size']}"],
            'Mean Difference': [m_comp.meandiffs[0]],
            'Lower CI': [m_comp.confint[0, 0]],
            'Upper CI': [m_comp.confint[0, 1]],
            'Reject Null': [m_comp.reject[0]],
        })
        # 将结果添加到总的数据框中
        global combined_results
        combined_results = pd.concat([combined_results, comparison_results], ignore_index=True)
 
# 运行 group1 的多重比较
run_tukey_comparison('group1', group1_sizes, {'10': 10, '20': 20, '30': 30})
 
# 运行 group2 的多重比较
run_tukey_comparison('group2', group2_sizes, {'10': 10, '20': 20, '30': 30})
 
# 调整表格显示顺序
combined_results = combined_results[['Group', 'Comparison', 'Mean Difference', 'Lower CI', 'Upper CI', 'Reject Null']]
 
# 打印合并结果
print("Combined Results:")
display(HTML(combined_results.to_html(index=False)))

事后检验结果：

Group	Comparison	Mean Difference	Lower CI	Upper CI	Reject Null
group1	10 vs 20	-9.533333	-13.918702	-5.147965	True
group1	10 vs 30	-30.266667	-35.333353	-25.199980	True
group1	20 vs 30	-20.733333	-25.577774	-15.888892	True
group2	10 vs 20	-9.533333	-13.918702	-5.147965	True
group2	10 vs 30	-30.266667	-35.333353	-25.199980	True
group2	20 vs 30	-20.733333	-25.577774	-15.888892	True

总的来说,这些结果表明:

1. 在 Group1 和 Group2 中,size 水平为 20 的组的平均值显著高于 size 水平为 10 的组。
2. 在 Group1 和 Group2 中,size 水平为 30 的组的平均值显著高于 size 水平为 10 的组。
3. 在 Group1 和 Group2 中,size 水平为 30 的组的平均值显著高于 size 水平为 20 的组。

   这些发现表明,不同的 size 水平对因变量有显著影响,并且 size 水平越高,因变量的值越高。此外,Group1 和 Group2 表现出相同的模式,这表明 group 因素可能对 size 水平之间的差异没有显著影响。