[[331958]]

在進行數據分析時，我們往往會遇到要對某個變量的影響因素進行分析的情況，而影響一事物的因素往往是很多的。比如在化工生產中，有溫度、壓力、劑量、反應時間等因素。每一因素的改變都有可能影響產品的數量和質量。我們往往要找出對產品質量有顯著影響的那些因素。而方差分析就是根據試驗的結果進行分析，鑒別各個有關因素對試驗結果影響的有效方法，本文主要講述如何用python中的兩種方法來進行方差分析。

首先，還是先簡介一下方差分析。

方差分析（Analysis of Variance，ANOVA）又稱“變異數分析”或“F檢驗”，是由羅納德·費舍爾（Ronald Aylmer Fisher）發明的，用于兩個及兩個以上樣本均數差別的顯著性檢驗，其原理是認為不同處理組的均數間的差別基本來源有兩個：

(1) 實驗條件，即不同的處理造成的差異，稱為組間差異。用變量在各組的均值與總均值之偏差平方和的總和表示，記作SSa，組間自由度dfa。

(2) 隨機誤差，如測量誤差造成的差異或個體間的差異，稱為組內差異，用變量在各組的均值與該組內變量值之偏差平方和的總和表示， 記作SSe，組內自由度dfe。

總偏差平方和 SSt = SSa + SSe。

組內SSe、組間SSa除以各自的自由度（組內dfe =n-m，組間dfa=m-1，其中n為樣本總數，m為組數），得到其均方MSe和MSa，一種情況是處理沒有作用，即各組樣本均來自同一總體，MSa/MSe≈1。另一種情況是處理確實有作用，組間均方是由于誤差與不同處理共同導致的結果，即各樣本來自不同總體。那么，MSa>>MSe（遠大于）。

MSa/MSe比值構成F分布。用F值與其臨界值比較，推斷各樣本是否來自相同的總體。

然后，我們再說明一下數據集。

數據集非常簡單，只有5組數值，每組數值有4個，共20個數字。分別命名為group1、group2、group3、group4和group5，數值都是隨意設置的，沒有什么要求，這里大家也可以根據自己的意愿設置數據。在這里，筆者專門將數據量設置得比較小，這樣方便觀察數據的之間的差異，我們的重點是方差分析的方法，而這里我們主要講的是單因素方差分析法。

group1 = [29.6, 24.3, 28.5, 32.0]

group2 = [27.3, 32.6, 30.8, 34.8]

group3 = [5.8, 6.2,11.0, 8.3]

group4 = [21.6, 17.4, 18.3, 19.0]

group5 = [29.2, 32.8, 25.0, 24.2]

設u1、u2、u3、u4和u5分別是這5個樣本所屬總體的均值，我們用單因素方差分析來檢驗下面的假設。

H0：u1=u2=u3=u4=u5

H1：u1、u2、u3、u4和u5不全相等

為了能更直觀了解這5組數據，我們首先手工計算一下這些數據的相關參數。這5組數據的總體情況如圖1所示。

圖1. 所用數據的基本情況

在圖1中，每列數據就是一個水平，這是一個統計學用語，水平和就是每組4個數值的總和，每組數據平均值分別是a1=28.6，a2=31.375，a3=7.825，a4=19.075，a5=27.8，全部20個數據的平均值為A=(a1+a2+a3+a4+a5)/5=114.675/5=22.935。所以總偏差平方和為ST=1616.65，此值為20個數據中每個數據與A的差的平方的總和，誤差平方和為SE=135.82，此值為每組數據中每個數據與這組數據的平均值的差的平方之和，效應平方和為SA=1480.83，此值為每組數據的平均值與A的差的平方之和，也等于ST減去SE的差。由此我們可以得出本例的方差分析表，如圖2所示。

圖2. 方差分析表

圖2中的因素就是各組數據間的差異，這個可以是隨機的，也可以是人為的，而誤差就是每組數據的之間差異。我們可以看到本例中得到的F值為40.8848，遠大于查表得到的F值F0.05(4,15)，其值為3.06，至于F0.05(4,15)的值我們同樣可以用python得出，后面會有講解。

以上就是這個例子的手工計算過程，下面我們用python來計算一下該例。

方法1：scipy

方法1用的庫是scipy，這是python中科學計算最常用的庫，其代碼如下，記得輸入前面的5組數據。 

