python进行方差分析

  方差分析(Analysis of Variance, ANOVA)用于两个及两个以上样本均值差别的显著性检验。 方差分析的原理通俗的解释就是将试验数据的总离散分解为来源于不同因素(一种是不可控的随机因素；一种是研究中施加的会对结果形成影响的可控因素)的离散，并作出数据估计，从而发现各个因素在总离散中所占的重要程度。方差分析中的基本概念主要有：

• 因变量: 又称响应变量。一般指待分析的指标。
• 因素: 又称因子。一般指调查或试验中需要考察的、可以控制的条件或影响因素。因素所处的不同状态称为水平。方差分析的目的就是为了验证因素的不同水平对试验结果造成的影响是否显著。

  每个因素每一个水平之下调查结果或实验观察结果可以称为一个“组”，可以计算各“组”内部变量的均值于方差，还可以计算因素之下全部观察结果的总平均以及组与组之间的方差。方差分析的基本思想就是从不同角度计算出有关的均值和方差，然后通过组内方差与组间方差的对比，在一定统计理论指导下分析条件误差和随机误差，进而分解或判断出调查或实验观察数据中必然因素(因子)和偶然因素(随机)的影响大小。
  按影响分析指标的因素个数多少的不同，可分为单因素分析、双因素分析和多因素分析。这里只介绍单因素方差分析和双因素方差分析。

1 单因素方差分析

1.1 基本模型

  单因素方差分析只考虑一个因素$A$对因变量的影响，设因素$A$的$r$个水平记为：$A_1,A_2,\dots ,A_r$,在水平$A_i$下进行$n_i$次重复试验，可获得观察指标的$n_i$个数据：$x_{i1},x_{i2},\dots ,x_{i{n}_i}(i=1,2,\dots ,r)$。第$A_i(i=1,2,\dots ,r)$水平之下各观察指可视作来自总体$X_i(i=1,2,\dots ,r)$，且$X_i$服从正态分布，$X_i$与$X_j(i\ne j)$相互独立(等方差)。即$X_i\sim N({\overline{X}}_i,S^2)$。$X_{i1},X_{i2},\dots ,X_{i{n}_i}$表示从总体$X_i$中抽取的样本。单因素方差分析的基本统计模型就是把因变量观察值$X_{ij}$分解为三项内容之和，如下：
