常用數據分析方法:方差分析及實現!

↑↑↑關注後"星標"Datawhale

作者：吳忠強，Datawhale優秀學習者，東北大學

一個複雜的事物,其中往往有許多因素互相制約又互相依存。方差分析是一種常用的數據分析方法，其目的是通過數據分析找出對該事物有顯著影響的因素、各因素之間的交互作用及顯著影響因素的最佳水平等。

本文介紹了方差分析的基礎概念，詳細講解了單因素方差分析、雙因素方差分析的原理，並且給出了它們的python實踐代碼。

本文大綱：

關於方差分析的基礎概念

單因素方差分析原理及python實現

雙因素方差分析原理及python實現

一、關於方差分析的基礎概念在科學試驗和生產實踐中，影響某一事物的因素往往很多，比如化工生產中，像原料成分，劑量，反應溫度，壓力等等很多因素都會影響產品的質量，有些因素影響較大，有些影響較小， 為了使生產過程穩定， 保證優質高產， 就有必要找出對產品質量有顯著影響的因素。如何找到影響因素呢？就需要試驗， 方差分析就是根據試驗的結果進行分析， 鑑別各個有關因素對試驗結果影響程度的有效方法。而根據涉及到的因素個數的不同， 又可以把方差分析分為單因素方差分析、多因素方差分析等。

下面我們先重點研究單因素方差分析， 通過一個例子，引出方差分析中的幾個概念：

某保險公司想了解一下某險種在不同的地區是否有不同的索賠額。於是他們就搜集了四個不同地區一年的索賠額情況的記錄如下表：

嘗試判斷一下， 地區這個因素是否對與索賠額產生了顯著的影響？

這個問題就是單因素方差分析的問題， 具體解決方法後面會說， 首先先由這個例子弄清楚幾個概念：

試驗指標：方差分析中， 把考察的試驗結果稱為試驗指標， 上面例子裡面的「索賠額」。

因素：對試驗指標產生影響的原因稱為因素， 如上面的「地區」

水平：因素中各個不同狀態， 比如上面我們有A1， A2， A3， A4四個狀態， 四個水平。

這個類比的話， 就類似於y就是試驗指標， 某個類別特徵x就是因素， 類別特徵x的不同取值就是水平。那麼通過方差分析， 就可以得到某個類別特徵對於y的一個影響程度了吧， 這會幫助分析某個類別特徵的重要性！

二、單因素方差分析原理及python實現所謂單因素方差分析， 就是僅考慮有一個因素
注意這裡的

基於上面的分析， 我們就可以把方差分析也看成一個檢驗假設的問題， 並有了原假設和備擇假設：

那麼這個假設檢驗的問題怎麼驗證呢？我們得先分析一下， 為啥各個

這個叫做總偏差平方和，如果這個越大， 就表示
關於

綜合所有的水平， 就可以得到誤差平方和的公式如下：

由於

通過上面的分析，我們會得到以下三點結論：

後面這部分的加和如果除以
我們上面說
這時候構造出了F統計量。在原假設成立的時候，

基於上面的分析，會得到一個單因素試驗方差分析表：

這個表就把上面所有的分析都給總結好了。但實際使用中，我們肯定是不會手算的，並且一般也不看F的值，我們是看p值的。

下面就用python實現一下上面的那個索賠額的例子， 看看單因素方差分析是怎麼做的：