from scipy import stats  
F, p = stats.f_oneway(group1, group2, group3, group4, group5)  
F_test = stats.f.ppf((1-0.05), 4, 15) 
print('F值是%.2f，p值是%.9f' % (F,p))  
print('F_test的值是%.2f' % (F_test))  
if F>=F_test:  
print('拒絕原假設，u1、u2、u3、u4、u5不全相等')  
else:  
print('接受原假設，u1=u2=u3=u4=u5') 

結果如圖3所示。

圖3. 方法1的計算結果

scipy的單因素方差分析比較簡單，只要調用stats模塊的f_oneway方法即可，在f_oneway中輸入各組數據，然后會自動返回兩個數值F與p，第一個數值F就表示我們算出的F值，和圖2中的F值一樣，而第二個值p就是這個F值所對應的概率，也就是假設檢驗問題中，由檢驗統計量的樣本觀察值得出的原假設可被拒絕的最小顯著性水平。在這里我們既可以通過F值來判斷，也可以通過p值來判斷，因為F大于F_test，落入了拒絕域，所以拒絕原假設，而p值也遠小于α分位數（這里為0.05），所以也拒絕原假設。而這里的F_test就是圖2中的F0.05(4,15)，計算方法就是用stats.f.ppf((1-0.05), 4, 15)，這里ppf的意思是Percent point function，也就是百分點函數，它是Cumulative distribution function（累積分布函數）的逆運算，這里需要注意的是ppf的第一個參數要輸入1-0.05，0.05也就是我們設定的顯著性水平α，其值通常取0.05，而第二個和第三個參數是兩個自由度，這兩個自由度分別是4和15，其求法在前面原理部分已經講過。

方法2：statsmodels

方法2用的是python的另一個統計學庫statsmodels，其代碼如下。 

import statsmodels.api as sm  
import pandas as pd  
from statsmodels.formula.api import ols  
num = sorted(['g1', 'g2', 'g3','g4', 'g5']*4)  
data = group1 + group2 + group3 + group4 + group5  
df = pd.DataFrame({'num':num, 'data': data})  
mod = ols('data ~ num', data=df).fit()          
ano_table = sm.stats.anova_lm(mod, typ=2)  
print(ano_table) 

結果如圖4所示。

圖4. 方法2的計算結果

從圖4中我們可以看到，得出的結果和前面手算以及scipy的結果一樣（部分小數精度問題可以忽略不計），圖中sum_sq列就表示平方和，df列就代表了自由度，這里還給出了p值就是PR(>F)列，信息比scipy要豐富一些。

從代碼上來看，statsmodels也同樣很簡單，只比scipy稍微復雜了一點，但卻提供了更多的信息。這里有幾點要注意的。一是我們生成了一個名為num的變量和一個名為data的變量，這兩個都是list類型，又用二者生成了名為df的pandas.DataFrame變量，這樣做的原因是statsmodels中普遍使用DataFrame數據格式，如果使用list類型會更麻煩一些。而data是把前面group1到group5中的數據放在了一個list中，num則是存放每個數據所對應的數據組信息，g1就代表這個數值屬于group1，g2則是對應group2，以此類推。這里還有一點要注意，就是num中數據格式最好是字符格式的，比如’a1’、‘num3’這樣的，不要是數字格式的，比如1、3、6.9這樣的，因為數字格式的數據很有可能會參與計算，最終的結果可能會出錯。第二點是mod = ols('data ~ num', data=df).fit()中的公式data ~ num，很多人對這一點很困惑，這種公式的使用方法來自于python的另一個庫patsy，其主要用于描述統計模型（尤其是線性模型），符號~前面的部分代表了y軸數據，后面的部分代表了x軸數據，根據這二者生成一個線性模型，ols中第二個參數data則是要輸入的數據源，一般是DataFrame格式，前面公式中符號~前后的名稱都要是data中的列名，這種方法確實有些奇怪，部分原因是patsy借鑒了R語言的一些用法。第三點是ano_table = sm.stats.anova_lm(mod, typ=2)中，typ=2的意思是DataFrame，typ共有3個值，分別是1、2和3，其中2代表了DataFrame格式。

總結

對比scipy和statsmodels這兩種方法，可以說是各有優勢。scipy是一個通用型庫，其包含了科學計算的多種模塊，統計分析只是其中一部分，而statsmodels是一個專門進行統計分析的庫，二者在功能上有一些差別，statsmodels在統計分析上更專業一些。而scipy的語法更符合python常用的語法，statsmodels的語法有些接近于R語言，對初學者可能有些陌生。所以大家可以根據自己的需要來選擇合適的方法。