使用Python实现方差分析

结合方差分析的内容，使用Python实现方差分析的处理代码。

01 样例说明

给出如下的样例数据：

**数据背景：**有A、B、C、D四个地区，不同地区的销售量不一样，现抽取了不同时间段内每个地区的销售量，试解决：

1、每个地区间的销售量是否相同？

2、不同月份的销售量是否相同？

3、不同时间与地区的销售量是否相同？

02 数据读取

首先使用pandas将数据读入DataFrame中。

# 使用index_col指定第一列为索引
df = pd.read_csv('/content/sales.csv',index_col=0)
df.index.name='sampledate'

03 判断每个地区间的销售量是否相同

（1）提出假设

地区对销售量没有显著影响

不全相等 地区对销售量有显著影响

（2）计算总平方和

# 计算总平方和
area_mean = df.mean()
entire_mean = df.mean().mean()
SST = ((df - entire_mean)**2).sum().sum()
SST
280581604.25

（3）计算组间平方和

SSA = 0
for col, colitems in df.iteritems():
  cnt = len(colitems)
  SSA += cnt*((colitems.mean() - entire_mean)**2)

SSA
57894573.52272727

（4）计算组内平方和

SSE = 0
for col, colitems in df.iteritems():
  col_mean = colitems.mean()
  SSE += ((colitems - col_mean)**2).sum()

SSE
222687030.72727272

可以验证 。

SST - SSA - SSE
0.0

（5）计算统计量

SST的自由度为n-1，其中n为全部观测值的个数。

SSA的自由度为k-1，其中k为因素水平（总体）的个数。

SSE的自由度为n-k。

MSA = SSA/(len(df.columns)-1)
MSA
19298191.174242426

MSE = SSE/(df.size-len(df.columns))
MSE
5567175.768181818

服从F分布的统计量为：

F = MSA/MSE
F
3.4664239064514066

（6）统计决策

若去显著性水平, F分布的临界值为

from scipy.stats import f
df1 = len(df.columns)-1
df2 = df.size-len(df.columns)

f.ppf(0.95, df1, df2)
2.8387453980206443

因为F大于，属于小概率事件，拒绝原假设，即不同地区间的销售量存在显著差异。

04 判断不同月份的销售量是否相同

（1）数据预处理

# 提取月份
df['Date'] = pd.to_datetime(df.index,format='%Y/%m/%d')
df['month'] = df['Date'].dt.month
# 按月求和并转置
df2 = df.groupby('month').sum()
df2.transpose()

month	1	3	12
A	26194	5548	24913
B	17089	7046	23494
C	25318	21568	35151
D	21181	19446	20849

（2） 计算统计量

结合上一小节的计算步骤，将单因素方差分析写成一个单一的函数，进行调用

'''单因素方差分析'''
def single_anova(dfx):
  col_mean = dfx.mean()
  entire_mean = dfx.mean().mean()
  SST = ((dfx - entire_mean)**2).sum().sum()

  SSA = 0
  for col, colitems in dfx.iteritems():
    cnt = len(colitems)
    SSA += cnt*((colitems.mean() - entire_mean)**2)

  SSE = SST - SSA

  # 计算自由度
  df1 = len(dfx.columns)-1
  df2 = dfx.size-len(dfx.columns)

  MSA = SSA/df1
  MSE = SSE/df2

  F = MSA/MSE
  print(F)
  F_edge = f.ppf(0.95, df1, df2)
  print(F_edge)

  if F>F_edge:
    print('假设不成立，各因素之间存在显著差异。')
  else:
    print('假设成立，各因素之间无显著差异。')

调用函数，直接可知结果

single_anova(df2.transpose())
4.096597704737729
4.25649472909375
假设成立，各因素之间无显著差异。

05 判断不同时间与地区的销售量是否相同

首先因为时间和地区这两个维度对销售量的影响是相互独立，因此可以采用无交互作用的双因素方差进行分析。

根据计算公式，直接给出计算函数。

'''双因素方差分析'''
def double_anova(dfy, alpha=0.05):
  col_mean = dfy.mean()
  row_mean = dfy.mean(1)

  entire_mean = dfy.mean().mean()
  print('entire_mean is {:f}'.format(entire_mean))

  # 总样本的误差平方和
  SST = ((dfy - entire_mean)**2).sum().sum()
  
  # 行因素产生的误差平方和
  SSR = ((row_mean - entire_mean)**2).sum() * len(dfy.columns)

  # 列因素产生的误差平方和
  SSC = ((col_mean - entire_mean)**2).sum() * len(dfy.index)

  # 随机误差平方和
  SSE = SST - SSR - SSC

  print('SST is {:f}'.format(SST))
  print('SSR is {:f}'.format(SSR))
  print('SSC is {:f}'.format(SSC))
  print('SSE is {:f}'.format(SSE))

  # 自由度
  f_col = len(dfy.columns) - 1
  f_row = len(dfy.index) - 1
  f_err = f_row * f_col

  # 均方差
  MSR = SSR/f_row
  MSC = SSC/f_col
  MSE = SSE/f_err

  print('MSR is {:f}'.format(MSR))
  print('MSC is {:f}'.format(MSC))
  print('MSE is {:f}'.format(MSE))

  # F统计量
  F_r = MSR/MSE
  F_c = MSC/MSE
  
  F_r_edge = f.ppf(1-alpha, f_row, f_err)
  F_c_edge = f.ppf(1-alpha, f_col, f_err)
  print(F_r_edge,F_c_edge)

  if F_r>F_r_edge:
    print('假设不成立，各行因素之间存在显著差异。')
  else:
    print('假设成立，各行因素之间无显著差异。')

  if F_c>F_c_edge:
    print('假设不成立，各列因素之间存在显著差异。')
  else:
    print('假设成立，各列因素之间无显著差异。')
  

读入数据进行方差分析。

double_anova(df)
entire_mean is 5631.750000
SST is 280581604.250000
SSR is 77678968.000000
SSC is 57894573.522727
SSE is 145008062.727273
MSR is 7767896.800000
MSC is 19298191.174242
MSE is 4833602.090909
2.164579917125473 2.9222771906450378
假设成立，各行因素之间无显著差异。
假设不成立，各列因素之间存在显著差异。

因此，不能认为时间对销售量没有显著影响，地区对销售量有显著影响。

06 总结

1. 本文按照方差分析的计算公式使用python实现了方差分析的函数，传入dataframe后可以输出分析结果；
2. 其中，双因素方差分析是SSR和SSC计算中要注意计算平方差后要乘以对应的行数或列数，这一步容易遗漏；
3. 另外，scipy.stats
   和statsmodels
   包中，有对应的ANOVA函数可以调用。