方差分析
一、基本概念
方差分析(Analysis of variance, ANOVA) 主要研究分类变量作为自变量时，对因变量的影响是否是显著的。
方差分析的方法是由20世纪的统计学家Ronald Aylmer Fisher在1918年到1925年之间提出并陆续完善起来的，该方法刚开始是用于解决田间实验的数据分析问题，因此，方差分析的学习是和实验设计、实验数据的分析密不可分的。
实验设计和方差分析都有自己相应的语言。因此，在这里我们通过一个焦虑症治疗的实例，先了解一些术语，并且思考一下，方差分析主要用于解决什么样的问题。
以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR，设计方案如表1-1所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。
| CBT | EMDR |
|---|---|
| s1 | s6 |
| s2 | s7 |
| s3 | s8 |
| s4 | s9 |
| s5 | s10 |
在这个实验设计中，治疗方案是两水平（CBT、EMDR）的组间因子。之所以称其为组间因子，是因为每位患者都仅被分配到一个组别中，没有患者同时接受CBT和EMDR。表中字母s代表受试者（患者）。STAI是因变量，治疗方案是自变量。由于在每种治疗方案下观测数相等，因此这种设计也称为均衡设计（balanced design）；若观测数不同，则称作非均衡设计（unbalanced design）。
因为仅有一个类别型变量，表1的统计设计又称为单因素方差分析（one-way ANOVA），或进一步称为单因素组间方差分析。方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明五周后两种疗法的STAI得分均值不同。
假设你只对CBT的效果感兴趣，则需将10个患者都放在CBT组中，然后在治疗五周和六个月后分别评价疗效，设计方案如表1-2所示。
| 时间 | ||
|---|---|---|
| 患者 | 5周 | 6个月 |
| s1 | ||
| s2 | ||
| s3 | ||
| s4 | ||
| s5 | ||
| s6 | ||
| s7 | ||
| s8 | ||
| s9 | ||
| s10 |
此时，时间（time）是两水平（五周、六个月）的组内因子。因为每位患者在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析；又由于每个受试者都不止一次被测量，也称作重复测量方差分析。当时间的F检验显著时，说明患者的STAI得分均值在五周和六个月间发生了改变。
现假设你对治疗方案差异和它随时间的改变都感兴趣，则将两个设计结合起来即可：随机分配五位患者到CBT，另外五位到EMDR，在五周和六个月后分别评价他们的STAI结果（见表1-3)。
| 时间 | |||
|---|---|---|---|
| 疗法 | 患者 | 5周 | 6个月 |
| CBT | s1 | ||
| s2 | |||
| s3 | |||
| s4 | |||
| s5 | |||
| EMDR | s6 | ||
| s7 | |||
| s8 | |||
| s9 | |||
| s10 |
疗法（therapy）和时间（time）都作为因子时，我们既可分析疗法的影响（时间跨度上的平均）和时间的影响（疗法类型跨度上的平均），又可分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。
当设计包含两个甚至更多的因子时，便是因素方差分析设计，比如两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。若因子设计包括组内和组间因子，又称作混合模型方差分析，当前的例子就是典型的双因素混合模型方差分析。
本例中，你将做三次F检验：疗法因素一次，时间因素一次，两者交互因素一次。若疗法结果显著，说明CBT和EMDR对焦虑症的治疗效果不同；若时间结果显著，说明焦虑度从五周到六个月发生了变化；若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五周到六个月的改变程度在两种疗法间是不同的）。
现在，我们对上面的实验设计稍微做些扩展。众所周知，抑郁症对病症治疗有影响，而且抑郁症和焦虑症常常同时出现。即使受试者被随机分配到不同的治疗方案中，在研究开始时，两组疗法中的患者抑郁水平就可能不同，任何治疗后的差异都有可能是最初的抑郁水平不同导致的，而不是由于实验的操作问题。抑郁症也可以解释因变量的组间差异，因此它常称为混淆因素（confounding factor）。由于你对抑郁症不感兴趣，它也被称作干扰变数（nuisance variable）。
假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。本案例中，BDI为协变量，该设计为协方差分析（ANCOVA）。
以上设计只记录了单个因变量情况（STAI），为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价）。当因变量不止一个时，设计被称作多元方差分析（MANOVA）， 若协变量也存在， 那么就叫多元协方差分析（MANCOVA）。
二、实例
这里采用Pokemon的种族值数据进行方差分析，数据源请移步https://github.com/datawhalechina/joyful-pandas/tree/master/data
数据结构如下：
表中包含了1~6代的Pokemon种族值数据，首先我们想确定每一代Pokemon的种族值总和（Total）是否有显著差异。
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.multicomp import MultiComparison
df=pd.read_csv('E:\pandas\source2\joyful-pandas-master\data\Pokemon.csv',encoding='ANSI')
##统计代数对宝可梦种族值是否有显著影响
df1=df.drop(index=(df.loc[(df['Generation']=='FALSE')].index)) #数据清洗，删除Generation有误的数据
g=df1['Generation'].unique() #统计总共有几代
args=[]
for i in list(g): #将数据按代数进行分类
args.append(df1[df1['Generation']==i]['Total'])
w,p=stats.levene(*args) #方差齐次性检验，方差齐次才能进行方差分析
if p<=0.05:
print('方差齐次性不成立')
else:
print('方差齐次，可进行下一步操作')
f,p1=stats.f_oneway(*args) #方差分析
if p1<=0.05:
print('代数对种族值有影响')
else:
print('没啥影响')
输出结果为：
方差齐次，可进行下一步操作
没啥影响
可以得出，代数对种族值没有影响，说明每代Pokemon的种族值之间没有显著差异。看来这款游戏的平衡做的还不错，没有出现一代更比一代强的情况。。
接下来看看属性对种族值有没有影响，这里只考虑第一属性。（Type 1）
#统计属性对宝可梦种族值是否有显著影响
df1=df.drop(index=(df.loc[(df['Generation']=='FALSE')].index))
g=df1['Type 1'].unique()
args=[]
for i in list(g):
args.append(df1[df1['Type 1']==i]['Total'])
w,p=stats.levene(*args)
if p<=0.05:
print('方差齐次性不成立')
else:
print('方差齐次，可进行下一步操作')
f,p1=stats.f_oneway(*args)
if p1<=0.05:
print('属性对种族值有影响')
else:
print('没啥影响')
输出结果为：
方差齐次，可进行下一步操作
属性对种族值有影响
看来属性对种族值是有影响的，那么究竟哪个属性高人一等呢？还需进一步分析。。
df2=df1.groupby(df1['Type 1'])
print(df2['Total'].mean().sort_values(ascending=False))
mc = MultiComparison(df1['Total'], df1['Type 1'])
result1 = mc.tukeyhsd()
print(result1)
得到以下结果：
Type 1
Dragon 550.531250
Steel 487.703704
Flying 485.000000
Psychic 475.947368
Fire 458.076923
Rock 453.750000
Dark 445.741935
Electric 443.409091
Ghost 439.562500
Ground 437.500000
Ice 433.458333
Water 430.455357
Grass 421.142857
Fighting 416.444444
Fairy 413.176471
Poison 403.814815
Normal 401.683673
Bug 378.927536
Tukeyhsd()函数的部分结果如下：
Multiple Comparison of Means - Tukey HSD,FWER=0.05
=====================================================
group1 group2 meandiff lower upper reject
-----------------------------------------------------
Bug Dark 66.8144 -20.6521 154.2809 False
Bug Dragon 171.6037 85.0853 258.1221 True
Bug Electric 64.4816 -13.5618 142.5249 False
Bug Fairy 34.2489 -75.2847 143.7826 False
Bug Fighting 37.5169 -54.3114 129.3452 False
Bug Fire 79.1494 4.8622 153.4366 True
Bug Flying 106.0725 -101.9711 314.116 False
Bug Ghost 60.635 -25.8834 147.1534 False
Bug Grass 42.2153 -26.4095 110.8401 False
Bug Ground 58.5725 -27.9459 145.0909 False
Bug Ice 54.5308 -41.3339 150.3955 False
Bug Normal 22.7561 -40.8162 86.3284 False
Bug Poison 24.8873 -66.941 116.7156 False
Bug Psychic 97.0198 24.6145 169.4252 True
Bug Rock 74.8225 -3.2208 152.8658 False
Bug Steel 108.7762 16.9479 200.6045 True
Bug Water 51.5278 -10.3811 113.4367 False
Dragon Electric -107.1222 -201.1059 -13.1384 True
Dragon Fairy -137.3548 -258.7624 -15.9472 True
Dragon Fighting -134.0868 -239.7969 -28.3767 True
Dragon Fire -92.4543 -183.3431 -1.5655 True
Dragon Flying -65.5312 -280.064 149.0015 False
Dragon Ghost -110.9688 -212.1005 -9.837 True
Dragon Grass -129.3884 -215.7108 -43.066 True
Dragon Ground -113.0312 -214.163 -11.8995 True
Dragon Ice -117.0729 -226.3076 -7.8382 True
Dragon Normal -148.8476 -231.2104 -66.4848 True
Dragon Poison -146.7164 -252.4265 -41.0063 True
Dragon Psychic -74.5839 -163.9412 14.7734 False
Dragon Rock -96.7812 -190.765 -2.7975 True
Dragon Steel -62.8275 -168.5376 42.8825 False
Dragon Water -120.0759 -201.1616 -38.9901 True
Reject列为False,并且lower<0<upper，则说明两个属性的种族值没有显著差别。反之则说明有显著差别。
那么结果表明，龙属性的pokemon与绝大多数其他属性的种族值都有显著差异，并且其平均值也是最高的。而虫属性则荣获种族值最后一名！并且和几个种族值高的属性都有显著差异，太惨了。。我盖诺赛克特表示不服！