单双因素方差分析及Python实现
概述
- 方差分析 (analysis of variance, ANOVA) 在20世纪20年代发展起来的一种统计方法, 由英国统计学家费希尔在进行试验设计时为解释试验数据而首先引入.
- 方差分析就是通过检验个总体的均值是否相等来判断分类型自变量对数值型因变量是否有显著影响。
- 从表面上看, 方差分析是检验多个总体均值是否相等的方法, 但本质上是它所研究的是分类型自变量对数值型因变量的影响. 例如, 变量之间有没有关系, 关系的强度如何等.
- 根据所分析的分类型自变量多少, 可以分为单因素方差分析 (one-way analysis of variance) 和 双因素方差分析 (two-way analysis of variance) .
方差分析的基本思想和原理
- 在判断均值之间是否有差异时需要借助方差, 即通过对数据误差来源的分析来判断不同总体的均值是否相等, 进而分析自变量对因变量是否有显著影响.
- 在方差分析中, 数据的误差使用平方和来表示.
方差分析中的基本假定
- 每个总体都应该服从正态分布
- 各个总体的方差$\sigma^2$必须相同.
- 观测值是独立的.
问题一般提法
- 设因素有$k$个水平, 每个水平的均值分别用$\mu_1, \mu_2, \dotsb, \mu_k$表示, 要检验$k$个水平 (总体) 的均值是否相等, 需要提出如下假设:
- $H_0: \mu_1 = \mu_2 = \dotsb = \mu_k$ 自变量对因变量没有显著影响
- $H_1: \mu_1 , \mu_2 , \dotsb , \mu_k$不全相等 自变量对因变量没有显著影响
单因素方差分析
- 单因素方差分析研究的是一个分类型自变量对一个数值型因变量的影响.
分析步骤
提出假设
构造检验的统计量
- 计算各样本的均值
$$\bar{x_i} = \frac{\sum_{j=1}^{n_i} x_{ij}}{n_i}, i =1,2,\dotsb, k$$
式中, $n_i$为第$i$个总体的样本量; $x_{ij}$为第$i$个总体的第$j$个观测值. - 计算全部观测值的总均值
$$\bar{\bar{x}} = \frac{\sum_{i=1}^k \sum_{j=1}^{n_i} x_{ij}}{n} = \frac{\sum_{i=1}^k n_i \bar{x_i}}{n}$$
式中, $n= n_1 + n_2 + \dotsb + n_k$ - 计算各误差平方和
- 总平方和 (sum of squares for total) 记为SST. 它是全部观测值$x_{ij}$与总均值$\bar{\bar{x}}$的误差平方和, 公式为: $$SST = \sum_{i=1}^k\sum_{n=1}^{n_i} (x_{ij} - \bar{\bar{x}})^2$$
- 组间平方和 (sum of squares for factor A), 记为SSA. 他是各组均值$\bar{x_i}$与总均值$\bar{\bar{x}}$的误差平方和, 反映各样本均值之间的差异程度, 又称为因素平方和. 公式为: $$SSA = \sum_{i=1}^k n_i(\bar{x_i} - \bar{\bar{x}})^2$$
- 组内平方和 (sum of squares for error), 记为SSE. 它是每个水平或组的各个样本数据与其组均值的误差平方和, 反映每个样本个观测值得离散状况, 有称为误差平方和, 反映随机误差的大小. 公式为: $$SSE = \sum_{i=1}^j\sum_{j=1}^{n_j} (x_{ij} - \bar{x_i})^2$$
- 上述三个平方和之间关系为 $$SST = SSA + SSE$$
- 计算统计量
MSE: SSE的均方称为组内均方或组内方差. 公式为: $$MSE = \frac{SSE}{n - k}$$MSA: SSA的均方称为组间均方或组间方差. 公式为: $$MSA = \frac{SSA}{k - 1}$$检验统计量F: MSA与MSE的比值, 当$H_0$为真时, 二者的比值服从分子自由度为$k - 1$, 分母自由度为$n - k$的$F$分布, 即 $$F = \frac{MSA}{MSE} \sim F(k-1, n-k)$$
做出统计决策
- 若原假设$H_0$成立, 则表明没有系统误差, 组间方差MSA与组内方差MSE的比值差异就不会太大; 若组间方差大于组内方差, 说明各总体之间显然不仅有随机误差, 还有系统误差.
- 将统计量的值$F$与给定的显著水平$\alpha$下的临界值$F_\alpha$进行比较:
- 若$F > F_\alpha$或者$P < \alpha$, 则拒绝原假设$H_0$, 即检测因素对观测值有显著影响.
- 若$F < F_\alpha$或者$P > \alpha$, 则不拒绝原假设$H_0$, 即检测因素对观测值没有显著影响.
关系强度的测量
- 可以用组间平方和SSA占总平方和的比例大小来反映, 这一比例记为$R^2$, 公式为: $$R^2 = \frac{SSA}{SST}$$
方差中的多重比较
- 多重比较方法 (multiple comparison procedures) 通过对总体均值之间的配对比较来进一步检验到底哪些均值之间存在差异.
最小显著差异方法 (least significant different, LSD)
- 由费希尔提出, 具体步骤
- 提出假设: $H_0: \mu_1 = \mu_2; H_1: \mu_1 \ne \mu_2$
- 计算检验统计量: $\bar{x_i} - \bar{x_j}$
- 计算LSD, 公式为:
$$LSD = t_{\alpha/2}\sqrt{MSE(\frac{1}{n_i} + \frac{1}{n_j})}$$
式中, $t_{\alpha/2}$为t分布的临界值, 自由度为n-k, k为因素中总体的个数; MSE为组内方差; $n_i$和$n_j$分别是第$i$个样本和第$j$个样本的样本量. - 根据显著性水平$\alpha$做出决策
- 若$|\bar{x_i} - \bar{x_j}| > LSD$, 则拒绝$H_0$
- 若$|\bar{x_i} - \bar{x_j}| < LSD$, 则不拒绝$H_0$
python实现
import pandas as pd
import numpy as np
from scipy import stats
def f_oneway(df, f_col, v_col):
df_tmp = df[[f_col, v_col]].copy()
x_bar = df_tmp[v_col].mean()
# 计算SST
ss_total = np.sum((df_tmp[v_col] - x_bar)**2)
# 计算SSA和SSE
ss_error = 0
ss_a = 0
groups = df_tmp.groupby(f_col)
for val, group in groups:
ss_error += np.sum((group[v_col] - group[v_col].mean())**2)
ss_a += len(group) * (group[v_col].mean() - x_bar)**2
# 计算MSA和MSE
n = len(df_tmp)
k = len(groups)
ms_a = ss_a / (k - 1)
ms_error = ss_error / (n - k)
# 计算F统计量和P值
F = ms_a / ms_error
P = stats.f.sf(F, k - 1, n - k)
if P < 0.05:
print('f_col: %s, v_col: %s, P值: %s' % (f_col, v_col, P))
双因素方差分析
无交互作用的双因素方差分析
分析步骤
- 提出假设
- 对行因素提出假设:
- $H_0: \mu_1 = \mu_2 = \dotsb = \mu_k$ 行因素对因变量没有显著影响;
- $H_1: \mu_1 , \mu_2 , \dotsb , \mu_k$不全相等 行因素对因变量没有显著影响
- 对列因素提出假设:
- $H_0: \mu_1 = \mu_2 = \dotsb = \mu_k$ 列因素对因变量没有显著影响;
- $H_1: \mu_1 , \mu_2 , \dotsb , \mu_k$不全相等 列因素对因变量没有显著影响
- 对行因素提出假设:
- 构造检验统计量
- 全部观测值$x_{ij}$与总均值$\bar{\bar{x}}$的误差平方和SST, 公式为: $$SST = \sum_{i=1}^k\sum_{n=1}^{n_i} (x_{ij} - \bar{\bar{x}})^2$$
- 行因素所产生的误差平方和SSR, 公式为: $$SSR = \sum_{i=1}^k\sum_{n=1}^{n_i} (\bar{x}_{i\cdot} - \bar{\bar{x}})^2$$
- 列因素所产生的误差平方和SSC, 公式为: $$SSC = \sum_{i=1}^k\sum_{n=1}^{n_i} (\bar{x}_{\cdot j} - \bar{\bar{x}})^2$$
- 随机误差平方和SSE, 公式为: $$SSE = \sum_{i=1}^k\sum_{n=1}^{n_i} (x_{ij} - \bar{x}{i\cdot} - \bar{x}{\cdot j} - \bar{\bar{x}})^2$$
- 上述各平方和关系: $$SST = SSR + SSC + SSE$$
- 行因素的均方MSR: $$MSR = \frac{SSR}{k-1}$$
- 列因素的均方MSC: $$MSC = \frac{SSC}{r-1}$$
- 随机误差项的均方MSE: $$MSE = \frac{SSE}{(k-1)(r-1)}$$
- 为检验行因素对因变量的影响是否显著, 采用统计量$F_R$, 公式为: $$F_R = \frac{MSR}{MSE} \sim F(k-1, (k-1)(r-1))$$
- 为检验列因素对因变量的影响是否显著, 采用统计量$F_C$, 公式为: $$F_C = \frac{MSC}{MSE} \sim F(r-1, (k-1)(r-1))$$
- 做出统计决策
- 若$F_R > F_{\alpha}$, 则拒绝行原假设, 即行因素对观测值有显著影响
- 若$F_C > F_{\alpha}$, 则拒绝列原假设, 即列因素对观测值有显著影响
关系强度的测量
- 行平方和与列平方和之和度量了两个自变量对应变量的联合效应, 其与总平方和的比值$R^2$, 反映了两个自变量与因变量的关系强度, 公式为: $$R^2 = \frac{SSR + SSC}{SST}$$
有交互作用的双因素方差分析
分析步骤
- 提出假设
- 构造检验统计量
- 总平方和SST: $$SST = \sum_{i=1}^k \sum_{j=1}^r \sum_{l=1}^m (x_{ijl} - \bar{\bar{x}})^2$$
- 行平方和SSR: $$SSR = rm\sum_{i=1}^k(\bar{x}_{i\cdot} - \bar{\bar{x}})^2$$
- 列平方和SSC: $$SSC = km\sum_{j=1}^r(\bar{x}_{\cdot j} - \bar{\bar{x}})^2$$
- 交互作用平方和SSRC: $$SSRC = m\sum_{i=1}^k \sum_{j=1}^r(\bar{x}{ij} - \bar{x}{i\cdot} - \bar{x}_{\cdot j} + \bar{\bar{x}})^2$$
- 误差平方和SSE: $$SSE = SST - SSR - SSC - SSRC$$
- 行因素均方MSR: $$MSR = \frac{SSR}{k-1}$$
- 列因素均方MSC: $$MSC = \frac{SSC}{r-1}$$
- 交互作用均方MSRC: $$MSRC = \frac{SSRC}{(k-1)(r-1)}$$
- 误差均方MSE: $$MSE = \frac{SSE}{kr(m - 1)}$$
- 做出统计决策
Python双因素方差实现
- 无交互与有交互
import pandas as pd
from scipy import stats
def f_twoway(df_c, col_fac1, col_fac2, col_sta, interaction=False):
df = df_c.copy()
list_fac1 = df[col_fac1].unique()
list_fac2 = df[col_fac2].unique()
k = len(list_fac1) # 行数
r = len(list_fac2) # 列数
# 全部kkr个样本数据的总平均值
x_bar = df[col_sta].mean()
# 计算总样本的误差平方和
SST = ((df[col_sta] - x_bar)**2).sum()
list_qa = []
list_qb = []
# 计算行因素误差平方
for i in list_fac1:
series_i = df[df[col_fac1]==i][col_sta]
xi_bar = series_i.mean()
list_qa.append((xi_bar - x_bar)**2)
# 计算列因素误差平方
for j in list_fac2:
series_j = df[df[col_fac2]==j][col_sta]
xj_bar = series_j.mean()
list_qb.append((xj_bar - x_bar)**2)
df_res = pd.DataFrame(columns=['方差来源','平方和','自由度','均方','F值','P值'])
if interaction == False:
SSR = r * sum(list_qa)
SSC = k * sum(list_qb)
# 随机误差平方和
SSE = SST - SSR - SSC
MSR = SSR / (k - 1)
MSC = SSC / (r - 1)
MSE = SSE / ((r - 1) * (k - 1))
# 计算F统计量
F_r = MSR / MSE
F_c = MSC / MSE
# 计算P值
P_r = stats.f.sf(F_r, k - 1, (r - 1) * (k -1 ))
P_c = stats.f.sf(F_c, r - 1, (r - 1) * (k -1 ))
df_res['方差来源'] = [col_fac1, col_fac2, '误差','总和']
df_res['平方和'] = [SSR, SSC, SSE , SST]
df_res['自由度'] = [k - 1, r - 1, (r - 1)*(k - 1), r * k - 1]
df_res['均方'] = [MSR, MSC, MSE, '-']
df_res['F值'] = [F_r, F_c, '-', '-']
df_res['P值'] = [P_r, P_c, '-', '-']
return df_res
elif interaction == True:
list_qe = []
m = len(df[(df[col_fac1] == list_fac1[0]) & (df[col_fac2] == list_fac2[0])])
for i in list_fac1:
for j in list_fac2:
series_ij = df[(df[col_fac1]==i) & (df[col_fac2]==j)][col_sta]
list_qe.append(((series_ij - series_ij.mean())**2).sum())
# 计算各误差平方和
SSR = r * m * sum(list_qa)
SSC = k * m * sum(list_qb)
SSE = sum(list_qe)
SSRC = SST - SSR - SSC -SSE
# 计算各均方误差
MSR = SSR / (k - 1)
MSC = SSC / (r - 1)
MSRC = SSRC / ((k - 1)*(r - 1))
MSE = SSE / (k * r * (m - 1))
# 计算F统计量
F_r = MSR / MSE
F_c = MSC / MSE
F_rc = MSRC / MSE
# 计算P值
P_r = stats.f.sf(F_r, k - 1, k * r * (m - 1))
P_c =stats.f.sf(F_c, r - 1, k * r * (m - 1))
P_rc = stats.f.sf(F_rc, (k - 1)*(r - 1), k * r * (m - 1))
df_res['方差来源'] = [col_fac1, col_fac2, col_fac1 + '*' + col_fac2, '误差', '总和']
df_res['平方和'] = [SSR, SSC, SSRC, SSE, SST]
df_res['自由度']=[k - 1, r - 1, (k - 1)*(r - 1), k * r * (m - 1), k * r * m - 1]
df_res['均方']=[MSR, MSC, MSRC, MSE, '-']
df_res['F值']=[F_r, F_c, F_rc, '-', '-']
df_res['值']=[P_r, P_c, P_rc, '-', '-']
return df_res
else:
return 'interaction参数错误'