python统计分析——单因素方差分析

参考资料：用python动手学统计学

方差分析：analysis of variance，缩写为ANOVA

1、背景知识

1.1 要使用方差分析，数据的总体必须服从正态分布，而且各个水平内部的方差必须相等。

1.2 反复检验导致显著性结果更易出现的问题叫作多重假设检验问题：设显著性水平为0.05，则出现第一类错误的概率为5%。如果连续进行2次检验，则每次检验的显著性水平都为0.05，在这种情况下，出现第一类错误的概率就是1-0.95*0.95=0.0975，约为10%，超过了5%。检验次数越多，越容易拒绝零假设，也就越容易出现第一类错误。而方差分析可以有效的解决这一问题，通过一次检验完成显著性判断。

1.3 方差分析的F值=效应的方差/误差的方差。如果F值大，就认为效应比误差的影响大。当总体服从同方差正态分布，F值的样本分布就叫作F分布。通过F分布的累计分布函数计算p值，当p值小于0.05时，就拒绝零假设。

1.4 如下面的小提琴图所示：小提琴之间的高度差，即效应的大小，叫作组间差异。各个小提琴间的高度，即误差的大小，叫作组内差异。

2、导入库

# 导入库
# 用于数据计算的库
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
# 用于绘图的库
from matplotlib import pyplot as plt
import seaborn as sns
sns.set
# 用于估计统计模型的库
import statsmodels.formula.api as smf
import statsmodels.api as sm

3、数据准备

# 数据准备
weather=['cloudy','cloudy','rainy','rainy','sunny','sunny']
beer=[6,8,2,4,10,12]
data=pd.DataFrame({'beer':beer,
                  'weather':weather})
print(data)

4、数据展示

由于样本量很小，所以我们绘制箱线图，而非小提琴图。

箱线图和小提琴图的绘制的参数设置参考：

python统计分析——箱线图（sns.boxplot）_python sns.boxplot-CSDN博客

python统计分析——箱线图（plt.boxplot）_python boxplot-CSDN博客

python统计分析——箱线图（df.boxplot）_df.boxplot()-CSDN博客

python统计分析——小提琴图（sns.violinplot）-CSDN博客

python统计分析——小提琴图（plt.violinplot）_plt violin-CSDN博客

sns.boxplot(x='weather',y='beer',data=data)

5、查看均值

print(data.groupby('weather').mean())

6、常规步骤进行方差分析

excel操作步骤可参考：excel统计分析——单因素方差分析_处理间平方和简化公式-CSDN博客

6.1 计算组间平方和和组内平方和

# 天气的影响
effect=[7,7,3,3,11,11]
# 组间平方和
mu_effect=np.mean(effect)
squares_model=np.sum((effect-mu_effect)**2)
print(squares_model)

# 从原始数据中减掉效应就是误差
resid=data.beer-effect
# 组内平方和
squares_resid=np.sum(resid**2)
print(squares_resid)

6.2 计算组间方差和组内方差

# 组间差异的自由度
df_model=2
# 组内差异的自由度
df_resid=6-1-2
# 组间方差
variance_model=squares_model/df_model
# 组内方差
variance_resid=squares_resid/df_resid
print('组间方差：',variance_model)
print('组内方差：',variance_resid)

6.3 计算F值和p值

f_ratio=variance_model/variance_resid
print('F值：',f_ratio)
# 用F分布的累计分布函数计算p值
p_value=1-stats.f.cdf(x=f_ratio,dfn=df_model,dfd=df_resid)
print('p值：',p_value)

p值小于0.05，认为weather对beer的影响显著。注：原则上我们应该使用容量更大的样本。

7、解释变量为分类变量的正态线性模型

本例的线性模型如下：

$beer\sim N(\beta_0+\beta_1\times rainy+\beta_2\times sunny, \sigma^2)$

变量“rainy”在天气为rainy时为1，在其余天气时为0；变量sunny同理。参数β1代表rainy的影响程度，参数β2代表sunny的影响程度。当然还有cloudy天气，当rainy和sunny均为0时，只剩下β0，它代表cloudy的影响程度。

8、使用statsmodels进行方差分析

8.1 方差分析表

# 无论解释变量是连续变量还是分类变量，都要使用smf.ols定义
anova_model=smf.ols('beer~weather',data=data).fit()
# 模型拟合完成后，就可以通过sm.stats.anova_lm函数进行方差分析
print(sm.stats.anova_lm(anova_model,typ=2))

8.2 模型系数

# 输出模型的系数
anova_model.params

beer=7.0-4.0*rainy+4.0*sunny，解释为：当weather为cloudy时，beer=7.0-4.0*0+4.0*0=7.0；当weather为rainy时，beer=7.0-4.0*1+4.0*0=3.0；当weather为sunny时，beer=7.0-4*0+4*1=11。

# 输出拟合值
print('拟合值：')
print(anova_model.fittedvalues)
# 输出残差
print("残差：")
print(anova_model.resid)

8.3 回归模型的方差分析

正态线性模型中广泛应用了方差分析。当解释变量为连续变量时，方差分析依然有效。此时组间称为回归，组内称为残差。

# 准备连续型解释变量
data=pd.DataFrame({
    'beer':np.array([45.3, 59.3, 40.4, 38. , 37. , 40.9, 60.2, 63.3, 51.1, 44.9, 47. ,
                     53.2, 43.5, 53.2, 37.4, 59.9, 41.5, 75.1, 55.6, 57.2, 46.5, 35.8,
                     51.9, 38.2, 66. , 55.3, 55.3, 43.3, 70.5, 38.8]),
    'temp':np.array([20.5, 25. , 10. , 26.9, 15.8,  4.2, 13.5, 26. , 23.3,  8.5, 26.2,
                     19.1, 24.3, 23.3,  8.4, 23.5, 13.9, 35.5, 27.2, 20.5, 10.2, 20.5,
                     21.6,  7.9, 42.2, 23.9, 36.9,  8.9, 36.4,  6.4])
})
# 估计模型
lm_model=smf.ols('beer~temp',data=data).fit()
# 打印方差分析表
print(sm.stats.anova_lm(lm_model,typ=2))