$$\{\begin{array}{l}{X}_{ij}=\overline{X}+{\alpha }_i+e_{ij},(i=1,2,\dots ,r;j=1,2,\dots ,n_i),e_{i}j\sim N(0,S^2)\\ {\sum }_{i=1}^r{n}_i{\alpha }_i=0\end{array}$$其中$\overline{X}=\frac{1}{n}{\sum }_{i=1}^r{n}_i{X}_i(n={\sum }_{i=1}^r{n}_i)$代表一般水平; ${\alpha }_i$表示因素$A$的不同水平导致的因变量差异，即$a_i={\overline{X}}_i-\overline{X}(i=1,2,\dots ,r)$；$e_{ij}=X_{ij}-{\overline{X}}_i$极为随机误差，通常假设其服从均值为0方差为$S^2$的正态分布。要说明因素$A$的不同水平对因变量有影响，只需证明$a_i$不全为0即可。

1.2 偏差平方和分解公式

  上文中提到的$X$、${\overline{X}}_i$、$S^2$、${\alpha }_i$和$e_{ij}$可由样本观察数据$\overline{x}$、${\overline{x}}_i$、${\alpha }_i$和$e_{ij}$来表示，即$$\overline{x}=\sum _{i=1}^r\sum _{j=1}^{n_i}{x}_{ij}/\sum _{i=1}^r{n}_i，{\overline{x}}_i=\frac{1}{n_i}\sum _{j=1}^{n_i}{x}_{ij},{\alpha }_i={\overline{x}}_i-\overline{x},e_{i}j=x_{ij}-{\overline{x}}_i，$$于是就有如下公式成立：$$\sum _{i=1}^r\sum _{j=1}^{n_i}(x_{ij}-\overline{x}{)}^2=\sum _{i=1}^r{n}_i({\overline{x}}_i-\overline{x}{)}^2+\sum _{i=1}^r\sum _{j=1}^{n_i}(x_{ij}-{\overline{x}}_i{)}^2$$上式可简记为$S{S}_T=S{S}_A+S{S}_E$。其中$S{S}_T={\sum }_{i=1}^r{\sum }_{j=1}^{n_i}(x_{ij}-\overline{x}{)}^2$是所有数据到总样本均值的距离平方和，是因变量对中心位置的变化的总度量；$S{S}_A={\sum }_{i=1}^r{n}_i({\overline{x}}_i-\overline{x}{)}^2$是组内样本均值与总平均值的离差加权平方和，是反映因子的不同水平造成因变量变化的总度量；$S{S}_E={\sum }_{i=1}^r{\sum }_{j=1}^{n_i}(x_{ij}-{\overline{x}}_i{)}^2$为组内离差平方和，是随机误差造成试验指标变化的总度量。

1.3 显著性检验

  单因子方差分析实质上是多个正态总体差异性统计检验，所以其原假设为：
$$H_0:{\alpha }_1={\alpha }_2=\cdots ={\alpha }_r=0$$上述原假设等价于：$$H_0:{\overline{X}}_1={\overline{X}}_2=\cdots ={\overline{X}}_r$$因为当$H_0$成立时，有$$\frac{S{S}_T}{S^2}=\sum _{i=1}^r\sum _{j=1}^{n_i}(\frac{x_{ij}-\overline{x}}{S}{)}^2\sim {\chi }^2(n-1)$$$$\frac{S{S}_A}{S^2}=\sum _{i=1}^r{n}_i(\frac{{\overline{x}}_i-\overline{x}}{S^2}{)}^2\sim {\chi }^2(r-1)$$$$\frac{S{S}_E}{S^2}=\frac{{\sum }_{i=1}^r{\sum }_{j=1}^{n_j}(x_{ij}-{\overline{x}}_j{)}^2}{S^2}\sim {\chi }^2(n-r)$$其中$n={\sum }_{i=1}^r{n}_i$。
  当$S{S}_T$固定不变时，如果不同因素水平之下因变量差异越大，即$S{S}_A$越大时$S{S}_E$越小。所以若$S{S}_A$显著大于$S{S}_E$，则说明各总体${\overline{X}}_i$之间的差异越显著，那么$H_0$不成立。所以使用$S{S}_A/S{S}_E$来构造检验统计量。具体如下：$$F=\frac{S{S}_A/(r-1)}{S{S}_E/(n-r)}\sim F(r-1,n-r)$$在给定显著性水平$\alpha$下，查$F$分布表值$F_{1-\alpha }(r-1,n-r)$。若统计量$F$值超过这一临界点，则拒绝原假设，认为样本均值不完全相等，否则不能拒绝原假设。

2 多因素方差分析

  设由$A,B$两个因素影响试验结果指标。因素$A$有$r$个水平，因素$B$有$s$个水平。依据因素$A,B$之间是否有相互作用，双因素方差分析又可以分为：无交互作用的双因素方法分析和有相互作用的双因素方法分析。

2.1 无交互作用的双因素方差分析

  当因素$A,B$之间没有交互作用时，因素$A$,$B$的不同水平的每种组合都只做一次试验。假设$x_{ij}(i=1,2,\dots ,r;j=1,2,\dots ,s)$之间相互独立，且$x_{ij}\sim N({\overline{X}}_{ij},S^2)$,则$x_{ij}={\overline{X}}_{ij}+e_{ij}(i=1,2,\dots ,r;j=1,2,\dots ,s)$，其中$e_{ij}\sim N(0,S^2)$。此时双因素方法分析的基本模型如下：$$\{\begin{array}{l}{X}_{ij}=\overline{X}+{\alpha }_i+{\beta }_j+e_{ij}\\ {\sum }_{i=1}^r{\alpha }_i=0,{\sum }_{j=1}^s{\beta }_j=0\end{array}$$其中$\overline{X}=1/rs{\sum }_{i=1}^r{\sum }_{j=1}^s{\overline{X}}_{ij}$为总平均值，${\alpha }_i={\overline{X}}_{i\bullet }-\overline{X}，X_{i\bullet }=1/s{\sum }_{j=1}^s{\overline{X}}_{ij}$为因素$A$在水平$i$下对因变量的效应值，${\beta }_j={\overline{X}}_{\bullet j}-\overline{X},{\overline{X}}_{\bullet j}=1/r{\sum }_{i=1}^r{\overline{X}}_{ij}$为因素$B$在水平$j$下对因变量的效应值。以上各项可有样本数据代替，具体如下：
$$\overline{x}=\frac{1}{rs}\sum _{i=1}^r\sum _{j=1}^s{x}_{ij},x_{i\bullet }=\frac{1}{s}\sum _{j=1}^s{x}_{ij},x_{\bullet j}=\frac{1}{r}\sum _{i=1}^r{x}_{ij}$$则有如下公式成立：$$\sum _{i=1}^r\sum _{j=1}^s(x_{ij}-\overline{x}{)}^2=\sum _{i=1}^r\sum _{j=1}^s({\overline{x}}_{i\bullet }-\overline{x}{)}^2+\sum _{i=1}^r\sum _{j=1}^s({\overline{x}}_{\bullet j}-\overline{x}{)}^2+\sum _{i=1}^r\sum _{j=1}^s(x_{ij}+\overline{x}-{\overline{x}}_{i\bullet }-{\overline{x}}_{\bullet j}{)}^2$$上述公式可以简记为$S{S}_T=S{S}_A+S{S}_B+S{S}_E$。
  在总体分布的正态等方差性条件满足时，有：
$$\frac{S{S}_T}{S^2}\sim {\chi }^2(rs-1),\frac{S{S}_A}{S^2}\sim {\chi }^2(r-1),\frac{S{S}_B}{S^2}\sim {\chi }^2(s-1),\frac{S{S}_E}{S^2}\sim {\chi }^2((r-1)(s-1))$$  此时方差分析的主要任务是：分析因素$A$和因素$B$对试验观察指标的影响大小，这一问题等价于以下两种假设：
(1) “因素$A$对因变量影响不显著”其对应原假设等价为：
$$H_{01}:{\overline{X}}_{1\bullet }={\overline{X}}_{2\bullet }=\cdots ={\overline{X}}_{r\bullet }$$其相应的检验统计量为：$$F_A=\frac{S{S}_A/(r-1)}{S{S}_E/(rs-r-s+1)}\sim F(r-1,rs-r-s+1)$$(2) “因素$B$对因变量影响不显著”其对应原假设等价为：
$$H_{02}:{\overline{X}}_{\bullet 1}={\overline{X}}_{\bullet 2}=\cdots ={\overline{X}}_{\bullet s}$$其相应的检验统计量为：$$F_B=\frac{S{S}_B/(s-1)}{S{S}_E/(rs-r-s+1)}\sim F(s-1,rs-r-s+1)$$

2.2 有交互作用的双因素方差分析

  当因素$A,B$之间存在相互作用时，两个因素的不同水平之下的组合都应该有若干个样本观察值。假设在每一对的因素水平组合$(A_i,B_j)$中，样本容量相同均为$n$。$x_{ijk}$表示因素$A、B$在第$i(i=1,2,\dots ,r)，j(j=1,2,\dots ,s)$个水平状态下第$k(k=1,2,\dots ,n)$个样本观测值。假设$x_{ijk}\sim N({\overline{X}}_{ij},S^2)(k=1,2,\dots ,n)$则此时的基本模型为：$$\{\begin{array}{l}{x}_{ijk}=\overline{X}+{\alpha }_i+{\beta }_j+{\gamma }_{ij}+e_{ij}\\ {\sum }_{i=1}^r{\alpha }_i=0,{\sum }_{j=1}^s{\beta }_j=0,{\sum }_{i=1}^r{\gamma }_{ij}=0,{\sum }_{j=1}^s{\gamma }_{ij}=0\end{array}$$表现在样本上则有：$$\begin{gathered} \sum _{i=1}^r\sum _{j=1}^s\sum _{k=1}^n(x_{ijk}-\overline{x}{)}^2=ns\sum _{i=1}^r({\overline{x}}_{i\bullet \bullet }-\overline{x}{)}^2+nr\sum _{j=1}^s({\overline{x}}_{\bullet j\bullet }-\overline{x}{)}^2+n\sum _{i=1}^r\sum _{j=1}^s({\overline{x}}_{ij\bullet }+\overline{x}-{\overline{x}}_{i\bullet \bullet }-{\overline{x}}_{\bullet j\bullet }{)}^2 \\ +\sum _{i=1}^r\sum _{j=1}^s\sum _{k=1}^n(x_{ijk}-{\overline{x}}_{ij\bullet }{)}^2 \end{gathered}$$其中:$$\overline{x}=\frac{1}{rsn}\sum _{i=1}^r\sum _{j=1}^s\sum _{k=1}^n{x}_{ijk},{\overline{x}}_{i\bullet \bullet }=\frac{1}{ns}\sum _{j=1}^s\sum _{k=1}^n{x}_{ijk}(i=1,2,\dots ,r)$$$${\overline{x}}_{\bullet j\bullet }=\frac{1}{rn}\sum _{i=1}^r\sum _{k=1}^n{x}_{ijk}(j=1,2,\dots ,s),{\overline{x}}_{ij\bullet }=\frac{1}{n}\sum _{k=1}^n{x}_{ijk}(i=1,2,\dots ,r;j=1,2,\dots ,s)$$上述公式可以简记为$S{S}_T=S{S}_A+S{S}_B+S{S}_{A\times B}+S{S}_E$。在总体分布的正态等方差性条件满足时，有$$\frac{S{S}_T}{S^2}\sim {\chi }^2(rns-1),\frac{S{S}_A}{S^2}\sim {\chi }^2(r-1),\frac{S{S}_B}{S^2}\sim {\chi }^2(s-1)$$$$\frac{S{S}_{A\times B}}{S^2}\sim {\chi }^2((r-1)(s-1)),\frac{S{S}_E}{S^2}\sim {\chi }^2(rsn-rs)$$  此时方差分析的主要任务是：分析因素$A$和因素$B$对因变量指标的影响大小，以及因素$A$和因素$B$的交互作用对因变量影响大小。这等价于以下三种统计假设。
(1)“因素A对因变量影响不显著”的原假设等价于：$$H_{01}:{\alpha }_1={\alpha }_2=\cdots ={\alpha }_r=0$$此时其检验统计量为： $$F_A=\frac{S{S}_A/(r-1)}{S{S}_E/(rsn-rs)}\sim F(r-1,rsn-rs)$$(2)“因素B对因变量影响不显著”的原假设等价于：$$H_{02}:{\beta }_1={\beta }_2=\cdots ={\beta }_s=0$$此时其检验统计量为：$$F_B=\frac{S{S}_B/(s-1)}{S{S}_E/(rsn-rs)}\sim F(s-1,rsn-rs)$$(3)“因素A和因素B交互作用对因变量影响不显著”的原假设等价于：$$H_{03}:{\gamma }_{ij}=0(i=1,2,\dots ,r;j=1,2,\dots ,s)$$此时其检验统计量为：$$F_{A\times B}=\frac{S{S}_{A\times B}/(rs-r-s+1)}{S{S}_E/(rsn-rs)}\sim F(rs-r-s+1,rsn-rs)$$

3.实验

import pandas as pd
import numpy as np
from scipy import stats

#1. 单因素方差分析
#背景：产品包装方式是否会对产品销量产生影响
data=pd.DataFrame(np.array([[152,208,182],
                           [188,256,198],
                           [238,300,268],
                           [192,280,220],
                           [180,270,200],
                           [115,210,128],
                           [125,185,110],
                           [100,165,105]]),columns=['A1','A2','A3'])

data_mean=data.sum().sum()/(data.shape[0]*data.shape[1])
SSA=(data.count()*(data.mean()-data_mean)**2).sum()
SSE=sum([((data[col]-data[col].mean())**2).sum() for col in data.columns])

SSA_degree=data.shape[1]-1
SSE_degree=data.count().sum()-data.shape[1]
F_stat=(SSA/SSA_degree)/(SSE/SSE_degree)
p_val=stats.f.sf(F_stat,SSA_degree,SSE_degree)
#p_val=0.02<0.05 拒绝原假设。则认为包装会对销量产生影响。

#2. 无交互作用的双因素方差分析
data_1=pd.DataFrame(np.array([[4152,6852,9852],
                              [4852,5112,8912],
                              [3546,5786,9978],
                              [5456,6105,9105],
                              [3426,3998,15918],
                              [6124,10124,16124],
                              [5846,9789,10100]]),columns=['A1','A2','A3'],
                             index=['B1','B2','B3','B4','B5','B6','B7'])

data_mean_1=data_1.sum().sum()/(data_1.shape[0]*data_1.shape[1])
SSA=data_1.shape[0]*sum((data_1.mean()-data_mean_1)**2)
SSA_degree=data_1.shape[1]-1
SSB=data_1.shape[1]*sum((data_1.mean(axis=1)-data_mean_1)**2)
SSB_degree=data_1.shape[0]-1
SSE=sum([(data_1.loc[idx,col]+data_mean_1-data_1.loc[idx,:].mean()-data_1.loc[:,col].mean())**2 
         for idx in data_1.index for col in data_1.columns])
SSE_degree=data_1.shape[0]*data_1.shape[1]-data_1.shape[1]-data_1.shape[0]+1

F_A_stat=(SSA/SSA_degree)/(SSE/SSE_degree)
p_A=stats.f.sf(F_A_stat,SSA_degree,SSE_degree)

F_B_stat=(SSB/SSB_degree)/(SSE/SSE_degree)
p_B=stats.f.sf(F_B_stat,SSB_degree,SSE_degree)

#3. 有交互作用的双因素方差分析
data_2=pd.DataFrame(np.array([[4152,4852,3546,5456,3426,6124,5846],
                              [3968,4568,5541,4879,3895,6123,5680],
                              [4944,4795,5687,4598,3678,6856,5982],
                              [4852,4132,4752,5123,4894,7102,6892],
                              [4454,5423,5124,6811,3987,8951,8123],
                              [5914,5684,6412,5132,4102,9102,8741],
                              [4745,3654,5612,5456,4243,9581,8210],
                              [5981,3986,5781,4925,4410,9451,8652],
                              [6852,5112,5786,6105,3998,10124,9789],
                              [7167,5418,8841,5912,4213,12563,11752],
                              [5144,5958,7682,5109,4318,10986,11012],
                              [4888,5702,5556,6200,5209,12107,13785],
                              [5784,6582,5879,6941,4125,13958,14843],
                              [6584,6511,6987,6100,4942,14129,12874],
                              [6475,6845,6555,7812,5109,13891,13290],
                              [7581,6124,7001,5261,5097,14121,14589],
                              [9852,8912,9978,9105,15918,16124,10100],
                              [10165,9741,11049,11591,15983,17569,10222],
                              [9142,9898,10054,9751,15719,14902,9055],
                              [9888,10095,9941,8912,15303,16139,11019],
                              [9912,11081,9875,8754,14789,13988,11980],
                              [10109,10055,9002,8288,14105,14373,11274],
                              [11075,11286,9915,10112,15666,14590,12290],
                              [12781,11029,9989,9123,16444,16879,11589]]),
                    columns=['B1','B2','B3','B4','B5','B6','B7'],
                    index=pd.MultiIndex.from_product([['A1','A2','A3'],
                                                      [1,2,3,4,5,6,7,8]]))
data_mean_2=data_2.sum().sum()/data_2.size

SSA=0
A_level=set(data_2.index.get_level_values(0))
for idx in A_level:
    tmp=data_2.loc[idx]
    tmp_mean=tmp.sum().sum()/tmp.size
    SSA+=(tmp.size*((tmp_mean-data_mean_2)**2))
SSA_degree=len(A_level)-1

SSB=data_2.shape[0]*((data_2.mean()-data_mean_2)**2).sum()
SSB_degree=data_2.shape[1]-1

SSAB=0
for idx in A_level:
    for col in data_2.columns:
        tmp_mean=data_2.loc[idx,col].mean()
        idx_mean=(data_2.loc[idx].sum().sum())/(data_2.loc[idx].size)
        col_mean=data_2[col].mean()
        SSAB+=data_2.loc[idx].shape[0]*(tmp_mean+data_mean_2-idx_mean-col_mean)**2
SSAB_degree=(len(A_level)-1)*(data_2.shape[1]-1)

SSE=0
for idx in A_level:
    for col in data_2.columns:
        tmp=data_2.loc[idx,col]
        SSE+=((tmp-tmp.mean())**2).sum()
SSE_degree=data_2.size-len(A_level)*data_2.shape[1]
        
F_A_stat=(SSA/SSA_degree)/(SSE/SSE_degree)
p_A=stats.f.sf(F_A_stat,SSA_degree,SSE_degree)
        
F_B_stat=(SSB/SSB_degree)/(SSE/SSE_degree)
p_B=stats.f.sf(F_B_stat,SSB_degree,SSE_degree)

F_AB_stat=(SSAB/SSAB_degree)/(SSE/SSE_degree)
p_AB=stats.f.sf(F_AB_stat,SSAB_degree,SSE_degree)

参考资料

1. https://zhuanlan.zhihu.com/p/49479042
2. 《统计学》