python 方差分析

Posted 2022-10-09

方差分析Analysis of Variance，ANOVA是对变量之间关系的定性分析方法，其探究的是一个因子变量对反应变量的影响。

1.分类

类别	因子数量	举例	备注
单因素方差分析	1	失业率、股票收益率
多因素方差分析	>=2	施肥量和灌溉量对于粮食产量的影响	探究每个因子的影响，而不是因子整体的影响
析因方差分析	>=2	酗酒与年龄段	因子之间也存在着影响

补充：

2.实验的影响因素

1. 不可控的随机因素 指那些无法规避的影响，如相同的土地里，相同的种子，获取的产量却有所差异。
2. 可控的实验因素 通过分析不同来源的波动（不可控随机VS 可控因素）对总波动（反应变量的总体变化）的贡献大小，从而确定可控因素（因子）对反应变量影响力的大小。

3.三个重要概念

以单因素为例， 说明方差分析的假设检验过程。

3.1 离差平方和

假设现在因子变量共有 $M$个水平，每个水平下试验或观测对象有$N_j$个$(j=1,2,...,M)。$令$Y_ij$表示第$j$个水平组别下反应变量的均值，$μ_0$代表所有反应变量的均值。若因子水平对反应变量无影响，则不同因子水平下反应变量的均值是相同的，这就是方差分析之原假设∶ $$H_0: μ_1 = μ_2 = ... = μ_M = μ_0$$ 现在我们观测到不同因子水乎下的样本数据$y_ij(j=1,2,...M,i=1,2,...,N_j)$,这样第$j$组的样本均值: $$\\baryj = \\fracy1j+y_2j+...y_N_jjN_j=\\frac1N_j\\sum_i=1^N_jy_ij, j =1,2,...,M$$ 而全样本的平均值为:

$$\\bary = \\frac1N\\sum_j=1^MN_j \\baryj$$ 其中 $N = \\sumj=1^MN_j \\baryj$ 为全样本数量。方差分析实质就是检验$\\bary是否与\\baryj相异$，现在样本观测值$yij$ 与全样本均值$\\bary$之间的偏差可以分为两个部分： $$yij - \\bary = y_ij - \\baryj + \\baryj - \\bary$$ 其中，$yij - \\baryj$ 被称为组内偏差, $\\baryj - \\bary$ 称为组内偏差。将上式两边加总，就得到反映样本数据波动情况的指标——总离差平方 Total Sum of Squares TSS $$\\sumj=1^M\\sumi=1^N_j(yij-\\bary) ^2= \\sum_j=1^M\\sum_i=1^N_j(y_ij-\\baryj)^2 + \\sumj=1^M\\sum_i=1^N_j(\\baryj-\\bary) ^2=\\sumj=1^M\\sum_i=1^N_j(y_ij-\\baryj)^2 + \\sumj=1^MN_j(\\bary_j-\\bary) ^2$$

其中等号右边第一项为误差平方和 Error Sum of Square ESS;第二项为因子平方和 Factor Sum of Squares FSS

3.2自由度

自由度是指当以样本的统计量来估计总体的参数时，样本中能够独立或自由变动的数据的个数。比如，样本方差的计算公式为∶ $$S = \\frac1n-1\\sum_i=1^n(x_i-\\barx)^2$$ 其中分母$n-1$就是自由度。

总结起来，我们进行方差分析的对象共有 $N$ 个样本观测值，分布在 $M$ 个组中，第 $j$ 个组的样本量为 $N_j$。

• $TSS$ 是衡量的是$N$个样本的总波动水平，这里所有的$N$个样本并不独立，它们满足一个约束条件（均值为$\\bary$）: $$\\sum_j=1^M\\sum_i=1^N_j(y_ij-\\bary) = 0 $$ 故真正独立的变量只有N-1个，TSS的自由度为N-1。

• $FSS$ 衡量的是由于因子水平变化导致的反应变量取值的波动。但是，$M$ 个因子组别的均值并不独立，$\\baryj, j= 1,...,M$ 满足一个约束条件∶ $$\\sumj=1^MN_j(\\baryij-\\bary) = 0 $$ 因此也丢失一个自由度，$FSS$ 的自由度是 $M-1$，其平均数组间均方差为: $$MSF=\\fracFSSM-1=\\frac1M-1\\sumj=1^MN_j(\\bary_j-\\bary)^2$$

