一個複雜的事物,其中每每有許多因素互相制約又互相依存。方差分析是一種經常使用的數據分析方法,其目的是經過數據分析找出對該事物有顯著影響的因素、各因素之間的交互做用及顯著影響因素的最佳水平等。python
本文介紹了方差分析的基礎概念,詳細講解了單因素方差分析、雙因素方差分析的原理,而且給出了它們的python實踐代碼。
api
本文大綱:數據結構
關於方差分析的基礎概念app
單因素方差分析原理及python實現ide
雙因素方差分析原理及python實現函數
在科學試驗和生產實踐中,影響某一事物的因素每每不少,好比化工生產中,像原料成分,劑量,反應溫度,壓力等等不少因素都會影響產品的質量,有些因素影響較大,有些影響較小, 爲了使生產過程穩定, 保證優質高產, 就有必要找出對產品質量有顯著影響的因素。如何找到影響因素呢?就須要試驗, 方差分析就是根據試驗的結果進行分析, 鑑別各個有關因素對試驗結果影響程度的有效方法。而根據涉及到的因素個數的不一樣, 又能夠把方差分析分爲單因素方差分析、多因素方差分析等。工具
下面咱們先重點研究單因素方差分析, 經過一個例子,引出方差分析中的幾個概念:學習
某保險公司想了解一下某險種在不一樣的地區是否有不一樣的索賠額。因而他們就蒐集了四個不一樣地區一年的索賠額狀況的記錄以下表:spa
嘗試判斷一下, 地區這個因素是否對與索賠額產生了顯著的影響?code
這個問題就是單因素方差分析的問題, 具體解決方法後面會說, 首先先由這個例子弄清楚幾個概念:
試驗指標:方差分析中, 把考察的試驗結果稱爲試驗指標, 上面例子裏面的「索賠額」。
因素:對試驗指標產生影響的緣由稱爲因素, 如上面的「地區」
這個類比的話, 就相似於y就是試驗指標, 某個類別特徵x就是因素, 類別特徵x的不一樣取值就是水平。那麼經過方差分析, 就能夠獲得某個類別特徵對於y的一個影響程度了吧, 這會幫助分析某個類別特徵的重要性!
所謂單因素方差分析, 就是僅考慮有一個因素對試驗指標的影響。假如因素有個水平, 分別在第個水平下進行屢次獨立的觀測, 所獲得的試驗指標數據以下:
注意這裏的不必定同樣,上面的例子。各整體間相互獨立,所以咱們會有下面的模型:簡單解釋一下上面這個在說啥:就是第個水平的第個觀測值,上面例子裏面就是第個地區第次的索賠額。表示第個水平的理論均值, 後面的表示的隨機偏差, 假設這個服從正態。第一個等式的意思就是某個觀測值能夠用某水平下的均值加一個偏差來表示。若是咱們想判斷某個因素對於試驗指標是否有顯著影響, 很直觀的就是咱們看看因素不一樣的水平下試驗指標的理論均值是否有顯著差別, 即理論均值是否徹底相同, 若是有顯著差別, 就說明不一樣的水平對試驗指標影響很大, 即對試驗指標有顯著影響。這也是方差分析的目標, 故把問題轉換成了比較不一樣水平下試驗指標的均值差別。顯著在這裏的意思是差別達到的某種程度。
基於上面的分析, 咱們就能夠把方差分析也當作一個檢驗假設的問題, 並有了原假設和備擇假設:
:
那麼這個假設檢驗的問題怎麼驗證呢?咱們得先分析一下, 爲啥各個會有差別?從上面的模型中, 咱們能夠看到, 因此第一個可能就是可能有差別, 好比, 那麼很容易就大於。另外一個可能就是隨機偏差的存在。在這樣的啓發下,咱們得找一個衡量所有之間差別的量, 就是下面這個了:
這個叫作總誤差平方和,若是這個越大, 就表示之間的差別就越大。這裏的表示總的觀測值個數:接下來,咱們把這個平方和分解開爲兩部分:一部分是因爲因素引發的差別, 這個叫作效應平方和, 另外一部分是因爲隨機偏差引發的差別,這個叫作偏差平方和
關於, 先固定一個, 此時對應的全部觀測值, 他們之間的差別與每一個水平的理論平均值就沒有關係了, 而是取決於隨機偏差, 反應這些觀察值差別程度的量其中
綜合全部的水平, 就能夠獲得偏差平方和的公式以下:
而上面二者相減, 就會獲得效應平方和:
經過上面的分析,咱們會獲得如下三點結論:
(這個分解式爲上面模型的方差分析)
這裏是由於:
後面這部分的加和若是除以的話會服從自由度爲的卡方(具體看第二篇卡方分佈的定義),那麼前面又一個r水平的累加,根據卡方分佈的可加性可得這個東西服從的卡方
基於上面的分析,會獲得一個單因素試驗方差分析表:
這個表就把上面全部的分析都給總結好了。但實際使用中,咱們確定是不會手算的,而且通常也不看F的值,咱們是看p值的。
下面就用python實現一下上面的那個索賠額的例子, 看看單因素方差分析是怎麼作的:
import pandas as pdimport numpy as np上面這段程序應該很容易懂, 首先前面是把數據構造出來, 而後進行一個方差的齊性檢驗, 這個用stats.levene函數, 這個的做用是要保證方差在每一個水平上某種程度上(顯著水平)是一致的, 這時候才能進行後面的均值分析, 由於方差分析的實質是檢驗多個水平的均值是否有顯著差別,若是各個水平的觀察值方差差別太大,只檢驗均值之間的差別就沒有意義了,因此要進行方差齊性檢驗。
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.f_oneway函數就能夠直接算出檢驗假設的值和值。咱們這裏關注的是值, 拿值和給出的(通常是0.05)比, 若是,咱們就接受原假設,不然拒絕原假設,這個例子中是0.134,大於,故接受原假設,認爲不一樣的地區的索賠額沒有顯著差別。
因此單因素方差這塊通常是懂了原理以後,用軟件去分析,能看懂就算入門了。固然這個若是手算的話,思路就是須要先求,而後根據上面的公式計算,計算完了以後除以自由度而後相除獲得值,而後比較和的大小,當,拒絕原假設,不然接受原假設。必定要注意這個值和值的比較標準是不一樣的。由於這是兩種假設檢驗的方法,值比較的這種是基於值法,而的那種是臨界值法。 上面的例子咱們還能夠進行那種單因素方差表的顯示格式:首先改一下數據的格式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
在不少狀況下, 只考慮一個指標對觀察值的影響顯然是不夠的, 這時就會用到多因素方差分析。雙因素方差分析和多因素方差分析原理上一致, 下面給出一種兩個因素之間有交互的一種形式寫法做爲補充。所謂雙因素方差分析, 就是有兩個因素做用於試驗的指標, 因素有個水平, 因素有個水平. 現對因素的水平的每對組合都做次試驗,也會獲得一個表:
並設 這裏的獨立, 類比着單因素方差分析那裏, 咱們就會先有下面的數學模型: 各獨立 這裏的表示的是第個因素第個因素下的第個觀測值。是組合下的全部觀測值的平均數(平均效應)。是隨機偏差, 這個其實和單因素那裏的理解是一個意思, 上面的單因素的那個表格放在雙因素這裏就至關於這裏的其中一個小格子了。那麼就開始引入一些新的公式, 由於既然每一個格子裏面有平均, 那麼每一行的格子和每一列的格子也會有平均, 總體上也會有平均, 因此下面就定義三個公式:
這個就是雙因素試驗方差分析的數學模型。對於這個模型, 咱們就會有三個假設檢驗的問題了:
因素A對於試驗結果是否帶來了顯著影響
因素B對於試驗結果是否帶來了顯著影響
二者的組合對於試驗結果是否帶來了顯著影響
與單因素的狀況相似, 咱們依然是採用平方和分解的方式進行驗證。首先咱們得先計算四個平均值:
因素A的水平因素B的水平的平均值:
因素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.一、無交互做用的狀況
因爲不考慮交互做用的影響,對每個因素組合 只需進行一次獨立試驗,稱爲無重複試驗。數據:考慮三種不一樣形式的廣告和五種不一樣的價格對某種商品銷量的影響。選取某市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.二、有交互做用的狀況
因爲因素有交互做用,須要對每個因素組合 分別進行 次 重複試驗,稱這種試驗爲等重複試驗。
數據:機率論課本上的那個例子, 火箭的射程與燃料的種類和推動器的型號有關,現對四種不一樣的燃料與三種不一樣型號的推動器進行試驗,每種組合各發射火箭兩次,測得火箭的射程結果以下(設顯著性水平爲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因素各個水平之間無顯著差別。
兩個實驗到這裏就結束了, 這裏再補充兩點別的知識:
1. ols函數裏面公式的寫法
'射程~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, 而且此時燃料和推動器的數據都發生了變化。
方差分析這塊到這裏就結束了, 隨着這篇文章的結束也意味着機率統計的知識串聯也到了尾聲, 簡單的回顧一下本篇的內容, 這篇文章主要是在實踐的角度進行的分析, 方差分析在統計中仍是很經常使用的, 比較適合類別因素對於數值指標的影響程度:
首先從單因素方差分析入手, 這個只考慮了一個因素對於指標的影響, 先分析了原理,而後基於python進行了實現。實際應用中,通常是會點原理,而後使用工具實現方差分析,會看結果,這樣就算入門了。
實際應用中, 或許能夠經過這種方法去分析類別特徵的重要性或者關聯性,以及類別和類別特徵之間的交互做用等。