Python方差分析

昨天介绍了方差分析，以及方差分析之前的正态性检验和方差齐性检验，然后使用Excel和R语言两种工具分别来进行了方差的分析操作，具体的可以点击《方差分析 in R语言 and Excel》跳转，今天就使用Python来实现一下。

数据还是昨天那些数据，复制链接到浏览器就可以直接下载了↓

https://linss.oss-cn-beijing.aliyuncs.com/Excel/ANOVA_data.xlsx

数据的业务逻辑大概是，有四个省份，然后进行了一次培训。培训前后不同人员的销售金额的变化，需要看每个省份培训后是否带来销售额的提升，再决定后续是否加大培训，以及全国其他省份都进行培训的推广。

更多详细的可以去看上一篇。然后这一篇使用Python还是演示一下完整的过程，包括正态性检验，方差齐性检验和最后的方差分析。

【正态性检验】

首先加载相关的包，然后导入数据↓

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
from scipy.stats import norm
data = pd.read_excel("数据分布.xlsx",sheet_name="正态")

然后把数据绘制一个直方图和正态拟合曲线，直观上看一下整体数据的分布情况↓

# 绘制直方图
plt.hist(data['销量'], bins=30, density=True, alpha=0.5, color='#0494c4', edgecolor='black', label='直方图')
# 计算正态分布的参数
mu, sigma = norm.fit(data['销量'])
# 绘制正态拟合曲线
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, sigma)
plt.plot(x, p, 'k', linewidth=2, label='正态分布拟合曲线')


# 添加图例和标签
plt.legend()
plt.title('正态分布直方图 and 拟合曲线')
plt.xlabel('销量')
plt.ylabel('概率分布')
plt.show()

从图形上看整体分布还是挺正态的，但是不能相信眼睛，还是需要通过统计检验才行。这里还是和昨天R语言一样，使用两种检验方法。

【Shapiro-Wilk检验】

这是检验数据正态性的一种非常流行的方法，特别适用于样本量较小（通常n < 50）的情况。如果p值小于显著性水平（通常是0.05），则拒绝正态性的假设。

from scipy.stats import shapiro
import numpy as np
# 进行Shapiro-Wilk检验
statistic, p_value = shapiro(data['销量'])
print("Shapiro-Wilk检验统计量：", statistic)
print("p值：", p_value)


alpha = 0.05
if p_value > alpha:
    print("p值大于显著性水平，接受原假设，数据可能服从正态分布。")
else:
    print("p值小于显著性水平，拒绝原假设，数据不服从正态分布。")

上面是检验的代码，主要是scipy里面的stats库的使用，结果如下，可以看到虽然我们眼睛看上去数据是正态的，但结果却不服从正太分布↓

• Shapiro-Wilk检验统计量：0.9919911026954651
• p值：0.0002547771146055311
• p值小于显著性水平，拒绝原假设，数据不服从正态分布。

【Kolmogorov-Smirnov检验】

这个检验适用于样本量较大的情况，通过比较样本分布与正态分布之间的差异。同样，如果p值小于显著性水平，则拒绝数据来自正态分布的假设。

from scipy.stats import kstest
import numpy as np


# 计算样本数据的均值和标准差
sample_mean = np.mean(data['销量'])
sample_std = np.std(data['销量'])


# 进行Kolmogorov-Smirnov检验
statistic, p_value = kstest(data['销量'], 'norm', args=(sample_mean, sample_std))


print("Kolmogorov-Smirnov检验统计量：", statistic)
print("p值：", p_value)
alpha = 0.05
if p_value > alpha:
    print("p值大于显著性水平，接受原假设，数据可能服从正态分布。")
else:
    print("p值小于显著性水平，拒绝原假设，数据不服从正态分布。")

• Kolmogorov-Smirnov检验统计量：0.04446265579894415
• p值：0.08207083590829756
• p值大于显著性水平，接受原假设，数据可能服从正态分布。

从两种结果也可以看出，第一个Shapiro-Wilk检验是小样本，我们这里有800个样本，所以结果是第二个对大样本检验的看KS检验通过了。

不过这只是个参考，我们还是单独对每个省每次的数据进行检验，只要这部分通过了检验就可以了↓

#按省份和次数进行正态性检验
for pro in ['四川省','河南省','山东省','江西省']:
    for times in ['第1次','第2次']:
        pro_df = data[(data['省份'] == pro) & (data['次数'] == times)]
        sample_mean = np.mean(pro_df['销量'])
        sample_std = np.std(pro_df['销量'])
        statistic, p_value = shapiro(pro_df['销量'])
        statistic, p_value1 = kstest(pro_df['销量'], 'norm', args=(sample_mean, sample_std))
        print(f"Shapiro-Wilk检验统计量，{pro}_{times}_p值：", p_value)
        print(f"Kolmogorov-Smirnov检验，{pro}_{times}_p值：", p_value1)
    print()

这里写两个循环就行了，然后分别使用两种检验方法进行检验，最终的输出结果如下↓

结果是很完美的，全部通过检验，于是可以进行下一步的方差齐性检验了。

【方差齐性检验】

方差齐性（Homoscedasticity）指的是在不同的群组、条件或时间点上，变量的方差保持一致的性质。这一假设在许多统计分析方法中都非常重要，特别是在进行线性回归分析、方差分析（ANOVA）等时，方差齐性是关键的假设之一。