• $ESS$ 反应的是由于样本与其所处因子水平的组别均值的偏差而产生的波动，需要满足 $M$ 个约束条件 $$\\sum_i=1^N_j(y_ij-\\baryj) = 0 , j=1,...,N$$ 从而失去了 $M$ 个自由度，所以 $ESS$ 的自由度是 $N-M$ ，其平均数组内均方差为∶ $$MSE=\\frac1N-M\\sumj=1^M\\sum_i=1^N_j(y_ij-\\bary_j)^2$$

• $TSS$、$FSS$ 和 $ESS$ 的自由度满足如下关系: $$N-1=（M-1）+（N一M）$$

3.3显著性检验

假设反应变量$Y_ij$满足条件: 根据因子水平划分的任一 $j$ 组，$Y_ij（i=1，2…，N_j）$为一 组独立同分布变量，且服从正态分布，即$Y_ij～N（μ_j, \\sigma_0^2）$。基于这个假设，可证明出组间均方差和组内均方差的期望值满足下列公式∶ $$E(MSF)=\\sigma_0^2 + \\frac1M-1\\sum_j=1MN_j(\\mu_j-\\mu_0)^2$$ $$E(MSE)=\\sigma_0^2$$ 在原假设$H_0: \\mu_1=\\mu_1=...=\\mu_M=\\mu_0$, $E(MSF)=E(MSE)=\\sigma_0^2$,而且方差分析的统计量： $$\\varphi = \\fracMSFMSE=\\fracFSS/(M-1)ESS/(N-M)$$ 服从$F（M-1，N-M）$分布。

• $\\varphi$ 统计量越大，说明组间均方差$MSF$与组内均方差$MSE$差异很大，且$MSF>MSE$，故$MSF$ 成为样本总波动的主要贡献，因子影响十分显著;
• $φ$ 统计量很小时，说明组间随机方差$MSE$是主要的方差来源，因子影响不显著。

我们可以查阅$F分布$的临界值表，或者计算$p值$来判断该统计量是否显著。

4. 方差分析的python实现

导入statsmodel模块

import pandas as pd
import statsmodels.stats.anova as anova
from statsmodels.formula.api import ols

4.1 单因素方差分析

进行分析的数据如下：

• 代码实现

model = ols(收益率 ~ C(行业), data=year_return.dropna()).fit()
table1 = anova.anova_lm(model)
print(table1)

• 输出结果

>>>                df     sum_sq    mean_sq        F          PR(>F)
>>> C(行业)    2.0  0.073032  0.036516  1.193163  0.455962
>>> Residual  2.0  0.061209  0.030604      NaN        NaN

4.2多因素方差分析

多因素方差分析的实现也很简单，我们只需要在线性回归模型里加入要研究的因素即可。现在用其来探讨婚姻状况和受教育水平对个人收入的影响，数据如下（部分）：

• 代码实现

model = ols(earnings ~ C(married)+C(educatn), data=year_return.dropna()).fit()
table2 = anova.anova_lm(model)
print(table2)

• 输出结果

>>>                df      sum_sq        mean_sq        F         PR(>F)
>>> C(married)     6.0  1.956487e+10  3.260811e+09  15.551238   9.355695e-18
>>> C(educatn)    19.0  2.082990e+11  1.096311e+10  52.284500  9.947527e-180
>>> Residual      4829.0  1.012553e+12  2.096818e+08     NaN         NaN

4.3析因方差分析

析因方差分析与多元素方差分析差不多，仅是多了一个因子的乘项。比如在上面的例子中，可以添加 married与educatn的乘项，以检验这两者对收入的影响是否与另一个因子的水平有关。

>>>                           df        sum_sq  ...          F         PR(>F)
>>> C(married)                6.0  1.956487e+10  ...  15.476346   1.163717e-17
>>> C(educatn)               19.0  2.082990e+11  ...  52.032706  1.368466e-178
>>> C(married):C(educatn)   114.0  2.202005e+10  ...   0.916761   7.246151e-01
>>> Residual               4745.0  9.997547e+11  ...        NaN            NaN