方差分析(Anova)

1.單因素方差分析：

單因素方差分析：只有一個因素A對實驗指標有影響，假設因素A有r個水平，分別在第i個水平下進行多次獨立的觀察，所得到的實驗指標數據如下：

A_1：N(μ₁,σ²)   X₁₁  X₁₂   ... X_1n1

A_2：N(μ₂,σ²)   X₂₁  X₂₂   ... X_2n2

Ar：N(μ_r,σ²)    X_r1   X_r2   ... X_rnr

注意：每個水平的觀測次數不一定一樣

各總體間相互獨立，因此有下面的模型：

 X_ij就是第i個水平的第j個觀測值，μ_{i就是第i個水平的理論均值，εi}顯示隨機誤差（誤差服從正態分布）

分析因素A對於實驗指標是否有顯著影響，可以看因素A不同水平的均值是否有顯著差異，因此有如下假設：

原假設：H₀:μ₁=μ₂=...μ_r

備選假設 H1：既是均值不全相等

X_ij有偏差，要不就是由於不同水平的均值不同，又或者是隨機誤差的存在，因此全部X_ij之間的差異的公式如下：

上面這個叫總偏差平方和 

有A因素引起的 差異叫效應平方和S_{A （反應的是在因素A的不同水平下，樣本均值和總體數據均值差異的平方和），}隨機誤差引起的差異，叫做誤差平方和S_{E （反應是在因素A的各個取值下，每組觀察數據與這組數據均值的平方誤差之和，反應的是隨機誤差的大小）}

首先計算誤差平方和 ，這樣個體之間的差異的每個水平的均值沒有關系，因此有如下：

 綜合上述表達，得到：

 總偏差平方和減去誤差平方和，得到

 S_E如果除以σ²則會符合自由度為n_i-1的卡方分布

當H0為真的時候，但是我們不知道σ^2，因此為了抵消這個未知量，我們構造的檢驗統計量為：

 

 我們最終只會關系p值，如果p>0.05則接受原假設，否則拒絕原假設

例子：

import pandas as pd
import numpy as np

from scipy import stats
from statsmodels.formula.api import ols
from 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)      #  stats.f_oneway函數就可以直接算出檢驗假設的f值和p值

為什么要做方差記性檢驗：方差齊性檢驗是方差分析的重要前提，是方差可加性原則應用的一個條件

方差的齊性檢驗,如果p<0.05則拒絕原假設，即是方差不齊性 

 如果手動去計算：

#首先將數據改成DataFrame形式
values = A1.copy()
groups = []
for i in range(1, len(data)):
    values.extend(data[i])  #extend() 函數用於在列表末尾一次性追加另一個序列中的多個值

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})


#單因素分析
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
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

 2.雙因素方差分析：

雙因素方差分析和多因素方差分析在原理上是一致的，

雙因素方差分析就是在因素A,B作用下試驗的指標，因素A有r個水平，因素B有s個水平，在A,B的不同水平下得到的試驗結果如下：

 並設有條件

X_ijk獨立，數學模型如下：

每一個格子都有一個平均值，每一行每一列也有平均值，這里先定義均值：

μ是總的均值，再定義兩個公式：

αi為水平Ai上的效應，βj為水平Bj的效應 ，很顯然

 將其代入到前面的公式里面，得到;

這個模型就會得到三個假設檢驗問題

因素A對於實驗結果是否帶來了顯著效果

因素B對於實驗結果是否帶來了顯著效果

兩者組合是否帶來了顯著效果

因素A的i水平和因素B的j水平的平均值;

 因素A的i水平上的平均值：

 因素B的j水平均值：

 總的均值：

 總偏差平方和：

 

 其中S_E是誤差平方和，S_A和S_B分別是因素A和B的效應平方和，S_AxB是A和B的組合效應平方和 

ST的自由度是rst-1，SE的自由度是rs(t-1)，S_A的自由度是r-1，S_B自由度是s-1

當H₀₁為真時：

 這時候取顯著水平α，得到的拒絕域為：

 同理H₀₂拒絕域為：

H₀₃的拒絕域為：

 

 導入雙因素分析使用到的包：

import pandas as pd
import numpy as np

from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

# 這三個交互效果的可視化畫圖
from statsmodels.graphics.api import interaction_plot
import matplotlib.pyplot as plt
from pylab import mpl      # 顯示中文

# 這個看某個因素各個水平之間的差異
from statsmodels.stats.multicomp import pairwise_tukeyhsd

2.1無交互作用的情況：即是對每一個組合因素只進行一次獨立實驗，每一格只有一個值，稱為無重復實驗

下面進行雙因素方差分析，簡要流程是，先用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=['廣告','價格','銷量'])

進行方差分析

# 方差分析
price_lm = ols('銷量~C(廣告)+C(價格)', data=df_t2).fit()
table = sm.stats.anova_lm(price_lm, typ=2)

 

 即是不同價格和廣告都會對銷量有顯著差異

fig = interaction_plot(df_t2['廣告'],df_t2['價格'], df_t2['銷量'],
                        ylabel='銷量', xlabel='廣告')

 

 

# 廣告與銷量的影響 
print(pairwise_tukeyhsd(df_t2['銷量'], df_t2['廣告'], alpha=0.05)) # 第一個必須是銷量， 也就是我們的指標

 

2.2有交互作用的情況：

即是每個格子有不止一個值，也稱為重復試驗

#先構造數據

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=['燃料','推進器','射程'])

#方差分析

moore_lm = ols('射程~燃料+推進器+燃料:推進器', data=df_t3).fit()
table = sm.stats.anova_lm(moore_lm, typ=1)

 

 

fig = interaction_plot(df_t3['燃料'],df_t3['推進器'], df_t3['射程'],
                        ylabel='射程', xlabel='燃料')

 

 從這個圖里面可以看出， (A4, B1)和(A3, B2)組合的進程最好

print(pairwise_tukeyhsd(df_t3['射程'], df_t3['燃料']))

 

 都是False， 說明A因素各個水平之間無顯著差異

print(pairwise_tukeyhsd(df_t3['射程'], df_t3['推進器']))

 

 B也沒有差異