python统计分析单因素方差分析
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 sm3、数据准备
# 数据准备
weather=['cloudy','cloudy','rainy','rainy','sunny','sunny']
beer=[6,8,2,4,10,12]
data=pd.DataFrame({'beer':beer,
'weather':weather})
print(data)4、数据展示
由于样本量很小，所以我们绘制箱线图，而非小提琴图。
箱线图和小提琴图的绘制的参数设置参考：
sns.boxplot(x='weather',y='beer',data=data)
5、查看均值
print(data.groupby('weather').mean())
6、常规步骤进行方差分析
excel操作步骤可参考：
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、解释变量为分类变量的正态线性模型
本例的线性模型如下：
变量“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.0rainy+4.0sunny，解释为：当weather为cloudy时，beer=7.0-4.00+4.00=7.0；当weather为rainy时，beer=7.0-4.01+4.00=3.0；当weather为sunny时，beer=7.0-40+41=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))
temp的归回效应显著大于残差效应，即存在线性回归关系。利用summary函数可以得到相关统计量的结果：
