用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模块进行计算】
1 # 引入数据 2 import pandas as pd 3 data_value = { '无促销':[23,19,17,26,28,23,24,30], 4 '被动促销':[26,22,20,30,36,28,30,32], 5 '主动促销':[30,23,25,32,48,40,41,46]}# 因变量 6 da = pd.DataFrame( data_value ).stack() 7 da.columns = ['水平','观测值'] 8 9 # ols模块进行分析 10 11 from statsmodels.formula.api import ols 12 from statsmodels.stats.anova import anova_lm 13 14 formula = '{} ~ {}'.format(da.columns[1], da.columns[0]) 15 model = ols( formula, da ).fit() 16 anovat = anova_lm(model) 17 print(anovat)
输出结果：
【python分析：用自定义函数进行计算】
1 def ANOVA_oneway( df, a = 0.05 ): 2 from scipy.stats import f 3 ''' 4 进行单因素方差分析 5 输入值：df - pd.DataFrame，第一列为水平，第二列为观测值；a - 显著性水平，默认为0.05 6 返回类型：字典 7 返回值：方差分析相关数据 8 ''' 9 res = { 'SSA':0, 'SST':0 } 10 mu = df[df.columns[1]].mean() 11 da = df.groupby( df.columns[0] ).agg( {df.columns[1]:['mean','count']}) 12 da.columns = ['mean','count'] 13 res['df_A'] = len(list(da.index)) - 1 # 自由度 14 # 组间误差平方和 15 for row in da.index: 16 res['SSA'] += (da.loc[row,'mean'] - mu )**2 * da.loc[row,'count'] 17 # 总误差平方和 18 for e in df[df.columns[1]].values: 19 res['SST'] += (e - mu )**2 20 res['SSE'] = res['SST'] - res['SSA'] # 组内误差平方和 21 res['df_E'] = len(df) - res['df_A'] - 1 # 残差自由度 22 res['df_T'] = len(df) - 1 # 总和自由度 23 res['MSA'] = res['SSA'] / res['df_A'] # 组间均方 24 res['MSE'] = res['SSE'] / res['df_E'] # 组内均方 25 res['F'] = res['MSA'] / res['MSE'] # F值 26 res['p_value'] = 1 - f(res['df_A'],res['df_E'] ).cdf( res['F']) #p值 27 res['a'] = a 28 res['F_alpha'] = f(res['df_A'],res['df_E'] ).ppf( 1-a ) # 基于显著性水平a的F临界值 29 return res 30 31 def print_ANOVA_oneway( d, maxedg = 90 ): 32 ''' 33 打印单因素方差分析表 34 输入值：d - dict字典，包含分析表所需要的数据; maxedg - 打印输出时装饰分隔符的最大长度 35 ''' 36 title = '【单因素方差分析表】' 37 print( title.center( maxedg )) 38 print( '=' * maxedg ) 39 print( '{:^12s}|{:^16s}|{:^6s}|{:^16s}|{:^12s}|{:^10s}|'.format('误差来源','平方和','自由度','均方和','F','p值')) 40 print( '-' * maxedg ) 41 print( '{:8s}|{:>18,.4f} |{:>8d} |{:>18,.4f} |{:>11.6f} |{:>10.3%} |'.format( '组间（因子影响）',d['SSA'],d['df_A'],d['MSA'],d['F'],d['p_value'])) 42 print( '{:10s}|{:>18,.4f} |{:>8d} |{:>18,.4f} |'.format( '组内（误差）',d['SSE'],d['df_E'],d['MSE'])) 43 print( '{:14s}|{:>18,.4f} |{:>8d} |'.format( '总和',d['SST'],d['df_T'])) 44 print( '-' * maxedg ) 45 print('备注：显著性水平为 {:.2%} 时，F的临界值是 {:.6f}。'.format(d['a'],d['F_alpha'])) 46 47 48 p = 0.95 # 设定置信度水平 49 maxedg = 93 # 设定输出时装饰分隔符的最大长度 50 # 计算并输出单因素方差分析表 51 res = ANOVA_oneway( da, a = 1-p ) 52 print_ANOVA_oneway( res, maxedg = maxedg )
