方差分析
方差分析(analysis of variance, ANOVA)是利用样本数据检验两个或两个以上总体均值间是否有差异的一种方法。在研究单个变量时，它能够解决多个总体的均值是否相等的检验问题；在研究多个变量对不同总体的影响时，它也是分析各个自变量对因变量影响的方法。
基本原理
方差分析主要是通过方差比较的方式对不同总体参数进行假设检验的。
\(F\) 检验的统计量
\[F = \frac{S_1^2}{S_2^2}=\frac{\sum{(x_{1i}-\bar{x}_1)^2}/\nu_1}{\sum{(x_{2i}-\bar{x}_2)^2}/\nu_2} \]
其中\(\nu_1,\nu_2\) 分别表示总体1和总体2的自由度。
方差是衡量一个总体或样本数据离散程度的重要指标，代表了其所反映数据的差异程度，同时也包含了数据变动的信息。当 \(F=1\) 时，表示2个总体方差相等，没有差异；当 \(F \approx 1\) 时，表示2个总体没有显著差异；当 \(F \neq 1\) 时，表示2个总体有显著差异。
在方差分析中，通常把影响因变量的可控制的定性变量或离散型变量称之为因素，而各个因素具有的表现称之为水平。而影响因变量的定量或连续型变量称之为协变量。
对于多总体均值比较的方差分析问题实际上是一个假设检验问题，其研究本质就是要比较在因素不同水平下的因变量总体均值是否相等，可建立如下线性模型进行分析
\[x_{ij} = \bar{x}_i+\varepsilon_{ij} \]
\(x_{ij}\) 表示作为影响因素的第\(i\)个水平下因变量的第\(j\)个观测值，\(\bar{x}_i\) 表示第\(i\)个水平下因变量的均值；\(\varepsilon_{ij}\) 表示第\(i\)个水平下的因变量第\(j\) 个观测值与该水平下因变量均值之间的残差，也称之为随机扰动项，服从均值为0的一个正态分布。在实际问题中，残差表示除所考虑因素之外的其他因素或不可观测的随机因素（如天气、政策变动、不可抗力因素等）的影响。
对于上述模型可找出造成因变量差异的各种不同来源，并根据方差比较的方法，对这些不同来源的差异进行分析，找出对因变量影响较大的因素及其水平。然后根据合适的参数估计方法对该线性模型进行估计，得出因素水平变动对因变量变动的具体影响。
组间方差：也称为水平之间的方差，即组间离差平方和除以自由度 \((r-1)\)（r为组数或水平数）。也包括随机性因素（因为随机因素的影响在随机抽样过程中不可避免）
组内方差：也称为水平内部的方差，即组内离差平方和除以自由度 \((n-r)\)（n为样本容量）。水平内部防擦好仅包括随机性因素。
如果不同水平对因变量没有影响，那么在水平之间的方差中，就仅仅有随机因素的差异，而没有系统性差异，它与水平内部方差就应该近似，即
\[F = \frac{组间方差}{组内方差} \approx 1 \]
反之，不同水平对结果有显著影响，水平之间的方差就会大于水平内的方差，当这个比值达到某个程度，或者说达到某临界点，就可做出不同水平之间存在着显著差异的判断。
因此，方差分析就是通过不同方差的比较，做出拒绝原假设或不能拒绝原假设的判断。组间方差和组内方差之比时一个统计量，这个统计量在正态假定下服从第一自由度为 \(r-1\)，第二自由度为 \(n-r\) 的 \(F\) 分布。
在对数据进行方差分析时，应该满足如下两个前提条件：
- 各组的观察数据，要看作是来自于正态分布总体的随机样本。该条件通常情况下较容易满足
- 各组的观察数据，是从具有相同方差且相互独立的总体中抽取得到的，即具有同方差性。即由各因素水平所区分的各总体的方差应该相等，其原假设为：\(H_0:\sigma_1^2=\sigma_2^2=...=\sigma_r^2\)。在做方差分析之前应当进行该同方差性检验
方差分析根据所研究的因变量数目不同，可以分为一元方差分析和多元方差分析。在分析一个因变量时，根据影响因素数目不同，又可以进一步细分为单因素方差分析和多因素方差分析；如果在影响因素中具有协变量，称之为协防擦好分析。含有协变量的一元多因素方差分析在实际应用中非常常见
在方差分析过程中，不仅可以对因素是否对因变量产生显著影响进行检验，还可以通过对因素水平的均值进行多重比较来分析因素的哪个或哪些水平对因变量的影响最显著，也可以通过一般化线性模型估计因素的水平对因变量的具体影响。
一元方差分析
当所研究的因变量只有一个时的方差分析就是一元方差分析，但其可受一个或多个因素或协变量的影响
一元单因素
一元单因素方差分析（one-way ANOVA）主要研究单独一个因素对因变量的影响。通过因素的不同水平对因变量进行分组，计算组间和组内方差，利用方差比较的方法对各分组所形成的总体进行均值比较，从而对各总体均值相等的原假设进行检验。
流程
1.样本差异分析
2.方差同质性检验
3.方差来源分解及检验过程
4.多重比较检验
5.方差分析模型的参数估计和预测
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
dc_sales = pd.read_csv('./data/dc_sales.csv')
print(dc_sales.head(5))
# market pixel sales
# 0 1 1 70
# 1 1 2 101
# 2 1 3 114
# 3 1 4 120
# 4 1 5 132
# 探究不同像素下对售价对影响，\alpha=0.05
# 检验成像元器件的像素数对销售量的影响，把每一类不同的像素的数码相机总销售量分别看成不同的总体
# H_0:\mu_1=\mu_2=\mu+3=\mu_4;H_1:\mu_1,\mu_2,\mu_3,\mu_4不完全相等
# 对数据作标签化处理
dc_sales['pixel'] = dc_sales['pixel'].astype('category')
dc_sales['pixel'].cat.categories = ['<500像素', '500-600', '600~800', '800~1000', '>1000']
dc_sales['pixel'].cat.set_categories = ['<500', '500-600', '600~800', '800~1000', '>1000']
# 编制销售数据表
table_res = pd.pivot_table(dc_sales, index=['pixel'], columns=['market'], values=['sales'], aggfunc='sum')
print(table_res)
# sales
# market 1 2 3 4 5 6 7 8
# pixel
# <500 70 67 82 87 80 80 87 96
# 500-600 101 76 97 88 92 99 123 90
# 600~800 114 96 128 103 107 91 99 119
# 800~1000 120 98 132 128 132 132 131 119
# >1000 132 102 123 119 123 135 126 117
# 对数据进行分组
G = dc_sales['pixel'].unique() # G用于统计变量pixel的像素属性
args = [] # 列表args用于存储不同像素属性下的销售数据
for i in list(G):
args.append(dc_sales[dc_sales['pixel'] == i]['sales'])
# 绘制盒须图初步考察因素各水平对因变量的影响
dc_sales_plot = plt.boxplot(args, vert=True, patch_artist=True)
colors = ['pink', 'lightblue', 'lightgreen', 'cyan', 'lightyellow']
for patch, color in zip(dc_sales_plot['boxes'], colors):
patch.set_facecolor(color)
# 指定为黑体中文字体，防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False
fig = plt.gcf()
fig.set_size_inches(8, 5)
combinebox = plt.subplot(111)
combinebox.set_xticklabels(G)
plt.show()
# 显示出不同的像素级别在销售上确实存在着显著区别，但是这种区别仅限于样本数据，对于总体数据，需要进一步判断
# 1.方差同质性检验
# 一元方差分析常用Levene的检验，多元方差分析常用Bartlett的球形检验法
levene_res = stats.levene(*args)
print(levene_res) # LeveneResult(statistic=0.233384556281214, pvalue=0.9176929576341715)
# pvalue比较大，故认为满足方差齐性假设
# 2.方差来源分解及检验过程
one_res = stats.f_oneway(*args)
print(one_res) # F_onewayResult(statistic=19.57176228742291, pvalue=1.5491302153222814e-08)
# pvalue接近0，故拒绝H_0，即认为像素级别对销售有显著影响
# 更详细计算结果
from statsmodels.formula.api import ols
dc_sales_anova = sm.stats.anova_lm(ols('sales ~ C(pixel)', dc_sales).fit())
"""
因变量～自变量列表：是statsmodels模块中表示变量关系的最主要形式
osl()的参数：1：表示模型，2：用于指定所分析的数据
fit()表示对模型进行拟合或估计。
"""
print(dc_sales_anova)
# df sum_sq mean_sq F PR(>F)
# C(pixel) 4.0 10472.850 2618.2125 19.571762 1.549130e-08
# Residual 35.0 4682.125 133.7750 NaN NaN
# 组间离差平方和为10472.85,对应的组间方差为2618.2125=10472.85/4
# 组内离差平方和为468.125,对应的组内方差为133.7750=4682.125/35
# 用于判定组内方差和组间方差的差异的F=19.571762=2618.2125/133.7750
# 对应的p值为1.549130e-08<\alpha
# 3.多重比较检验
# 判定具体哪个水平对观察变量产生了显著影响，这就是单因素方差分析的均值多重比较检验
# 有多种方法：
# LSD法：最灵敏，会犯假阳性错误
# Sidak法：比LSD法保守
# Bonferroni法：比Sidak法更保守，常用
# Scheffe法：多用于进行比较的两组间样本含量不等时
# Dunnet法：常用于多个实验组与一个对照组的比较
# SNK法：寻找同质亚组的方法
# Turkey法：最迟钝，要求各组样本含量相同
# Duncan法：与Sidak法类似
from statsmodels.stats.multicomp import pairwise_tukeyhsd
dc_sales_anova_post = pairwise_tukeyhsd(dc_sales['sales'], dc_sales['pixel'],alpha=0.05)
# 系统自动将不同像素数进行两两比较，并在reject列给出是否应该决绝两组属性没有差异的检验结果
res = dc_sales_anova_post.summary()
print(res)
# Multiple Comparison of Means - Tukey HSD, FWER=0.05
# ==========================================================
# group1 group2 meandiff p-adj lower upper reject
# ----------------------------------------------------------
# 500-600 600~800 11.375 0.3029 -5.2516 28.0016 False
# 500-600 800~1000 28.25 0.001 11.6234 44.8766 True
# 500-600 <500像素 -14.625 0.1072 -31.2516 2.0016 False
# 500-600 >1000 26.375 0.001 9.7484 43.0016 True
# 600~800 800~1000 16.875 0.0453 0.2484 33.5016 True
# 600~800 <500像素 -26.0 0.001 -42.6266 -9.3734 True
# 600~800 >1000 15.0 0.0934 -1.6266 31.6266 False
# 800~1000 <500像素 -42.875 0.001 -59.5016 -26.2484 True
# 800~1000 >1000 -1.875 0.9 -18.5016 14.7516 False
# <500像素 >1000 41.0 0.001 24.3734 57.6266 True
# ----------------------------------------------------------
# <600的相机与中高像素相比，销售量明显较少，且差异最为显著；像素高的相机明显比像素低的销量大
# 4.参数估计和预测
# 方差分析实际上是对一般线性模型进行分析，还可以对用于方差分析的线性模型进行参数估计和假设检验
# 根据参数估计的结果，可以得出因素水平之间的变动状况对因变量产生的具体影响并据此进行模型预测
formula = 'sales~C(pixel)'
dc_sales_est = ols(formula, dc_sales).fit() # dc_sales_est是一个模型对象
print(dc_sales_est.summary2())
# 第1张表主要展示用于模型诊断的总体信息，如拟合度判定系数R^2方、F统计量、P值、AIC和BIC等信息指数等
# 第2张表可知各种像素数对销售量的具体影响：
# Intercept表示<500像素的对因变量的具体影响：在此像素水平下，相机销量为81.125台
# 其他水平对因变量的影响均以截距项为基准来衡量，其对应的参数估计值代表了各个水平对应变量影响与截距项对因变量影响的差距
# Results: Ordinary least squares
# ====================================================================
# Model: OLS Adj. R-squared: 0.656
# Dependent Variable: sales AIC: 314.0202
# Date: 2020-06-04 14:32 BIC: 322.4646
# No. Observations: 40 Log-Likelihood: -152.01
# Df Model: 4 F-statistic: 19.57
# Df Residuals: 35 Prob (F-statistic): 1.55e-08
# R-squared: 0.691 Scale: 133.77
# --------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# --------------------------------------------------------------------
# Intercept 81.1250 4.0892 19.8387 0.0000 72.8234 89.4266
# C(pixel)[T.500-600] 14.6250 5.7831 2.5289 0.0161 2.8848 26.3652
# C(pixel)[T.600~800] 26.0000 5.7831 4.4959 0.0001 14.2598 37.7402
# C(pixel)[T.800~1000] 42.8750 5.7831 7.4139 0.0000 31.1348 54.6152
# C(pixel)[T.>1000] 41.0000 5.7831 7.0897 0.0000 29.2598 52.7402
# --------------------------------------------------------------------
# Omnibus: 0.757 Durbin-Watson: 1.535
# Prob(Omnibus): 0.685 Jarque-Bera (JB): 0.172
# Skew: -0.090 Prob(JB): 0.917
# Kurtosis: 3.266 Condition No.: 6
# ====================================================================
#
# 估计不含截距项模型的参数，直接查看绝对数值
formula = 'sales~C(pixel)-1'
dc_sales_est = ols(formula, dc_sales).fit() # dc_sales_est是一个模型对象
print(dc_sales_est.summary2())
# 不含截距项的参数估计记过代表了因素水平对因变量的绝对影响
# Results: Ordinary least squares
# =====================================================================
# Model: OLS Adj. R-squared: 0.656
# Dependent Variable: sales AIC: 314.0202
# Date: 2020-06-04 16:08 BIC: 322.4646
# No. Observations: 40 Log-Likelihood: -152.01
# Df Model: 4 F-statistic: 19.57
# Df Residuals: 35 Prob (F-statistic): 1.55e-08
# R-squared: 0.691 Scale: 133.78
# ---------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------
# C(pixel)[<500像素] 81.1250 4.0892 19.8387 0.0000 72.8234 89.4266
# C(pixel)[500-600] 95.7500 4.0892 23.4151 0.0000 87.4484 104.0516
# C(pixel)[600~800] 107.1250 4.0892 26.1968 0.0000 98.8234 115.4266
# C(pixel)[800~1000] 124.0000 4.0892 30.3235 0.0000 115.6984 132.3016
# C(pixel)[>1000] 122.1250 4.0892 29.8650 0.0000 113.8234 130.4266
# ---------------------------------------------------------------------
# Omnibus: 0.757 Durbin-Watson: 1.535
# Prob(Omnibus): 0.685 Jarque-Bera (JB): 0.172
# Skew: -0.090 Prob(JB): 0.917
# Kurtosis: 3.266 Condition No.: 1
# =====================================================================
#
# 5.方差分析模型的预测
# 可以使用模型的参数估计值对因变量进行预测，确保预测较准确的前提是估计出的模型得依据统计理论进行模型诊断。
# 本例中R^2和F值均较大，可以预测
# 方法一：
print(dc_sales_est.fittedvalues)
# 0 81.125
# 1 95.750
# 3 124.000
# 4 122.125
# ...
# 38 124.000
# 39 122.125
# dtype: float64
# 方法二
dc_sales_influence = dc_sales_est.get_influence()
print(dc_sales_influence.summary_table())
# ==================================================================================================
# obs endog fitted Cook's student. hat diag dffits ext.stud. dffits
# value d residual internal residual
# --------------------------------------------------------------------------------------------------
# 0 70.000 81.125 0.030 -1.028 0.125 -0.389 -1.029 -0.389
# 1 101.000 95.750 0.007 0.485 0.125 0.183 0.480 0.181
# 2 114.000 107.125 0.012 0.635 0.125 0.240 0.630 0.238
# 3 120.000 124.000 0.004 -0.370 0.125 -0.140 -0.365 -0.138
# 4 132.000 122.125 0.024 0.913 0.125 0.345 0.911 0.344
# 5 67.000 81.125 0.049 -1.306 0.125 -0.493 -1.319 -0.499
# 6 76.000 95.750 0.095 -1.825 0.125 -0.690 -1.892 -0.199
# ...
# 37 119.000 107.125 0.034 1.098 0.125 0.415 1.101 0.416
# 38 119.000 124.000 0.006 -0.462 0.125 -0.175 -0.457 -0.173
# 39 117.000 122.125 0.006 -0.474 0.125 -0.179 -0.468 -0.177
# ==================================================================================================
一元多因素
当有两个或者两个以上的因素对因变量产生影响时，可以用多因素方差分析的方法来进行分析。多因素方差分析的原理与单因素方差分析基本一致，也是利用方差分析比较的方法，通过假设检验的过程来判断多个因素是否对因变量产生显著性影响。
在多因素方差分析中，由于影响因变量的因素有多个，其中某些因素除了自身对因变量产生影响之外，它们之间也有可能会共同对因变量产生影响。在多因素方差分析中，把因素单独对因变量产生的影响称之为主效应；把因素之间共同对因变量产生的影响，或者因素某些水平同时出现时，除了主效应之外的附加影响，称之为交互效应
多因素方差分析不仅要考虑每个因素的主效应，往往还要考虑因素之间的交互效应。此外，多因素方差分析往往假定因素与因变量之间的关系是线性关系。从这个方面来说，方差分析的模型是如下一个一般化线性模型的延续
因变量=因素1主效应+因素2主效应+...+因素n主效应+因素交互效应1+因素交互效应2+...+因素交互效应m+随机误差
多因素方差分析往往选用一般化线性模型进行参数估计。
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 学历、购房者所在单位类型、收入水平和户型对购房面积影响显著，显著性水平\alpha=0.05
house = pd.read_csv('./data/house.csv')
print(house.head(5))
# education unit income type space
# 0 1 3 2 2 75.0
# 1 1 5 1 6 55.0
# 2 3 1 2 8 56.0
# 3 2 6 2 4 51.0
# 4 2 4 1 5 60.0
# 数据预处理：加标签
house['education'] = house['education'].astype('category')
house['education'].cat.categories = ['初中及以下', '高中(中专)', '大学', '研究生及以上']
house['education'].cat.set_categories = ['初中及以下', '高中(中专)', '大学', '研究生及以上']
house['unit'] = house['unit'].astype('category')
house['unit'].cat.categories = ['国营企业', '行政事业单位', '大专院校科研单位', '私营企业', '失业', '其他']
house['unit'].cat.set_categories = ['国营企业', '行政事业单位', '大专院校科研单位', '私营企业', '失业', '其他']
house['income'] = house['income'].astype('category')
house['income'].cat.categories = ['<10000', '10000~25000', '25000~50000', '50000~75000', '>75000']
house['income'].cat.set_categories = ['<10000', '10000~25000', '25000~50000', '50000~75000', '>75000']
house['type'] = house['type'].astype('category')
house['type'].cat.categories = ['一室一厅', '二室一厅', '二室二厅', '三室一厅', '三室二厅', '三室三厅',
'四室二厅一卫', '四室二厅二卫', '四室三厅一卫', '四室三厅二卫', '更大户型']
house['type'].cat.set_categories = ['一室一厅', '二室一厅', '二室二厅', '三室一厅', '三室二厅', '三室三厅',
'四室二厅一卫', '四室二厅二卫', '四室三厅一卫', '四室三厅二卫', '更大户型']
# 1.只考虑主效应
# a.方差来源分解及检验过程
from statsmodels.formula.api import ols
formula = 'space ~ C(education)+C(unit)+C(income)+C(type)'
house_anova = sm.stats.anova_lm(ols(formula, data=house).fit(), typ=3)
"""
typ=3表示做方差分析type III型检验
typ III即平方和分解法，适用于平衡的ANOVA模型和非平衡的ANOVA模型，凡适用typ I和typ II的模型均可用该法
参数默认适用typ I
"""
print(house_anova)
# sum_sq df F PR(>F)
# Intercept 28663.946440 1.0 89.701281 1.591592e-19
# C(education) 1519.401986 3.0 1.584945 1.922550e-01
# C(unit) 886.371451 5.0 0.554764 7.347011e-01
# C(income) 10545.816026 4.0 8.250549 1.990132e-06
# C(type) 9604.428837 10.0 3.005621 1.093296e-03
# Residual 143477.460152 449.0 NaN NaN
# 由于education和unint因素主效应的影响十分不显著，所以剔除这两个因素
formula = 'space ~ C(income)+C(type)'
house_anova = sm.stats.anova_lm(ols(formula, data=house).fit(), typ=3)
print(house_anova)
# sum_sq df F PR(>F)
# Intercept 29277.490099 1.0 91.609792 6.682452e-20
# C(income) 12033.242560 4.0 9.413058 2.558367e-07
# C(type) 10553.485226 10.0 3.302204 3.765261e-04
# Residual 146052.214252 457.0 NaN NaN
# b.对因素进行多重比较检验
# 因素的哪些水平对因变量的影响最大
from statsmodels.stats.multicomp import pairwise_tukeyhsd
house_anova_post = pairwise_tukeyhsd(house['space'], house['income'], alpha=0.05)
res = house_anova_post.summary()
print(res)
# Multiple Comparison of Means - Tukey HSD, FWER=0.05
# ===============================================================
# group1 group2 meandiff p-adj lower upper reject
# ---------------------------------------------------------------
# 10000~25000 25000~50000 7.6175 0.0042 1.7029 13.5321 True
# 10000~25000 50000~75000 12.1263 0.1151 -1.665 25.9176 False
# 10000~25000 <10000 -4.6062 0.1992 -10.4621 1.2497 False
# 10000~25000 >75000 31.2777 0.001 15.812 46.7435 True
# 25000~50000 50000~75000 4.5088 0.9 -9.7758 18.7933 False
# 25000~50000 <10000 -12.2238 0.001 -19.1621 -5.2854 True
# 25000~50000 >75000 23.6602 0.001 7.753 39.5674 True
# 50000~75000 <10000 -16.7325 0.0122 -30.9929 -2.4722 True
# 50000~75000 >75000 19.1514 0.073 -1.0539 39.3567 False
# <10000 >75000 35.884 0.001 19.9985 51.7694 True
# ---------------------------------------------------------------
# 可以认为，低收入群体和高收入群体对购房面积最敏感，而中收入群体对购房面积敏感性不如高手如何低收入群体。
house_anova_post = pairwise_tukeyhsd(house['space'], house['type'], alpha=0.05)
res = house_anova_post.summary()
print(res)
# Multiple Comparison of Means - Tukey HSD, FWER=0.05
# =======================================================
# group1 group2 meandiff p-adj lower upper reject
# -------------------------------------------------------
# 一室一厅 三室一厅 -20.4681 0.5215 -50.9333 9.9971 False
# 一室一厅 三室三厅 -25.75 0.8722 -77.6303 26.1303 False
# 一室一厅 三室二厅 -15.3197 0.8552 -45.667 15.0275 False
# 一室一厅 二室一厅 -21.5359 0.4602 -52.238 9.1661 False
# 一室一厅 二室二厅 -20.0805 0.5676 -51.0323 10.8713 False
# 一室一厅 四室三厅一卫 -36.25 0.7825 -103.2272 30.7272 False
# 一室一厅 四室三厅二卫 -22.09 0.7656 -62.2763 18.0963 False
# 一室一厅 四室二厅一卫 -3.2714 0.9 -37.2351 30.6922 False
# 一室一厅 四室二厅二卫 -6.8017 0.9 -38.7539 25.1504 False
# 一室一厅 更大户型 6.6667 0.9 -27.9202 41.2536 False
# 三室一厅 三室三厅 -5.2819 0.9 -48.0056 37.4418 False
# 三室一厅 三室二厅 5.1484 0.4717 -2.2479 12.5446 False
# 三室一厅 二室一厅 -1.0678 0.9 -9.8065 7.6709 False
# 三室一厅 二室二厅 0.3876 0.9 -9.1918 9.9669 False
# 三室一厅 四室三厅一卫 -15.7819 0.9 -75.9458 44.382 False
# 三室一厅 四室三厅二卫 -1.6219 0.9 -28.9841 25.7403 False
# 三室一厅 四室二厅一卫 17.1967 0.0434 0.2474 34.1459 True
# 三室一厅 四室二厅二卫 13.6664 0.018 1.229 26.1037 True
# 三室一厅 更大户型 27.1348 0.001 8.9688 45.3007 True
# 三室三厅 三室二厅 10.4303 0.9 -32.2095 53.07 False
# 三室三厅 二室一厅 4.2141 0.9 -38.6789 47.107 False
# 三室三厅 二室二厅 5.6695 0.9 -37.4026 48.7416 False
# 三室三厅 四室三厅一卫 -10.5 0.9 -83.8699 62.8699 False
# 三室三厅 四室三厅二卫 3.66 0.9 -46.4612 53.7812 False
# 三室三厅 四室二厅一卫 22.4786 0.8721 -22.8063 67.7634 False
# 三室三厅 四室二厅二卫 18.9483 0.9 -24.8482 62.7447 False
# 三室三厅 更大户型 32.4167 0.4447 -13.3375 78.1708 False
# 三室二厅 二室一厅 -6.2162 0.359 -14.5345 2.1021 False
# 三室二厅 二室二厅 -4.7608 0.8297 -13.9582 4.4367 False
# 三室二厅 四室三厅一卫 -20.9303 0.9 -81.0346 39.174 False
# 三室二厅 四室三厅二卫 -6.7703 0.9 -34.0011 20.4606 False
# 三室二厅 四室二厅一卫 12.0483 0.4195 -4.6881 28.7847 False
# 三室二厅 四室二厅二卫 8.518 0.4604 -3.6276 20.6637 False
# 三室二厅 更大户型 21.9864 0.0042 4.0189 39.9539 True
# 二室一厅 二室二厅 1.4554 0.9 -8.8525 11.7634 False
# 二室一厅 四室三厅一卫 -14.7141 0.9 -74.9983 45.5702 False
# 二室一厅 四室三厅二卫 -0.5541 0.9 -28.1797 27.0716 False
# 二室一厅 四室二厅一卫 18.2645 0.0299 0.8931 35.636 True
# 二室一厅 四室二厅二卫 14.7342 0.0123 1.7274 27.7411 True
# 二室一厅 更大户型 28.2026 0.001 9.6422 46.7631 True
# 二室二厅 四室三厅一卫 -16.1695 0.9 -76.5813 44.2423 False
# 二室二厅 四室三厅二卫 -2.0095 0.9 -29.9125 25.8935 False
# 二室二厅 四室二厅一卫 16.8091 0.0846 -1.0001 34.6182 False
# 二室二厅 四室二厅二卫 13.2788 0.0622 -0.3071 26.8647 False
# 二室二厅 更大户型 26.7472 0.001 7.7764 45.7179 True
# 四室三厅一卫 四室三厅二卫 14.16 0.9 -51.464 79.784 False
# 四室三厅一卫 四室二厅一卫 32.9786 0.8012 -29.0303 94.9874 False
# 四室三厅一卫 四室二厅二卫 29.4483 0.8982 -31.4821 90.3786 False
# 四室三厅一卫 更大户型 42.9167 0.4882 -19.4357 105.2691 False
# 四室三厅二卫 四室二厅一卫 18.8186 0.6591 -12.3919 50.029 False
# 四室三厅二卫 四室二厅二卫 15.2883 0.8109 -13.7204 44.2969 False
# 四室三厅二卫 更大户型 28.7567 0.1213 -3.1309 60.6442 False
# 四室二厅一卫 四室二厅二卫 -3.5303 0.9 -23.0262 15.9656 False
# 四室二厅一卫 更大户型 9.9381 0.9 -13.6289 33.5051 False
# 四室二厅二卫 更大户型 13.4684 0.5551 -7.094 34.0308 False
# -------------------------------------------------------
# 三室一厅和两室一厅比较显著
# 结论：
# 收入和户型对消费者在拟购房面积决策上的影响非常显著，而学历和工作单位对该决策影响不显著
# 其中，低收入群体和高收入群体对购房面积的要求比较敏感，实用型的户型对面积影响较大
# c.利用方差分析模型进行参数估计
house_anova_est = ols(formula, data=house).fit()
res = house_anova_est.summary2()
# print(res)
# Results: Ordinary least squares
# ===========================================================================
# Model: OLS Adj. R-squared: 0.146
# Dependent Variable: space AIC: 4076.2755
# Date: 2020-06-05 15:13 BIC: 4138.6302
# No. Observations: 472 Log-Likelihood: -2023.1
# Df Model: 14 F-statistic: 6.752
# Df Residuals: 457 Prob (F-statistic): 1.38e-12
# R-squared: 0.171 Scale: 319.59
# ---------------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------------
# Intercept 86.1464 9.0005 9.5713 0.0000 68.4589 103.8339
# C(income)[T.10000~25000] 4.2072 2.1086 1.9953 0.0466 0.0635 8.3509
# C(income)[T.25000~50000] 9.9601 2.5641 3.8844 0.0001 4.9212 14.9990
# C(income)[T.50000~75000] 16.1291 5.1202 3.1501 0.0017 6.0670 26.1911
# C(income)[T.>75000] 29.3518 5.9372 4.9437 0.0000 17.6842 41.0193
# C(type)[T.二室一厅] -23.0808 9.1696 -2.5171 0.0122 -41.1005 -5.0610
# C(type)[T.二室二厅] -22.5007 9.2538 -2.4315 0.0154 -40.6859 -4.3155
# C(type)[T.三室一厅] -23.0634 9.1100 -2.5317 0.0117 -40.9660 -5.1608
# C(type)[T.三室二厅] -19.4782 9.0933 -2.1420 0.0327 -37.3481 -1.6083
# C(type)[T.三室三厅] -28.6264 15.5200 -1.8445 0.0658 -59.1259 1.8730
# C(type)[T.四室二厅一卫] -7.2395 10.1864 -0.7107 0.4776 -27.2575 12.7784
# C(type)[T.四室二厅二卫] -12.6134 9.5902 -1.3152 0.1891 -31.4598 6.2329
# C(type)[T.四室三厅一卫] -38.3536 20.0149 -1.9162 0.0560 -77.6864 0.9791
# C(type)[T.四室三厅二卫] -26.4948 12.0484 -2.1990 0.0284 -50.1719 -2.8176
# C(type)[T.更大户型] -4.5995 10.5095 -0.4377 0.6618 -25.2524 16.0533
# ---------------------------------------------------------------------------
# Omnibus: 88.385 Durbin-Watson: 1.372
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 157.668
# Skew: 1.083 Prob(JB): 0.000
# Kurtosis: 4.825 Condition No.: 47
# ===========================================================================
# 截距为income<10000的且户型为一室一厅
# 非截距
house_anova_est = ols('space ~ C(income)+C(type)-1', data=house).fit()
res = house_anova_est.summary2()
print(res)
# Results: Ordinary least squares
# =========================================================================
# Model: OLS Adj. R-squared: 0.146
# Dependent Variable: space AIC: 4076.2755
# Date: 2020-06-05 15:42 BIC: 4138.6302
# No. Observations: 472 Log-Likelihood: -2023.1
# Df Model: 14 F-statistic: 6.752
# Df Residuals: 457 Prob (F-statistic): 1.38e-12
# R-squared: 0.171 Scale: 319.59
# -------------------------------------------------------------------------
# Coef. Std.Err. t P>|t| [0.025 0.975]
# -------------------------------------------------------------------------
# C(income)[<10000] 86.1464 9.0005 9.5713 0.0000 68.4589 103.8339
# C(income)[10000~25000] 90.3536 9.0005 10.0387 0.0000 72.6661 108.0411
# C(income)[25000~50000] 96.1065 9.1886 10.4593 0.0000 78.0492 114.1637
# C(income)[50000~75000] 102.2754 10.2017 10.0253 0.0000 82.2273 122.3236
# C(income)[>75000] 115.4981 10.6282 10.8671 0.0000 94.6119 136.3844
# C(type)[T.二室一厅] -23.0808 9.1696 -2.5171 0.0122 -41.1005 -5.0610
# C(type)[T.二室二厅] -22.5007 9.2538 -2.4315 0.0154 -40.6859 -4.3155
# C(type)[T.三室一厅] -23.0634 9.1100 -2.5317 0.0117 -40.9660 -5.1608
# C(type)[T.三室二厅] -19.4782 9.0933 -2.1420 0.0327 -37.3481 -1.6083
# C(type)[T.三室三厅] -28.6264 15.5200 -1.8445 0.0658 -59.1259 1.8730
# C(type)[T.四室二厅一卫] -7.2395 10.1864 -0.7107 0.4776 -27.2575 12.7784
# C(type)[T.四室二厅二卫] -12.6134 9.5902 -1.3152 0.1891 -31.4598 6.2329
# C(type)[T.四室三厅一卫] -38.3536 20.0149 -1.9162 0.0560 -77.6864 0.9791
# C(type)[T.四室三厅二卫] -26.4948 12.0484 -2.1990 0.0284 -50.1719 -2.8176
# C(type)[T.更大户型] -4.5995 10.5095 -0.4377 0.6618 -25.2524 16.0533
# -------------------------------------------------------------------------
# Omnibus: 88.385 Durbin-Watson: 1.372
# Prob(Omnibus): 0.000 Jarque-Bera (JB): 157.668
# Skew: 1.083 Prob(JB): 0.000
# Kurtosis: 4.825 Condition No.: 35
# =========================================================================
# 本例模型的拟合优度0.171偏低，但是对于实际截面数据而言，该模型拟合程度勉强能接受。
# 如，考察年收入在50000-75000元，购买意向户型为三室二厅的拟购面积为：102.2754-19.4782（平米）
# 2.存在交互效应的多因素方差分析
# 当不确定因素之间是否存在交互效应时，可以首先考虑利用方差分析的群模型进行分析，即考虑主效应和交互效应。
formula = 'space ~ C(income)*C(type)'
"""
在statsmodels中，“：”表示因素之间的交互效应，“ * ”表示因素之间的全部效应，即变量间的交互效应+参加交互的变量的各自主效应。
即 'space ~ C(income)*C(type)' 等价于 'space ~ C(income)+C(type)+C(income):C(type)'
"""
house_anova_inter = sm.stats.anova_lm(ols(formula, data=house).fit())
print(house_anova_inter)
# df sum_sq mean_sq F PR(>F)
# C(income) 4.0 19655.107559 4913.776890 16.629186 1.077135e-12
# C(type) 10.0 10553.485226 1055.348523 3.571507 1.438294e-04
# C(income):C(type) 40.0 20470.983957 511.774599 1.731946 4.691272e-03
# Residual 436.0 128834.132740 295.491130 NaN NaN
r = ols(formula, data=house).fit().rsquared
print(r) # 拟合优度系数R^2=0.269>0.171
# 主效应及交互效应的p值均较小，同时R^2比之前大，故各因素的主效应与交互效应均非常显著
# 交互效应图
from statsmodels.graphics.api import interaction_plot
# 指定为黑体中文字体，防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False
plt.figure(figsize=(12, 6))
fig = interaction_plot(np.array(house['income']), np.array(house['type']), np.array(house['space']), ax=plt.gca())
fig_adj = plt.subplot(111)
plt.legend(prop={'family': 'Heiti TC', 'size': 10.5}, loc='upper left', frameon=False)
fig_adj.set_xticklabels(house['income'].unique())
plt.show()
协方差分析
将那些难以控制的因素当做协变量，在排除协变量影响的条件下，分析可控因素对因变量的影响，从而更加准确地对可控因素进行评价。协变量之间没有交互效应，且与因素变量之间也没有交互效应。
考虑协变量的方差分析模型的一般形式如下：
因变量 = 因数主效应+因素间交互效应+协变量+随机变量
对于该模型，同样可以用一般化线性模型的形式进行方差分析
示例
import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# 研究卖场、保修、提成对笔记本销售影响情况，显著性水平\alpha=0.05
sale_points = pd.read_csv('./data/sale_points.csv')
print(sale_points.head(5))
# market warranty sales points
# 0 1 1 26.0 1.8
# 1 1 1 22.0 1.1
# 2 1 1 21.8 0.9
# 3 1 1 33.1 2.2
# 4 2 1 22.0 2.0
# sales为因变量，因素：market，warranty，协变量：points
# 数据预处理：加标签
sale_points['market'] = sale_points['market'].astype('category')
sale_points['market'].cat.categories = ['market1', 'market2', 'market3']
sale_points['market'].cat.set_categories = ['market1', 'market2', 'market3']
sale_points['warranty'] = sale_points['warranty'].astype('category')
sale_points['warranty'].cat.categories = ['1 year', '3 years']
sale_points['warranty'].cat.set_categories = ['1 year', '3 years']
# 1. 协方差分析
from statsmodels.formula.api import ols
# 使用statsmodels进行协方差分析，在定义formula的时候列示出协变量即可
formula = 'sales ~ points+C(market)*C(warranty)'
sale_points_anova_cov = sm.stats.anova_lm(ols(formula, data=sale_points).fit())
print(sale_points_anova_cov)
# df sum_sq mean_sq F PR(>F)
# C(market) 2.0 593.160833 296.580417 56.984051 2.903413e-08
# C(warranty) 1.0 512.450417 512.450417 98.460650 1.734504e-08
# C(market):C(warranty) 2.0 167.155833 83.577917 16.058404 1.211601e-04
# points 1.0 196.523934 196.523934 37.759505 1.079351e-05
# Residual 17.0 88.478566 5.204622 NaN NaN
# 个效应对因变量的影响均非常显著
# 2. 多重比较检验
from statsmodels.stats.multicomp import pairwise_tukeyhsd
market_anova_post = pairwise_tukeyhsd(sale_points['sales'], sale_points['market'], alpha=0.05)
res = market_anova_post.summary()
print(res)
# Multiple Comparison of Means - Tukey HSD, FWER=0.05
# =======================================================
# group1 group2 meandiff p-adj lower upper reject
# -------------------------------------------------------
# market1 market2 -4.8625 0.3423 -13.4016 3.6766 False
# market1 market3 7.2375 0.1065 -1.3016 15.7766 False
# market2 market3 12.1 0.0049 3.5609 20.6391 True
# -------------------------------------------------------
warranty__anova_post = pairwise_tukeyhsd(sale_points['sales'], sale_points['warranty'], alpha=0.05)
res = warranty__anova_post.summary()
print(res)
# Multiple Comparison of Means - Tukey HSD, FWER=0.05
# ====================================================
# group1 group2 meandiff p-adj lower upper reject
# ----------------------------------------------------
# 1 year 3 years 9.2417 0.0034 3.4056 15.0777 True
# ----------------------------------------------------
# 可以看到market2和market3对销售额有显著影响，一年担保和三年担保对销售额有显著影响
# 3. 参数估计
sale_points_anova_cov_est = ols(formula, data=sale_points).fit()
res = sale_points_anova_cov_est.summary()
print(res)
# OLS Regression Results
# ==============================================================================
# Dep. Variable: sales R-squared: 0.943
# Model: OLS Adj. R-squared: 0.923
# Method: Least Squares F-statistic: 47.05
# Date: Fri, 05 Jun 2020 Prob (F-statistic): 1.16e-09
# Time: 16:36:35 Log-Likelihood: -49.711
# No. Observations: 24 AIC: 113.4
# Df Residuals: 17 BIC: 121.7
# Df Model: 6
# Covariance Type: nonrobust
# ===============================================================================================================
# coef std err t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------------------------------------------------
# Intercept 12.8441 2.386 5.382 0.000 7.809 17.879
# C(market)[T.market2] -8.0349 1.707 -4.706 0.000 -11.637 -4.433
# C(market)[T.market3] 1.3456 1.615 0.833 0.416 -2.061 4.752
# C(warranty)[T.3 years] 3.0485 1.673 1.822 0.086 -0.481 6.578
# C(market)[T.market2]:C(warranty)[T.3 years] 5.4217 2.478 2.188 0.043 0.193 10.650
# C(market)[T.market3]:C(warranty)[T.3 years] 14.0594 2.338 6.014 0.000 9.127 18.991
# points 8.5873 1.397 6.145 0.000 5.639 11.536
# ==============================================================================
# Omnibus: 0.619 Durbin-Watson: 3.022
# Prob(Omnibus): 0.734 Jarque-Bera (JB): 0.690
# Skew: -0.305 Prob(JB): 0.708
# Kurtosis: 2.435 Cond. No. 16.5
# ==============================================================================
#
sale_points_anova_cov_est = ols('sales ~ points+C(market)*C(warranty)-1', data=sale_points).fit()
res = sale_points_anova_cov_est.summary()
print(res)
# OLS Regression Results
# ==============================================================================
# Dep. Variable: sales R-squared: 0.943
# Model: OLS Adj. R-squared: 0.923
# Method: Least Squares F-statistic: 47.05
# Date: Fri, 05 Jun 2020 Prob (F-statistic): 1.16e-09
# Time: 16:39:34 Log-Likelihood: -49.711
# No. Observations: 24 AIC: 113.4
# Df Residuals: 17 BIC: 121.7
# Df Model: 6
# Covariance Type: nonrobust
# ===============================================================================================================
# coef std err t P>|t| [0.025 0.975]
# ---------------------------------------------------------------------------------------------------------------
# C(market)[market1] 12.8441 2.386 5.382 0.000 7.809 17.879
# C(market)[market2] 4.8092 2.890 1.664 0.114 -1.288 10.906
# C(market)[market3] 14.1897 2.448 5.796 0.000 9.025 19.355
# C(warranty)[T.3 years] 3.0485 1.673 1.822 0.086 -0.481 6.578
# C(market)[T.market2]:C(warranty)[T.3 years] 5.4217 2.478 2.188 0.043 0.193 10.650
# C(market)[T.market3]:C(warranty)[T.3 years] 14.0594 2.338 6.014 0.000 9.127 18.991
# points 8.5873 1.397 6.145 0.000 5.639 11.536
# ==============================================================================
# Omnibus: 0.619 Durbin-Watson: 3.022
# Prob(Omnibus): 0.734 Jarque-Bera (JB): 0.690
# Skew: -0.305 Prob(JB): 0.708
# Kurtosis: 2.435 Cond. No. 18.5
# ==============================================================================
# 担保三年，market3的笔记本电脑下的拟销售额为：14.1897+3.0485+14.0594+8.5873
# 4.预测
sale_points_influence = sale_points_anova_cov_est.get_influence()
print(sale_points_influence.summary_table())
# obs endog fitted Cook's student. hat diag dffits ext.stud. dffits
# value d residual internal residual
# --------------------------------------------------------------------------------------------------
# 0 26.000 28.301 0.080 -1.192 0.284 -0.750 -1.208 -0.760
# 1 22.000 22.290 0.002 -0.153 0.310 -0.103 -0.149 -0.100
# 2 21.800 20.573 0.042 0.686 0.385 0.543 0.675 0.534
# 3 33.100 31.736 0.069 0.795 0.434 0.696 0.786 0.688
# 4 22.000 21.984 0.000 0.008 0.254 0.005 0.008 0.005
# 5 19.000 17.690 0.031 0.691 0.310 0.463 0.680 0.456
# 6 17.500 21.984 0.251 -2.275 0.254 -1.327 -2.647 -1.543
# 7 26.000 22.842 0.134 1.614 0.265 0.969 1.702 1.022
# 8 23.000 24.494 0.037 -0.781 0.296 -0.506 -0.771 -0.500
# 9 25.000 25.353 0.002 -0.182 0.273 -0.111 -0.176 -0.108
# 10 32.000 30.506 0.037 0.781 0.296 0.506 0.771 0.500
# 11 30.000 29.647 0.002 0.182 0.273 0.111 0.176 0.108
# 12 36.000 33.067 0.114 1.497 0.262 0.893 1.559 0.930
# 13 32.000 34.355 0.088 -1.227 0.291 -0.787 -1.246 -0.799
# 14 28.000 26.283 0.084 0.962 0.388 0.767 0.960 0.765
# 15 30.000 32.294 0.066 -1.164 0.253 -0.678 -1.177 -0.685
# 16 28.000 26.160 0.041 0.931 0.250 0.538 0.927 0.536
# 17 23.000 23.584 0.005 -0.304 0.290 -0.194 -0.296 -0.189
# 18 24.500 27.019 0.079 -1.277 0.252 -0.741 -1.303 -0.756
# 19 30.000 28.736 0.023 0.652 0.278 0.405 0.641 0.398
# 20 41.000 41.602 0.006 -0.312 0.284 -0.197 -0.304 -0.191
# 21 46.000 46.841 0.011 -0.436 0.285 -0.276 -0.425 -0.269
# 22 48.500 45.896 0.091 1.331 0.265 0.798 1.364 0.818
# 23 41.300 42.461 0.018 -0.594 0.265 -0.357 -0.582 -0.350
# ==================================================================================================