import pandas as pdimport numpy as np
from scipy import statsfrom statsmodels.formula.api import olsfrom statsmodels.stats.anova import anova_lm
# 這是那四個水平的索賠額的觀測值A1 = [1.6, 1.61, 1.65, 1.68, 1.7, 1.7, 1.78]A2 = [1.5, 1.64, 1.4, 1.7, 1.75]A3 = [1.6, 1.55, 1.6, 1.62, 1.64, 1.60, 1.74, 1.8]A4 = [1.51, 1.52, 1.53, 1.57, 1.64, 1.6]
data = [A1, A2, A3, A4]# 方差的齊性檢驗w, p = stats.levene(*data)if p < 0.05:    print('方差齊性假設不成立')
# 成立之後， 就可以進行單因素方差分析f, p = stats.f_oneway(*data)print(f, p)      #  2.06507381767795 0.13406910483160134
上面這段程序應該很容易懂， 首先前面是把數據構造出來， 然後進行一個方差的齊性檢驗， 這個用stats.levene函數， 這個的作用是要保證方差在每個水平上某種程度上(顯著水平)是一致的， 這時候才能進行後面的均值分析, 因為方差分析的實質是檢驗多個水平的均值是否有顯著差異，如果各個水平的觀察值方差差異太大，只檢驗均值之間的差異就沒有意義了，所以要進行方差齊性檢驗。
後面通過stats.f_oneway函數就可以直接算出檢驗假設的
所以單因素方差這塊一般是懂了原理之後，用軟體去分析，能看懂就算入門了。當然這個如果手算的話，思路就是需要先求
values = A1.copy()groups = []for i in range(1, len(data)):    values.extend(data[i])
for i, j in zip(range(4), data):    groups.extend(np.repeat('A'+str(i+1), len(j)).tolist())
df = pd.DataFrame({'values': values, 'groups': groups})df
數據長這個樣子了，也是我們一般見到的pandas的形式：
通過下面的方式做單因素方差分析：
anova_res = anova_lm(ols('values~C(groups)', df).fit())anova_res.columns = ['自由度', '平方和', '均方', 'F值', 'P值']anova_res.index = ['因素A', '誤差']anova_res        # 這種情況下看p值  >0.05 所以接受H0
結果如下：
這樣就會得到單因素方差分析表的格式。當然， 為了考慮的全面些， 我們應該評估檢驗的假設條件， 就是看看每個數據是不是真的服從正態。這裡就使用上一篇文章中學習到的判斷數據是不是服從正態的方法了Shapiro-Wilk test（小樣本情況下， 常用的正態檢驗方法）：
A1 = [1.6, 1.61, 1.65, 1.68, 1.7, 1.7, 1.78]A2 = [1.5, 1.64, 1.4, 1.7, 1.75]A3 = [1.6, 1.55, 1.6, 1.62, 1.64, 1.60, 1.74, 1.8]A4 = [1.51, 1.52, 1.53, 1.57, 1.64, 1.6]
data = [A1, A2, A3, A4]
from scipy.stats import shapiro
def normal_judge(data):    stat, p = shapiro(data)    if p > 0.05:        return 'stat={:.3f}, p = {:.3f}, probably gaussian'.format(stat,p)    else:        return 'stat={:.3f}, p = {:.3f}, probably not gaussian'.format(stat,p)
for d in data:    print(normal_judge(d))
結果如下：
stat=0.942, p = 0.660, probably gaussianstat=0.938, p = 0.655, probably gaussianstat=0.850, p = 0.096, probably gaussianstat=0.918, p = 0.489, probably gaussian
三、雙因素方差分析及python實現在很多情況下， 只考慮一個指標對觀察值的影響顯然是不夠的， 這時就會用到多因素方差分析。雙因素方差分析和多因素方差分析原理上一致， 下面給出一種兩個因素之間有交互的一種形式寫法作為補充。所謂雙因素方差分析， 就是有兩個因素
這裡的
那麼就開始引入一些新的公式， 因為既然每個格子裡面有平均， 那麼每一行的格子和每一列的格子也會有平均， 整體上也會有平均， 所以下面就定義三個公式：
我們稱
這個就是雙因素試驗方差分析的數學模型。對於這個模型， 我們就會有三個假設檢驗的問題了：
因素A對於試驗結果是否帶來了顯著影響
因素B對於試驗結果是否帶來了顯著影響
兩者的組合對於試驗結果是否帶來了顯著影響
與單因素的情況類似， 我們依然是採用平方和分解的方式進行驗證。首先我們得先計算四個平均值：
因素A的
因素A的
因素B的
總平均值：
有了上面的平均值， 我們就可以得到偏差平方和了， 總偏差平方和如下：
就得到了
其中
這裡也給出每個平方和的自由度，
和單因素方差分析那裡的思路是一樣的， 碰到具體問題的時候， 我們一般不會採用手算的形式， 如果手算的話， 思路和上面一樣， 就是先根據公式求四個平均值， 然後根據平均值求那四個平方和的東西， 求完了之後算三個
導入這次用到的包（依然是單因素分析時的ols和anova_lm）
import pandas as pdimport numpy as np
from scipy import statsfrom statsmodels.formula.api import olsfrom statsmodels.stats.anova import anova_lm
# 這三個交互效果的可視化畫圖from statsmodels.graphics.api import interaction_plotimport matplotlib.pyplot as pltfrom pylab import mpl      # 顯示中文
# 這個看某個因素各個水平之間的差異from statsmodels.stats.multicomp import pairwise_tukeyhsd
3.1、無交互作用的情況
由於不考慮交互作用的影響，對每一個因素組合 數據：考慮三種不同形式的廣告和五種不同的價格對某種商品銷量的影響。選取某市15家大超市，每家超市選用其中的一個組合，統計出一個月的銷量如下（設顯著性水平為0.05）：
下面進行雙因素方差分析，簡要流程是，先用pandas庫的DataFrame數據結構來構造輸入數據格式。然後用statsmodels庫中的ols函數得到最小二乘線性回歸模型。最後用statsmodels庫中的anova_lm函數進行方差分析。
dic_t2=[{'廣告':'A1','價格':'B1','銷量':276},{'廣告':'A1','價格':'B2','銷量':352},       {'廣告':'A1','價格':'B3','銷量':178},{'廣告':'A1','價格':'B4','銷量':295},       {'廣告':'A1','價格':'B5','銷量':273},{'廣告':'A2','價格':'B1','銷量':114},       {'廣告':'A2','價格':'B2','銷量':176},{'廣告':'A2','價格':'B3','銷量':102},       {'廣告':'A2','價格':'B4','銷量':155},{'廣告':'A2','價格':'B5','銷量':128},       {'廣告':'A3','價格':'B1','銷量':364},{'廣告':'A3','價格':'B2','銷量':547},       {'廣告':'A3','價格':'B3','銷量':288},{'廣告':'A3','價格':'B4','銷量':392},       {'廣告':'A3','價格':'B5','銷量':378}]df_t2=pd.DataFrame(dic_t2,columns=['廣告','價格','銷量'])df_t2
數據長這樣：
# 方差分析price_lm = ols('銷量~C(廣告)+C(價格)', data=df_t2).fit()table = sm.stats.anova_lm(price_lm, typ=2)table
結果如下：
可以發現這裡的p值都是小於0.05的， 所以我們要拒絕掉原假設， 即可認為不同的廣告形式， 不同的價格均造成商品銷量的顯著差異。
fig = interaction_plot(df_t2['廣告'],df_t2['價格'], df_t2['銷量'],                        ylabel='銷量', xlabel='廣告')
結果如下：
再來分析一下單因素各個水平之間的顯著差異：
# 廣告與銷量的影響  注意這個的顯著水平是0.01print(pairwise_tukeyhsd(df_t2['銷量'], df_t2['廣告'], alpha=0.01)) # 第一個必須是銷量， 也就是我們的指標
結果如下：	
這個可以得到的結論是在顯著水平0.01的時候， A2和A3的p值小於0.01， reject=True, 即認為A2和A3有顯著性差異。
3.2、有交互作用的情況
由於因素有交互作用，需要對每一個因素組合 
數據：概率論課本上的那個例子， 火箭的射程與燃料的種類和推進器的型號有關，現對四種不同的燃料與三種不同型號的推進器進行試驗，每種組合各發射火箭兩次，測得火箭的射程結果如下（設顯著性水平為0.01）：
dic_t3=[{'燃料':'A1','推進器':'B1','射程':58.2},{'燃料':'A1','推進器':'B1','射程':52.6},       {'燃料':'A1','推進器':'B2','射程':56.2},{'燃料':'A1','推進器':'B2','射程':41.2},       {'燃料':'A1','推進器':'B3','射程':65.3},{'燃料':'A1','推進器':'B3','射程':60.8},       {'燃料':'A2','推進器':'B1','射程':49.1},{'燃料':'A2','推進器':'B1','射程':42.8},       {'燃料':'A2','推進器':'B2','射程':54.1},{'燃料':'A2','推進器':'B2','射程':50.5},       {'燃料':'A2','推進器':'B3','射程':51.6},{'燃料':'A2','推進器':'B3','射程':48.4},       {'燃料':'A3','推進器':'B1','射程':60.1},{'燃料':'A3','推進器':'B1','射程':58.3},       {'燃料':'A3','推進器':'B2','射程':70.9},{'燃料':'A3','推進器':'B2','射程':73.2},       {'燃料':'A3','推進器':'B3','射程':39.2},{'燃料':'A3','推進器':'B3','射程':40.7},       {'燃料':'A4','推進器':'B1','射程':75.8},{'燃料':'A4','推進器':'B1','射程':71.5},       {'燃料':'A4','推進器':'B2','射程':58.2},{'燃料':'A4','推進器':'B2','射程':51.0},       {'燃料':'A4','推進器':'B3','射程':48.7},{'燃料':'A4','推進器':'B3','射程':41.4},]df_t3=pd.DataFrame(dic_t3,columns=['燃料','推進器','射程'])df_t3.head()
結果這樣：
下面是方差分析：
moore_lm = ols('射程~燃料+推進器+燃料:推進器', data=df_t3).fit()table = sm.stats.anova_lm(moore_lm, typ=1)table
這裡得到的結論就是燃料的P值是大於0.01的， 而推進器和兩者組合的p值都小於0.01， 並且兩者的組合非常小， 這就說明燃料對於火箭的射程沒有顯著影響， 而後兩者都有顯著影響，兩者的交互作用更是高度顯著。
fig = interaction_plot(df_t3['燃料'],df_t3['推進器'], df_t3['射程'],                        ylabel='射程', xlabel='燃料')
結果如下：
從這個圖裡面可以看出， (A4, B1)和(A3, B2)組合的進程最好。黃金搭檔。單因素差異性分析：
print(pairwise_tukeyhsd(df_t3['射程'], df_t3['燃料']))
結果：
	都是False， 說明A因素各個水平之間無顯著差異。
