[[331958]]

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

首先，還是先簡(jiǎn)介一下方差分析。

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

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

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

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

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

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

然后，我們?cè)僬f明一下數(shù)據(jù)集。

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

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]

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

H0：u1=u2=u3=u4=u5

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

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

圖1. 所用數(shù)據(jù)的基本情況

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

圖2. 方差分析表

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

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

方法1：scipy

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

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('拒絕原假設(shè)，u1、u2、u3、u4、u5不全相等')  
else:  
print('接受原假設(shè)，u1=u2=u3=u4=u5') 

結(jié)果如圖3所示。

圖3. 方法1的計(jì)算結(jié)果

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

方法2：statsmodels

方法2用的是python的另一個(gè)統(tǒng)計(jì)學(xué)庫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) 

結(jié)果如圖4所示。

圖4. 方法2的計(jì)算結(jié)果

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

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

總結(jié)

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