还是昨天的两种检验方法，我们还是先对整体做一个打样。

【Levene检验】

Levene检验是一种非参数的方法，它通过对各组数据的方差进行比较来检验方差齐性。该方法对数据的分布假设较为宽松，因此在数据不满足正态分布假设时也能有效地进行方差齐性检验。

# 进行Levene检验
statistic, p_value = levene(df_time1['销量'], df_time2['销量'])


print("Levene检验统计量：", statistic)
print("p值：", p_value)


alpha = 0.05
if p_value > alpha:
    print("p值大于显著性水平，接受原假设，数据具有方差齐性。")
else:
    print("p值小于显著性水平，拒绝原假设，数据不具有方差齐性。")

• Levene检验统计量：0.0326868006037593
• p值：0.8565744615389395
• p值大于显著性水平，接受原假设，数据具有方差齐性。

从结果可以看出，整体数据两次的方差齐性检验是很完美的通过了。

【Bartlett检验】

Bartlett检验是一种基于正态分布假设的方法，它要求各组数据满足正态分布。与Levene检验相比，Bartlett检验对数据的正态性要求更为严格。

from scipy.stats import bartlett
# 进行Bartlett检验
statistic, p_value = bartlett(df_time1['销量'], df_time2['销量'])


print("Bartlett检验统计量：", statistic)
print("p值：", p_value)


alpha = 0.05
if p_value > alpha:
    print("p值大于显著性水平，接受原假设，数据具有方差齐性。")
else:
    print("p值小于显著性水平，拒绝原假设，数据不具有方差齐性。")

• Bartlett检验统计量：0.638368357166069
• p值：0.4243022569404873
• p值大于显著性水平，接受原假设，数据具有方差齐性。

也是很完美的通过了检测，但是我们更关心的是每个省份是否两次的方差一致，所以还是需要分别对各省进行检验↓

for pro in ['四川省','河南省','山东省','江西省']:
    pro_df = data[(data['省份'] == pro)]
    df_pro_time1 = pro_df[pro_df['次数']=="第1次"]
    df_pro_time2 = pro_df[pro_df['次数']=="第2次"]


    statistic, p_value = levene(df_pro_time1['销量'], df_pro_time2['销量'])
    statistic, p_value1 = bartlett(df_pro_time1['销量'], df_pro_time2['销量'])
    print(f"Levene检验统计量，{pro}_p值：", p_value)
    print(f"Bartlett检验量，{pro}_p值：", p_value1)
    print()

结果还是很完美，两种检验都通过了方差齐性检验。最后我们终于可以进行方差分析了。

【方差分析】

首先导入数据↓

import numpy as np
import pandas as pd
from scipy.stats import f_oneway
import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
import os
os.chdir(r"G:\Code_s\PowerBI")
data = pd.read_excel("数据分布.xlsx",sheet_name="正态")

然后看一数据的分布情况，画一个箱线图直观看看变化↓

plt.figure(figsize=(10, 6))
sns.boxplot(data=data, x='省份', y='销量', hue='次数')
plt.title('四个省份两个阶段的数据分布',size=26)
plt.xlabel('省份', size = 16)
plt.ylabel('销量', size = 16)
plt.xticks(fontsize=14)
plt.show()

结果如上，基本上有个大概判断了，四川两次变化明显，其他省份变化并不明显。接下来通过方差分析，来验证我们的假想。首先看一下整体的两次变化↓

df_time1 = data[data['次数']=="第1次"]
df_time2 = data[data['次数']=="第2次"]
# 进行方差分析
f_statistic, p_value = f_oneway(df_time1['销量'], df_time2['销量'])
print("F统计量：", f_statistic)
print("p 值：", p_value)


alpha = 0.05
if p_value < alpha:
    print("拒绝原假设，不同组的均值存在显著差异")
else:
    print("接受原假设，不同组的均值没有显著差异")

• F统计量：2.663355928582429
• p 值：0.10307840303171252
• 接受原假设，不同组的均值没有显著差异

整体结果表示两次培训前后没有变化，但是我们更关心每个省份的情况，于是我们需要写一个循环，来看看各省的情况↓

for pro in ['四川省','河南省','山东省','江西省']:
    pro_df = data[(data['省份'] == pro)]
    df_pro_time1 = pro_df[pro_df['次数']=="第1次"]
    df_pro_time2 = pro_df[pro_df['次数']=="第2次"]


    f_statistic, p_value = f_oneway(df_pro_time1['销量'], df_pro_time2['销量'])
    print(f"{pro}_p值：", p_value)
    alpha = 0.05
    if p_value < alpha:
        print("拒绝原假设，不同组的均值存在显著差异")
    else:
        print("接受原假设，不同组的均值没有显著差异")
    print()

结果和我们预期一样，只有四川存在显著差异，其他几个省份培训前后都没什么变化。

链接是我使用PowerBI整合的历史文章，按类型分类，可以根据需求查询：Microsoft Power BI↓

https://app.powerbi.com/view?r=eyJrIjoiNjI2NWQ3NjktYjU0ZC00ZWZhLTgzMDgtMGI4ZTk1ZDlkODM3IiwidCI6IjI3NDQ3MWQ0LTM4ZDQtNDVlZS1hMmJkLWU1NTVhOTBkYzM4NiJ9

End