兩個實驗到這裡就結束了， 這裡再補充兩點別的知識：
'射程~C(燃料)+C(推進器)+C(燃料):C(推進器)' ：相當於射程是y(指標), 燃料和推進器是x(影響因素)， 三項加和的前兩項表示兩個主效應， 第三項表示考慮兩者的交互效應， 不加C也可。
'射程~C(燃料, Sum)*C(推進器, Sum)'和上面效果是一致的， 星號在這裡表示既考慮主效應也考慮交互效應*'銷量~C(廣告)+C(價格)'：這個表示不考慮交互相應
但是要注意， 考慮交互相應和不考慮交互相應導致的Se（殘差項）會不同， 所以會影響最終的結果。
stats.anova_lm(moore_lm, typ=1)這裡面的typ參數， 這個參數我嘗試還沒有完全搞明白到底是什麼意思， 這個參數有1,2，3 三個可選項， 分別代表著不同的偏差平方和的計算方法， 我在第二個實驗中嘗試過改這個參數，改成1的時候發現就加了一列mean_sq， 然後其他的沒變。
改成3的時候發現加一行Intercept， 並且此時燃料和推進器的數據都發生了變化。
四、寫到最後
方差分析這塊到這裡就結束了， 隨著這篇文章的結束也意味著概率統計的知識串聯也到了尾聲， 簡單的回顧一下本篇的內容， 這篇文章主要是在實踐的角度進行的分析， 方差分析在統計中還是很常用的， 比較適合類別因素對於數值指標的影響程度：
實際應用中， 或許可以通過這種方法去分析類別特徵的重要性或者關聯性，以及類別和類別特徵之間的交互作用等。
本文電子版 後臺回復 概率統計 獲取
「感謝你的分享，點讚，在看三連↓