方差分析
方差分析(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
# ==================================================================================================