一个复杂的事物,其中每每有许多因素互相制约又互相依存。方差分析是一种经常使用的数据分析方法,其目的是经过数据分析找出对该事物有显著影响的因素、各因素之间的交互做用及显著影响因素的最佳水平等。python
本文介绍了方差分析的基础概念,详细讲解了单因素方差分析、双因素方差分析的原理,而且给出了它们的python实践代码。
api
本文大纲:数据结构
关于方差分析的基础概念app
单因素方差分析原理及python实现ide
双因素方差分析原理及python实现函数
在科学试验和生产实践中,影响某一事物的因素每每不少,好比化工生产中,像原料成分,剂量,反应温度,压力等等不少因素都会影响产品的质量,有些因素影响较大,有些影响较小, 为了使生产过程稳定, 保证优质高产, 就有必要找出对产品质量有显著影响的因素。如何找到影响因素呢?就须要试验, 方差分析就是根据试验的结果进行分析, 鉴别各个有关因素对试验结果影响程度的有效方法。而根据涉及到的因素个数的不一样, 又能够把方差分析分为单因素方差分析、多因素方差分析等。工具
下面咱们先重点研究单因素方差分析, 经过一个例子,引出方差分析中的几个概念:学习
某保险公司想了解一下某险种在不一样的地区是否有不一样的索赔额。因而他们就搜集了四个不一样地区一年的索赔额状况的记录以下表:spa
尝试判断一下, 地区这个因素是否对与索赔额产生了显著的影响?code
这个问题就是单因素方差分析的问题, 具体解决方法后面会说, 首先先由这个例子弄清楚几个概念:
试验指标:方差分析中, 把考察的试验结果称为试验指标, 上面例子里面的“索赔额”。
因素:对试验指标产生影响的缘由称为因素, 如上面的“地区”
这个类比的话, 就相似于y就是试验指标, 某个类别特征x就是因素, 类别特征x的不一样取值就是水平。那么经过方差分析, 就能够获得某个类别特征对于y的一个影响程度了吧, 这会帮助分析某个类别特征的重要性!
所谓单因素方差分析, 就是仅考虑有一个因素对试验指标的影响。假如因素有个水平, 分别在第个水平下进行屡次独立的观测, 所获得的试验指标数据以下:
注意这里的不必定同样,上面的例子。各整体间相互独立,所以咱们会有下面的模型:简单解释一下上面这个在说啥:就是第个水平的第个观测值,上面例子里面就是第个地区第次的索赔额。表示第个水平的理论均值, 后面的表示的随机偏差, 假设这个服从正态。第一个等式的意思就是某个观测值能够用某水平下的均值加一个偏差来表示。若是咱们想判断某个因素对于试验指标是否有显著影响, 很直观的就是咱们看看因素不一样的水平下试验指标的理论均值是否有显著差别, 即理论均值是否彻底相同, 若是有显著差别, 就说明不一样的水平对试验指标影响很大, 即对试验指标有显著影响。这也是方差分析的目标, 故把问题转换成了比较不一样水平下试验指标的均值差别。显著在这里的意思是差别达到的某种程度。
基于上面的分析, 咱们就能够把方差分析也当作一个检验假设的问题, 并有了原假设和备择假设:
:
那么这个假设检验的问题怎么验证呢?咱们得先分析一下, 为啥各个会有差别?从上面的模型中, 咱们能够看到, 因此第一个可能就是可能有差别, 好比, 那么很容易就大于。另外一个可能就是随机偏差的存在。在这样的启发下,咱们得找一个衡量所有之间差别的量, 就是下面这个了:
这个叫作总误差平方和,若是这个越大, 就表示之间的差别就越大。这里的表示总的观测值个数:接下来,咱们把这个平方和分解开为两部分:一部分是因为因素引发的差别, 这个叫作效应平方和, 另外一部分是因为随机偏差引发的差别,这个叫作偏差平方和
关于, 先固定一个, 此时对应的全部观测值, 他们之间的差别与每一个水平的理论平均值就没有关系了, 而是取决于随机偏差, 反应这些观察值差别程度的量其中
综合全部的水平, 就能够获得偏差平方和的公式以下:
而上面二者相减, 就会获得效应平方和:
经过上面的分析,咱们会获得如下三点结论:
(这个分解式为上面模型的方差分析)
这里是由于:
后面这部分的加和若是除以的话会服从自由度为的卡方(具体看第二篇卡方分布的定义),那么前面又一个r水平的累加,根据卡方分布的可加性可得这个东西服从的卡方
基于上面的分析,会获得一个单因素试验方差分析表:
这个表就把上面全部的分析都给总结好了。但实际使用中,咱们确定是不会手算的,而且通常也不看F的值,咱们是看p值的。
下面就用python实现一下上面的那个索赔额的例子, 看看单因素方差分析是怎么作的:
import pandas as pdimport numpy as np上面这段程序应该很容易懂, 首先前面是把数据构造出来, 而后进行一个方差的齐性检验, 这个用stats.levene函数, 这个的做用是要保证方差在每一个水平上某种程度上(显著水平)是一致的, 这时候才能进行后面的均值分析, 由于方差分析的实质是检验多个水平的均值是否有显著差别,若是各个水平的观察值方差差别太大,只检验均值之间的差别就没有意义了,因此要进行方差齐性检验。
from scipy import statsfrom statsmodels.formula.api import olsfrom statsmodels.stats.anova import anova_lm
# 这是那四个水平的索赔额的观测值A1 = [1.6, 1.61, 1.65, 1.68, 1.7, 1.7, 1.78]A2 = [1.5, 1.64, 1.4, 1.7, 1.75]A3 = [1.6, 1.55, 1.6, 1.62, 1.64, 1.60, 1.74, 1.8]A4 = [1.51, 1.52, 1.53, 1.57, 1.64, 1.6]
data = [A1, A2, A3, A4]# 方差的齐性检验w, p = stats.levene(*data)if p < 0.05: print('方差齐性假设不成立')
# 成立以后, 就能够进行单因素方差分析f, p = stats.f_oneway(*data)print(f, p) # 2.06507381767795 0.13406910483160134
后面经过stats.f_oneway函数就能够直接算出检验假设的值和值。咱们这里关注的是值, 拿值和给出的(通常是0.05)比, 若是,咱们就接受原假设,不然拒绝原假设,这个例子中是0.134,大于,故接受原假设,认为不一样的地区的索赔额没有显著差别。
因此单因素方差这块通常是懂了原理以后,用软件去分析,能看懂就算入门了。固然这个若是手算的话,思路就是须要先求,而后根据上面的公式计算,计算完了以后除以自由度而后相除获得值,而后比较和的大小,当,拒绝原假设,不然接受原假设。必定要注意这个值和值的比较标准是不一样的。由于这是两种假设检验的方法,值比较的这种是基于值法,而的那种是临界值法。 上面的例子咱们还能够进行那种单因素方差表的显示格式:首先改一下数据的格式values = A1.copy()groups = []for i in range(1, len(data)): values.extend(data[i])
for i, j in zip(range(4), data): groups.extend(np.repeat('A'+str(i+1), len(j)).tolist())
df = pd.DataFrame({'values': values, 'groups': groups})df
数据长这个样子了,也是咱们通常见到的pandas的形式:
经过下面的方式作单因素方差分析:
anova_res = anova_lm(ols('values~C(groups)', df).fit())anova_res.columns = ['自由度', '平方和', '均方', 'F值', 'P值']anova_res.index = ['因素A', '偏差']anova_res # 这种状况下看p值 >0.05 因此接受H0
结果以下:
这样就会获得单因素方差分析表的格式。固然, 为了考虑的全面些, 咱们应该评估检验的假设条件, 就是看看每一个数据是否是真的服从正态。这里就使用上一篇文章中学习到的判断数据是否是服从正态的方法了Shapiro-Wilk test(小样本状况下, 经常使用的正态检验方法):
# 数据格式张这样A1 = [1.6, 1.61, 1.65, 1.68, 1.7, 1.7, 1.78]A2 = [1.5, 1.64, 1.4, 1.7, 1.75]A3 = [1.6, 1.55, 1.6, 1.62, 1.64, 1.60, 1.74, 1.8]A4 = [1.51, 1.52, 1.53, 1.57, 1.64, 1.6]
data = [A1, A2, A3, A4]
from scipy.stats import shapiro
def normal_judge(data): stat, p = shapiro(data) if p > 0.05: return 'stat={:.3f}, p = {:.3f}, probably gaussian'.format(stat,p) else: return 'stat={:.3f}, p = {:.3f}, probably not gaussian'.format(stat,p)
for d in data: print(normal_judge(d))
结果以下:
stat=0.942, p = 0.660, probably gaussianstat=0.938, p = 0.655, probably gaussianstat=0.850, p = 0.096, probably gaussianstat=0.918, p = 0.489, probably gaussian
在不少状况下, 只考虑一个指标对观察值的影响显然是不够的, 这时就会用到多因素方差分析。双因素方差分析和多因素方差分析原理上一致, 下面给出一种两个因素之间有交互的一种形式写法做为补充。所谓双因素方差分析, 就是有两个因素做用于试验的指标, 因素有个水平, 因素有个水平. 现对因素的水平的每对组合都做次试验,也会获得一个表:
并设 这里的独立, 类比着单因素方差分析那里, 咱们就会先有下面的数学模型: 各独立 这里的表示的是第个因素第个因素下的第个观测值。是组合下的全部观测值的平均数(平均效应)。是随机偏差, 这个其实和单因素那里的理解是一个意思, 上面的单因素的那个表格放在双因素这里就至关于这里的其中一个小格子了。那么就开始引入一些新的公式, 由于既然每一个格子里面有平均, 那么每一行的格子和每一列的格子也会有平均, 总体上也会有平均, 因此下面就定义三个公式:
这个就是双因素试验方差分析的数学模型。对于这个模型, 咱们就会有三个假设检验的问题了:
因素A对于试验结果是否带来了显著影响
因素B对于试验结果是否带来了显著影响
二者的组合对于试验结果是否带来了显著影响
与单因素的状况相似, 咱们依然是采用平方和分解的方式进行验证。首先咱们得先计算四个平均值:
因素A的水平因素B的水平的平均值:
因素A的水平上的平均值:
因素B的水平平均值:
总平均值:
有了上面的平均值, 咱们就能够获得误差平方和了, 总误差平方和以下:
就获得了
这里也给出每一个平方和的自由度,的自由度, 自由度是, 自由度是, 自由度, 自由度是。那么和单因素水平分析那样, 咱们能够获得每一个假设下面的拒绝域形式:
导入此次用到的包(依然是单因素分析时的ols和anova_lm)
import pandas as pdimport numpy as np
from scipy import statsfrom statsmodels.formula.api import olsfrom statsmodels.stats.anova import anova_lm
# 这三个交互效果的可视化画图from statsmodels.graphics.api import interaction_plotimport matplotlib.pyplot as pltfrom pylab import mpl # 显示中文
# 这个看某个因素各个水平之间的差别from statsmodels.stats.multicomp import pairwise_tukeyhsd
3.一、无交互做用的状况
因为不考虑交互做用的影响,对每个因素组合 只需进行一次独立试验,称为无重复试验。数据:考虑三种不一样形式的广告和五种不一样的价格对某种商品销量的影响。选取某市15家大超市,每家超市选用其中的一个组合,统计出一个月的销量以下(设显著性水平为0.05):
下面进行双因素方差分析,简要流程是,先用pandas库的DataFrame数据结构来构造输入数据格式。而后用statsmodels库中的ols函数获得最小二乘线性回归模型。最后用statsmodels库中的anova_lm函数进行方差分析。
dic_t2=[{'广告':'A1','价格':'B1','销量':276},{'广告':'A1','价格':'B2','销量':352}, {'广告':'A1','价格':'B3','销量':178},{'广告':'A1','价格':'B4','销量':295}, {'广告':'A1','价格':'B5','销量':273},{'广告':'A2','价格':'B1','销量':114}, {'广告':'A2','价格':'B2','销量':176},{'广告':'A2','价格':'B3','销量':102}, {'广告':'A2','价格':'B4','销量':155},{'广告':'A2','价格':'B5','销量':128}, {'广告':'A3','价格':'B1','销量':364},{'广告':'A3','价格':'B2','销量':547}, {'广告':'A3','价格':'B3','销量':288},{'广告':'A3','价格':'B4','销量':392}, {'广告':'A3','价格':'B5','销量':378}]df_t2=pd.DataFrame(dic_t2,columns=['广告','价格','销量'])df_t2
数据长这样:
# 方差分析price_lm = ols('销量~C(广告)+C(价格)', data=df_t2).fit()table = sm.stats.anova_lm(price_lm, typ=2)table
结果以下:
能够发现这里的p值都是小于0.05的, 因此咱们要拒绝掉原假设, 便可认为不一样的广告形式, 不一样的价格均形成商品销量的显著差别。下面还能够看一下交互影响效果:
fig = interaction_plot(df_t2['广告'],df_t2['价格'], df_t2['销量'], ylabel='销量', xlabel='广告')
结果以下:
再来分析一下单因素各个水平之间的显著差别:
# 广告与销量的影响 注意这个的显著水平是0.01print(pairwise_tukeyhsd(df_t2['销量'], df_t2['广告'], alpha=0.01)) # 第一个必须是销量, 也就是咱们的指标
结果以下:
这个能够获得的结论是在显著水平0.01的时候, A2和A3的p值小于0.01, reject=True, 即认为A2和A3有显著性差别。
3.二、有交互做用的状况
因为因素有交互做用,须要对每个因素组合 分别进行 次 重复试验,称这种试验为等重复试验。
数据:几率论课本上的那个例子, 火箭的射程与燃料的种类和推动器的型号有关,现对四种不一样的燃料与三种不一样型号的推动器进行试验,每种组合各发射火箭两次,测得火箭的射程结果以下(设显著性水平为0.01):第一步依然是先构造数据,
dic_t3=[{'燃料':'A1','推动器':'B1','射程':58.2},{'燃料':'A1','推动器':'B1','射程':52.6}, {'燃料':'A1','推动器':'B2','射程':56.2},{'燃料':'A1','推动器':'B2','射程':41.2}, {'燃料':'A1','推动器':'B3','射程':65.3},{'燃料':'A1','推动器':'B3','射程':60.8}, {'燃料':'A2','推动器':'B1','射程':49.1},{'燃料':'A2','推动器':'B1','射程':42.8}, {'燃料':'A2','推动器':'B2','射程':54.1},{'燃料':'A2','推动器':'B2','射程':50.5}, {'燃料':'A2','推动器':'B3','射程':51.6},{'燃料':'A2','推动器':'B3','射程':48.4}, {'燃料':'A3','推动器':'B1','射程':60.1},{'燃料':'A3','推动器':'B1','射程':58.3}, {'燃料':'A3','推动器':'B2','射程':70.9},{'燃料':'A3','推动器':'B2','射程':73.2}, {'燃料':'A3','推动器':'B3','射程':39.2},{'燃料':'A3','推动器':'B3','射程':40.7}, {'燃料':'A4','推动器':'B1','射程':75.8},{'燃料':'A4','推动器':'B1','射程':71.5}, {'燃料':'A4','推动器':'B2','射程':58.2},{'燃料':'A4','推动器':'B2','射程':51.0}, {'燃料':'A4','推动器':'B3','射程':48.7},{'燃料':'A4','推动器':'B3','射程':41.4},]df_t3=pd.DataFrame(dic_t3,columns=['燃料','推动器','射程'])df_t3.head()
结果这样:
下面是方差分析:
moore_lm = ols('射程~燃料+推动器+燃料:推动器', data=df_t3).fit()table = sm.stats.anova_lm(moore_lm, typ=1)table结果以下:
这里获得的结论就是燃料的P值是大于0.01的, 而推动器和二者组合的p值都小于0.01, 而且二者的组合很是小, 这就说明燃料对于火箭的射程没有显著影响, 然后二者都有显著影响,二者的交互做用更是高度显著。下面是交互效应效果:
fig = interaction_plot(df_t3['燃料'],df_t3['推动器'], df_t3['射程'], ylabel='射程', xlabel='燃料')
结果以下:
从这个图里面能够看出, (A4, B1)和(A3, B2)组合的进程最好。黄金搭档。单因素差别性分析:
print(pairwise_tukeyhsd(df_t3['射程'], df_t3['燃料']))
结果:
都是False, 说明A因素各个水平之间无显著差别。
两个实验到这里就结束了, 这里再补充两点别的知识:
1. ols函数里面公式的写法
'射程~C(燃料)+C(推动器)+C(燃料):C(推动器)' :至关于射程是y(指标), 燃料和推动器是x(影响因素), 三项加和的前两项表示两个主效应, 第三项表示考虑二者的交互效应, 不加C也可。
'射程~C(燃料, Sum)*C(推动器, Sum)'和上面效果是一致的, 星号在这里表示既考虑主效应也考虑交互效应*'销量~C(广告)+C(价格)':这个表示不考虑交互相应
可是要注意, 考虑交互相应和不考虑交互相应致使的Se(残差项)会不一样, 因此会影响最终的结果。
stats.anova_lm(moore_lm, typ=1)这里面的typ参数, 这个参数我尝试尚未彻底搞明白究竟是什么意思, 这个参数有1,2,3 三个可选项, 分别表明着不一样的误差平方和的计算方法, 我在第二个实验中尝试过改这个参数,改为1的时候发现就加了一列mean_sq, 而后其余的没变。
改为3的时候发现加一行Intercept, 而且此时燃料和推动器的数据都发生了变化。
方差分析这块到这里就结束了, 随着这篇文章的结束也意味着几率统计的知识串联也到了尾声, 简单的回顾一下本篇的内容, 这篇文章主要是在实践的角度进行的分析, 方差分析在统计中仍是很经常使用的, 比较适合类别因素对于数值指标的影响程度:
首先从单因素方差分析入手, 这个只考虑了一个因素对于指标的影响, 先分析了原理,而后基于python进行了实现。实际应用中,通常是会点原理,而后使用工具实现方差分析,会看结果,这样就算入门了。
实际应用中, 或许能够经过这种方法去分析类别特征的重要性或者关联性,以及类别和类别特征之间的交互做用等。