一、概述
（F檢驗）顯著性檢驗：檢測自變量是否真正影響到因變量的波動。
（t檢驗）回歸系數檢驗：單個自變量在模型中是否有效。
二、回歸模型檢驗
檢驗回歸模型的好壞常用的是F檢驗和t檢驗。F檢驗驗證的是偏回歸系數是否不全為0（或全為0），t檢驗驗證的是單個自變量是否對因變量的影響是顯著的（或不顯著）。
F檢驗和t檢驗步驟：
- 提出問題的原假設和備擇假設
- 在原假設的條件下，構造統計量
- 根據樣本信息，計算統計量的值
- 對比統計量的值和理論F分布的值，計算統計量的值超過理論值，則拒絕原假設，否則接受原假設
待檢驗數據集先構造成向量的模式：
可表示成如下形式：
其中β 為n×1的一維向量。
1) 假設
F檢驗假設：
t檢驗假設：
H0為原假設，H1為備擇假設。F檢驗拒絕原假設的條件為計算的F檢驗的值大於查到的理論F值。t檢驗可以通過P值和擬合優度判斷變量的顯著性及模型組合的優劣。
2）計算過程
F檢驗計算過程：
上圖為假設其中一個點所在的平面，由以上點計算出ESS(誤差平方和)，RSS（回歸離差平方和）,TSS（總的離差平方和）。
計算公式為：
其中ESS和RSS都會隨着模型的變化而發生變化（估計值變動）。而TSS衡量的是因變量和均值之間的離差平方和，不會隨着模型的變化而變化。
由以上公式構造F統計量：
（n為數據集向量的行數，p為列數（自變量的個數），n為離差平方和RSS的自由度，n-p-1為誤差平方和ESS的自由度），模型擬合越好，ESS越小，RSS越大，F值越大。
3）python中計算F值得方法：可以直接通過model.fvalue得到f值或者通過計算得到。
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
|
from
sklearn
import
model_selection
import
statsmodels.api as sm
import
pandas as pd
import
numpy as np
import
matplotlib.pyplot as plt
from
sklearn
import
preprocessing
Profit
=
pd.read_excel(r
'Predict to Profit.xlsx'
)
Profit.head()
dummies
=
pd.get_dummies(Profit.State,prefix
=
'State'
)
Profit_New
=
pd.concat([Profit,dummies],axis
=
1
)
Profit_New.drop(labels
=
[
'State'
,
'State_New York'
],axis
=
1
,inplace
=
True
)
train , test
=
model_selection.train_test_split(Profit_New,test_size
=
0.2
,random_state
=
1234
)
model
=
sm.formula.ols(
'Profit~RD_Spend+Administration+Marketing_Spend+State_California+State_Florida'
,data
=
train).fit()
ybar
=
train.Profit.mean()
p
=
model.df_model
#自變量個數
n
=
train.shape[
0
]
#行數，觀測個數
RSS
=
np.
sum
((model.fittedvalues
-
ybar)
*
*
2
)
#計算離差平方和 估計值model.fittedvalues
ESS
=
np.
sum
(model.resid
*
*
2
)
#計算誤差平方和，誤差表示model.resid
F
=
RSS
/
p
/
(ESS
/
(n
-
p
-
1
))
print
( F, model.fvalue)
|
結果如下：
求理論F值，使用Scipy庫子模塊stats中f.ppf
|
1
2
3
|
from
scipy.stats
import
f
F_Theroy
=
f.ppf(q
=
0.95
,dfn
=
p,dfd
=
n
-
p
-
1
)
F_Theroy
|
由於計算的F值遠遠大於理論F值，所以拒絕原假設，證明多元線性回歸方程是顯著的，偏回歸系數不全為0，即所有自變量聯合起來的組合確實對利潤有顯著性影響。
4) t檢驗手工計算比較復雜，套用python的驗計算方式，model.summary()
參數：R-squared為擬合優度（反映模型對樣本的擬合程度）,Adj.R-squared為修正的可決系數。
通過t檢驗得出，只有RD_Spend的P值小於0.05，其余變量均大於0.05，說明其余變量不是影響利潤Profit的重要因素。
三、回歸模型診斷
實際應用中，因變量若為數值型變量，可以考慮線性回歸模型。模型需滿足如下假設：
因變量（殘差）服從正態分布，自變量之間不存在多重共線性，自變量和因變量之間存在線性關系，用於建模的數據不存在異常值，殘差項滿足方差異性和獨立性。
- 正態性檢驗（直方圖，pp圖或qq圖，shapiro檢驗和K-S檢驗）
- 多重共線性檢驗（VIF>10，存在多重共線性，>100則變量間存在嚴重的多重共線性）
- 線性相關性（數據框的corrwith方法或者seaborn模塊的pairplot函數作圖觀察）|p|~>=0.8 高度相關 ，0.5《 |p|《0.8 中度相關，0.3《 |p|《0.5 弱相關，|p|《0.3 幾乎不相關
- 異常值檢驗（帽子矩陣、DFFITS准則、學生化殘差、Cook距離）
- 殘差獨立性檢驗（summary()方法觀察Durbin-Watson值），DW值在2左右，則殘差項之間不相關，偏離2較遠，則不滿足殘差獨立性假設
- 方差齊性檢驗（圖形法和BP檢驗）
具體方法過程：
1）正態性檢測（模型的前提假設是殘差服從（0，1）正態分布，實質上由方程y=Xβ+ε，因變量也應該服從正態分布），檢測的是因變量是否符合正態分布
- 直方圖 scipy.stats seaborn.distplot
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
import
seaborn as sns
import
scipy.stats as stats
from
pylab
import
*
mpl.rcParams[
'font.sans-serif'
]
=
[
'SimHei'
]
plt.rcParams[
'axes.unicode_minus'
]
=
False
sns.distplot(a
=
Profit_New.Profit,
bins
=
10
,
fit
=
stats.norm,
norm_hist
=
True
,
hist_kws
=
{
'color'
:
'green'
,
'edgecolor'
:
'black'
},
kde_kws
=
{
'color'
:
'black'
,
'linestyle'
:
'--'
,
'label'
:
'核密度曲線'
},
fit_kws
=
{
'color'
:
'red'
,
'linestyle'
:
':'
,
'label'
:
'正態密度曲線'
}
)
plt.legend()
plt.show()
|
通過觀察核密度曲線和理論正態密度曲線的擬合程度，2條曲線近似吻合，可認為因變量Profit服從正態分布。
- PP圖和QQ圖 ppplot qqplot
|
1
2
3
4
5
6
|
pq_plot
=
sm.ProbPlot(Profit_New.Profit)
pq_plot.ppplot(line
=
'45'
)
plt.title(
'PP圖'
)
pq_plot.qqplot(line
=
'q'
)
plt.title(
'QQ圖'
)
plt.show()
|
觀察到散點均均勻的分布在直線上，說明因變量近似服從正態分布。
- Shapiro檢驗和K-S檢驗（因變量數據量>5000 use K-S檢驗法，<5000 use Shapiro檢驗）
|
1
|
stats.shapiro(Profit_New.Profit)
|
輸出為
(0.9793398380279541, 0.537902295589447)，第一個值為shapiro檢驗的統計量值，第二個值為對應概率p值,p值大於0.05置信水平，利潤變量服從正態分布假設。
K-S檢驗 KSValue = stats.kstest(rvs = datay,args = (datay.mean(),datay.std(),cdf = 'norm'))
需要注意的是KS檢驗，必須指定args參數，來傳遞檢驗變量的均值和標准差。
2）多重共線性檢驗 （自變量之間的檢測）
檢測模型中自變量之間存在較高的線性相關關系，若存在會給模型帶來嚴重的后果。其檢驗可以通過方差膨脹因子VIF來鑒定。（VIF>10，存在多重共線性，>100則變量間存在嚴重的多重共線性）
消除多重共線性，可以通過不同自變量之間的組合觀察vif值，確定最佳組合。R語言中消除多重共線性，確定最優組合使用模型的step方法（逐步回歸）。
計算方式：
X1為第一個自變量。構造x1與其余自變量的線性回歸模型，得到回歸模型的判決系數R2，得到第一個自變量的膨脹因子，依次類推計算剩余變量的VIF。
Python 的statsmodels模塊可直接計算VIF值
|
1
2
3
4
5
|
from
statsmodels.stats.outliers_influence
import
variance_inflation_factor
X
=
sm.add_constant(Profit_New.ix[:,[
'RD_Spend'
,
'Marketing_Spend'
]])
vif
=
pd.DataFrame()
vif[
'features'
]
=
X.columns
vif[
"VIF Factor"
]
=
[variance_inflation_factor(X.values,i)
for
i
in
range
(X.shape[
1
])]
|
|
1
|
得到vif的值：<br><br>
|
計算所有的變量:
|
1
2
3
4
5
|
X
=
sm.add_constant(Profit_New.drop(
'Profit'
,axis
=
1
).ix[:])
vif
=
pd.DataFrame()
vif[
'features'
]
=
X.columns
vif[
"VIF Factor"
]
=
[variance_inflation_factor(X.values,i)
for
i
in
range
(X.shape[
1
])]
vif
|
|
1
|
vif值:
|
第二種模型得到的const的值大於10，說明構建模型的變量之間存在多重共線性，需要刪除變量重新選擇模型。第一種模式選擇2個自變量（'RD_Spend','Marketing_Spend'）是符合自變量之間無線性相關要求的。
3）線性相關性檢驗 （各變量之間線性相關性檢測）
- 圖形法 seaborn pariplot函數
|
1
2
|
sns.pairplot(Profit_New.ix[:,[
'RD_Spend'
,
'Administration'
,
'Marketing_Spend'
,
'Profit'
]])
#去除啞變量
plt.show()
|
散點圖矩陣圖示：
觀察散點圖分布得出，RD_Spend，Marketing_Spend兩個變量和因變量Profit之間存在線性關系，Administration與Profit之間散點圖呈水平趨勢，說明2者之間不相關。
- 數據框的corrwith方法(該方式優點是可以計算任意指定變量之間的相關系數)
判斷條件：
|
1
2
|
# Profit 和各變量之間的相關系數
Profit_New.drop([
'Profit'
,
'State_California'
,
'State_Florida'
],axis
=
1
).corrwith(Profit_New.Profit)
|
結果看出：自變量RD_Spend和Marketing_Spend與因變量之間存在較高的線性相關性，而其他變量與利潤之間線性相關性比較弱。可忽略其線性相關關系。
需要注意的是通過corrwith檢測變量之間不存在線性相關關系，結合seaborn.pariplot觀察，可能存在其他關系如2次方，3次方，倒數關系等等（在變量通過t檢驗被確認為對因變量影響是顯著的情況下，還需要對模型再次檢查）。
通過如下方式重新確定模型為： Profit = 51902.112471 + 0.785116RD_Spend + 0.019402Marketing_Spend
|
1
2
|
model3
=
sm.formula.ols(
'Profit~RD_Spend + Marketing_Spend '
,data
=
train).fit()
model3.params
|
4）異常值檢測 （帽子矩陣、DFFITS准則、學生化殘差、Cook距離）
確認方式：|學生化殘差| > 2，則對應的數據點為異常值。
對異常值得處理：異常值檢測 model.get_influence() ， 學生化殘差值 model.get_influence().resid_studentized_external
確認異常樣本的比例，比例<5%，考慮可以直接刪除異常值，>5%，可以生成啞變量，對於異常點，設置啞變量為1，否則為0。
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|
# 異常值檢測
outliers
=
model3.get_influence()
# 帽子矩陣
leverage
=
outliers.hat_matrix_diag
# dffits值
dffits
=
outliers.dffits[
0
]
# 學生化殘差
resid_stu
=
outliers.resid_studentized_external
# cook距離
cook
=
outliers.cooks_distance[
0
]
# 合並各種異常值檢驗的統計量值
contatl
=
pd.concat([pd.Series(leverage, name
=
'leverage'
),
pd.Series(dffits, name
=
'dffits'
),
pd.Series(resid_stu, name
=
'resid_stu'
),
pd.Series(cook, name
=
'cook'
)
],axis
=
1
)
train.index
=
range
(train.shape[
0
])
profit_outliers
=
pd.concat([train,contatl ],axis
=
1
)
profit_outliers.head()
|
使用非異常值的數據點進行建模：
|
1
2
3
4
5
6
7
|
# 計算異常值比例
outliers_ratio
=
np.
sum
(np.where((np.
abs
(profit_outliers.resid_stu)>
2
),
1
,
0
))
/
profit_outliers.shape[
0
]
outliers_ratio
# 0.02564
none_outliers
=
profit_outliers.ix[np.
abs
(profit_outliers.resid_stu)<
=
2
,]
model4
=
sm.formula.ols(
'Profit~RD_Spend+Marketing_Spend'
,data
=
none_outliers).fit()
model4.params
|
再次更新模型為：Profit = 51827.416821 + 0.797038RD_Spend + 0.01774Marketing_Spend
5）殘差獨立性檢驗
通過mdoel.summary函數觀察Durbin-Watson的值，在2左右說明殘差之間是不相關的，偏離2較遠，說明殘差之間是不獨立的。
model4.summary()，DW的值為2.065，說明模型殘差滿足獨立性假設。
6）殘差方差齊性檢驗 （圖形法和BP法）
目的： 殘差項的方差不隨自變量的變動而呈現某種趨勢，即殘差項的方差不受自變量變化而影響
圖形法：繪制殘差和自變量之間的散點圖
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
ax1
=
plt.subplot2grid(shape
=
(
2
,
1
) , loc
=
(
0
,
0
))
# 設置第一張子圖位置
# 散點圖繪制
# ax1.scatter(none_outliers.RD_Spend, none_outliers.resid_stu ) #學生化殘差與自變量散點圖
ax1.scatter(none_outliers.RD_Spend, (model4.resid
-
model4.resid.mean())
/
model4.resid.std())
#標准化殘差和自變量散點圖
# 添加水平參考線
ax1.hlines(y
=
0
,
xmin
=
none_outliers.RD_Spend.
min
(),
xmax
=
none_outliers.RD_Spend.
max
(),
color
=
'red'
,
linestyle
=
'--'
)
ax1.set_xlabel(
'RD_Spend'
)
ax1.set_ylabel(
'Std_Residual'
)
ax2
=
plt.subplot2grid(shape
=
(
2
,
1
) , loc
=
(
1
,
0
))
# ax1.scatter(none_outliers.Marketing_Spend, none_outliers.resid_stu )
ax2.scatter(none_outliers.Marketing_Spend, (model4.resid
-
model4.resid.mean())
/
model4.resid.std())
ax2.hlines(y
=
0
,
xmin
=
none_outliers.Marketing_Spend.
min
(),
xmax
=
none_outliers.Marketing_Spend.
max
(),
color
=
'magenta'
,
linestyle
=
'--'
)
ax2.set_xlabel(
'Marketing_Spend'
)
ax2.set_ylabel(
'Std_Residual'
)
# 調整2子圖之間距離
plt.subplots_adjust(hspace
=
0.6
,wspace
=
0.3
)
plt.show()
|
注意的是采用學生化殘差或標准化殘差繪制和自變量的散點圖均可以，2者圖形類似。
原因是學生化殘差計算公式為：
圖形結果如下：
根據結果可知標准化殘差並沒有隨着自變量變化而呈現某種趨勢，圖形中所有點均勻的分布在y=0兩側，可以認為殘差項滿足方差齊性假設。
方差齊性檢測完整代碼：
BP法（拉格朗日乘子檢驗）：statsmodels模塊het_breuschpagan函數 https://programtalk.com/python-examples/statsmodels.stats/
|
1
|
sm.stats.diagnostic.het_breuschpagan(model4.resid, exog_het
=
model4.model.exog)
|
結果得到4個值
(1.4675103668308342, 0.48010272699006384, 0.7029751237162462, 0.5019659740962872)
第一個值為LM統計量，第二個為p值，第三個為F統計量，第四個為F統計量的p值，通過觀察p值均大於0.05，證明接受殘差方差齊性原假設，即殘差不受自變量的影響而變化。
若果殘差不滿足方差齊性假設，觀察殘差和自變量的關系，通過模型變換或加權最小二乘法對模型加以處理。
四、回歸模型預測
|
1
2
3
4
5
6
7
8
9
10
|
pred4
=
model4.predict(exog
=
test.ix[:,[
'RD_Spend'
,
'Marketing_Spend'
]])
plt.scatter(x
=
test.Profit ,y
=
pred4)
plt.plot([ test.Profit.
min
(),test.Profit.
max
()],
[test.Profit.
min
(),test.Profit.
max
()],
color
=
'red'
,
linestyle
=
'--'
)
plt.xlabel(
'實際值'
)
plt.ylabel(
'預測值'
)
plt.show()
|
真實值和預測值之間差異：
|
1
|
print
(pd.DataFrame({
"prediction"
:pred4,
'real'
:test.Profit}))
|
圖形展示效果：
