用Python學分析 - 單因素方差分析

 

單因素方差分析(One-Way Analysis of Variance)

判斷控制變量是否對觀測變量產生了顯著影響

分析步驟

1. 建立檢驗假設

　　 - H0：不同因子水平間的均值無差異
　　- H1：不同因子水平間的均值有顯著差異
　　- 【注意】有差異，有可能是所有因子水平間都存在差異，也有可能只有兩個因子水平間的均值存在差異

2. 計算檢驗統計量F值

　　F = MSA / MSE
　　MSA = SSA / ( k - 1 )    MSA：組間均方, 對總體方差的一個估計
　　MSE = SSE / ( n - k )    MSE：組內均方,不論H0是否為真，MSE都是總體方差的一個無偏估計
　　SST = SSA + SSE        SST：總誤差平方和，反映全部觀測值的離散情況
       SSA:組間誤差平方和，也稱水平項誤差平方和，反映各因子水平（總體）的樣本均值之間的差異程度
       SSE: 組內誤差平方和

3. 確定P值

4. 方差分析表

5. 根據給定的顯著性水平，並作出決策

　　根據F值進行假設檢驗
　　根據選定的顯著性水平，F值大於臨界值時，將拒絕原假設
　　根據P值進行假設檢驗

 6. 進一步分析

     方差齊性檢驗

     多重比較檢驗

　　- 確定控制變量的不同水平對觀測變量的影響程度
　　- 哪個水平的作用明顯區別於其他水平
　　- 哪個水平的作用是不顯著
　　- 等等

【python分析：用ols模塊進行計算】

# 引入數據
import pandas as pd
data_value = { '無促銷':[23,19,17,26,28,23,24,30],
              '被動促銷':[26,22,20,30,36,28,30,32],
              '主動促銷':[30,23,25,32,48,40,41,46]}# 因變量
da = pd.DataFrame( data_value ).stack()
da.columns = ['水平','觀測值']

# ols模塊進行分析

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

formula = '{} ~ {}'.format(da.columns[1], da.columns[0])
model = ols( formula, da ).fit()
anovat = anova_lm(model)
print(anovat)

輸出結果：

【python分析：用自定義函數進行計算】

def ANOVA_oneway( df, a = 0.05 ):
    from scipy.stats import f
    '''
    進行單因素方差分析
    輸入值：df - pd.DataFrame，第一列為水平，第二列為觀測值；a - 顯著性水平，默認為0.05
    返回類型：字典
    返回值：方差分析相關數據
    '''
    res = { 'SSA':0, 'SST':0 }
    mu = df[df.columns[1]].mean()
    da = df.groupby( df.columns[0] ).agg( {df.columns[1]:['mean','count']})
    da.columns = ['mean','count']
    res['df_A'] = len(list(da.index)) - 1        # 自由度
    # 組間誤差平方和
    for row in da.index:
        res['SSA'] += (da.loc[row,'mean'] - mu )**2 * da.loc[row,'count']
    # 總誤差平方和
    for e in df[df.columns[1]].values:
        res['SST'] += (e - mu )**2         
    res['SSE'] = res['SST'] - res['SSA']         # 組內誤差平方和
    res['df_E'] = len(df) - res['df_A'] - 1      # 殘差自由度
    res['df_T'] = len(df) - 1                    # 總和自由度
    res['MSA'] = res['SSA'] / res['df_A']        # 組間均方
    res['MSE'] = res['SSE'] / res['df_E']        # 組內均方
    res['F'] = res['MSA'] / res['MSE']           # F值
    res['p_value'] = 1 - f(res['df_A'],res['df_E'] ).cdf( res['F'])  #p值
    res['a'] = a
    res['F_alpha'] = f(res['df_A'],res['df_E'] ).ppf( 1-a ) # 基於顯著性水平a的F臨界值 
    return res

def print_ANOVA_oneway( d, maxedg = 90 ):
    '''
    打印單因素方差分析表
    輸入值：d - dict字典，包含分析表所需要的數據; maxedg - 打印輸出時裝飾分隔符的最大長度
    '''
    title = '【單因素方差分析表】'
    print( title.center( maxedg ))
    print( '=' *  maxedg )
    print( '{:^12s}|{:^16s}|{:^6s}|{:^16s}|{:^12s}|{:^10s}|'.format('誤差來源','平方和','自由度','均方和','F','p值'))
    print( '-' *  maxedg )
    print( '{:8s}|{:>18,.4f} |{:>8d} |{:>18,.4f} |{:>11.6f} |{:>10.3%} |'.format( '組間（因子影響）',d['SSA'],d['df_A'],d['MSA'],d['F'],d['p_value']))
    print( '{:10s}|{:>18,.4f} |{:>8d} |{:>18,.4f} |'.format( '組內（誤差）',d['SSE'],d['df_E'],d['MSE']))
    print( '{:14s}|{:>18,.4f} |{:>8d} |'.format( '總和',d['SST'],d['df_T']))
    print( '-' *  maxedg )
    print('備注：顯著性水平為 {:.2%} 時，F的臨界值是 {:.6f}。'.format(d['a'],d['F_alpha']))


p = 0.95 # 設定置信度水平
maxedg = 93 # 設定輸出時裝飾分隔符的最大長度
# 計算並輸出單因素方差分析表
res = ANOVA_oneway( da, a = 1-p )
print_ANOVA_oneway( res, maxedg = maxedg )