從零開始學Python數據分析【22】--線性回歸診斷(第一部分)

作者：劉順祥

個人微信公眾號：每天進步一點點2015

前文傳送門：

從零開始學Python數據分析【1】--數據類型及結構

從零開始學Python數據分析【2】-- 數值計算及正則表達式

從零開始學Python數據分析【3】-- 控制流與自定義函數

從零開始學Python數據分析【4】-- numpy

從零開始學Python數據分析【5】-- pandas(序列部分)

從零開始學Python數據分析【6】-- pandas(數據框部分01)

從零開始學Python數據分析【7】-- pandas(數據框部分02)

從零開始學Python數據分析【8】-- pandas(數據框部分03)

從零開始學Python數據分析【9】-- pandas(數據框部分04)

從零開始學Python數據分析【10】-- matplotlib(條形圖)

從零開始學Python數據分析【11】-- matplotlib(餅圖)

從零開始學Python數據分析【12】-- matplotlib(箱線圖)

從零開始學Python數據分析【13】-- matplotlib(直方圖)

從零開始學Python數據分析【14】-- matplotlib(折線圖)

從零開始學Python數據分析【15】-- matplotlib(散點圖)

從零開始學Python數據分析【16】-- matplotlib(雷達圖)

從零開始學Python數據分析【17】-- matplotlib(面積圖)

從零開始學Python數據分析【18】-- matplotlib(熱力圖)

從零開始學Python數據分析【19】-- matplotlib(樹地圖)

從零開始學Python數據分析【20】--線性回歸（理論部分）

從零開始學Python數據分析【21】--線性回歸（實戰部分）

前言

      在上一期中，關於線性回歸模型的創建，我們對比了Python和R語言的具體代碼實現，受到了很多網友的關注。也有一些朋友問到，關於線性回歸模型的那些前提假設為什麼沒有作分享，這期和下期我們就來回答一下這個問題。由於線性回歸模型的偏回歸係數通過最小二乘法（OLS）實現的，關於最小二乘法的使用是有一些假設前提的，具體是：

      當然，除了上面提到的5點之外，我們還需要注意的是，線性回歸模型對異常值是非常敏感的，模型預測的結果非常容易受到異常值的影響，所以，我們還需要對原始數據進行異常點識別和處理。在本期中，我們先針對上面提到的前三個假設做判別，此外再分享一下關於異常點的相關內容。本次分享的數據集來自於UCI【高爐煤氣聯合循環發電(CCPP)數據集(數據來源：http://archive.ics.uci.edu/ml/datasets/combined+cycle+power+plant)】。

線性性檢驗

      線性回歸模型，顧名思義，首先要保證自變量與因變量之間存在線性關係。關於線性關係的判斷，我們可以通過圖形或Pearson相關係數來識別，具體Python代碼如下：

# 導入第三方包
import pandas as pd
import numpy as np
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm
import scipy.stats as stats
from sklearn.metrics import mean_squared_error
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

# 數據讀取
ccpp = pd.read_excel('CCPP.xlsx')ccpp.describe()
# AT:溫度
# V:壓力
# AP:相對溼度
# RH:排氣量
# PE:發電量

# 繪製各變量之間的散點圖
sns.pairplot(ccpp)plt.show()

從上面的散點圖來看，似乎AP(相對溼度)和RH(排氣量)與PE(發電量)之間並不存在明顯的線性關係，具體我們還需要看一下PE與其餘變量之間的Pearson相關係數值。

# 發電量與自變量之間的相關係數
ccpp.corrwith(ccpp.PE)

      

從返回的結果來看，PE(發電量)與AT(溫度)和V(壓力)之間的相關係數還是蠻高的，而PE與AP(相對溼度)和RH(排氣量)之間的相關係數就明顯小很多。一般情況下，當Pearson相關係數低於0.4，則表明變量之間存在弱相關關係；當Pearson相關係數在0.4~0.6之間，則說明變量之間存在中度相關關係；當相關係數在0.6以上時，則反映變量之間存在強相關關係。經過對比發現，PE與RH之間的為弱相關關係，故不考慮將該變量納入模型。當然，變量之間不存在線性關係並不代表不存在任何關係，可能是二次函數關係、對數關係等，所以一般還需要進行檢驗和變量轉換。

多重共線性檢驗

      先來了解一下，如果模型的自變量之間存在嚴重的多重共線性的話，會產生什麼後果呢?

導致OLS估計量可能無效；

增大OLS估計量的方差；

變量的顯著性檢驗將失去意義；

模型缺乏穩定性；

      所以，多重共線性檢驗就顯得非常重要了，關於多重共線性的檢驗可以使用方差膨脹因子(VIF)來鑑定，如果VIF大於10，則說明變量存在多重共線性。一旦發現變量之間存在多重共線性的話，可以考慮刪除變量和重新選擇模型（嶺回歸法）。

# 將因變量PE，自變量AT,V,AP和截距項（值為1的1維數組）以數據框的形式組合起來
y, X = dmatrices('PE~AT+V+AP', data = ccpp, return_type='dataframe')

# 構造空的數據框
vif = pd.DataFrame()vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]vif["features"] = X.columnsvif

結果顯示，所有自變量的VIF均低於10，說明自變量之間並不存在多重共線性的隱患。

異常點檢測

      在異常點檢測之前，我們需要對現有的數據，進行線性回歸模型的構造。具體操作如下，與上一期文章中回歸模型的建模一致：

# 構造PE與AT、V和AP之間的線性模型
fit = sm.formula.ols('PE~AT+V+AP',data = ccpp).fit()fit.summary()

# 計算模型的RMSE值
pred = fit.predict()np.sqrt(mean_squared_error(ccpp.PE, pred))

通過上面的建模結果來看，一切都顯得那麼的完美（模型的顯著性檢驗通過，偏回歸係數的顯著性檢驗通過；R方值也達到了0.9以上）。儘管這樣，我們還是需要基於這個模型，完成異常點的檢測。
      關於異常點的檢測方法，一般可以通過高槓桿值點（帽子矩陣）或DFFITS值、學生化殘差、cook距離和covratio值來判斷。這些值的具體計算腳本如下：

# 離群點檢驗
outliers = fit.get_influence()

# 高槓桿值點（帽子矩陣）
leverage = outliers.hat_matrix_diag
# dffits值
dffits = outliers.dffits[0]
# 學生化殘差
resid_stu = outliers.resid_studentized_external
# cook距離
cook = outliers.cooks_distance[0]
# covratio值
covratio = outliers.cov_ratio

# 將上面的幾種異常值檢驗統計量與原始數據集合併
contat1 = pd.concat([pd.Series(leverage, name = 'leverage'),pd.Series(dffits, name = 'dffits'),                     pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook, name = 'cook'),                     pd.Series(covratio, name = 'covratio'),],axis = 1)ccpp_outliers = pd.concat([ccpp,contat1], axis = 1)ccpp_outliers.head()

通過參考薛毅老師的《統計建模與R軟體》書可知，當高槓桿值點（或帽子矩陣）大於2(p+1)/n時，則認為該樣本點可能存在異常（其中p為自變量的個數，n為觀測的個數）；當DFFITS統計值大於2sqrt((p+1)/n)時 ，則認為該樣本點可能存在異常；當學生化殘差的絕對值大於2，則認為該樣本點可能存在異常；對於cook距離來說，則沒有明確的判斷標準，一般來說，值越大則為異常點的可能性就越高；對於covratio值來說，如果一個樣本的covratio值離數值1越遠，則認為該樣本越可能是異常值。這裡我們就以學生化殘差作為評判標準，因為其包含了帽子矩陣和DFFITS的信息。

# 計算異常值數量的比例
outliers_ratio = sum(np.where((np.abs(ccpp_outliers.resid_stu)>2),1,0))/ccpp_outliers.shape[0]outliers_ratio

# 刪除異常值
ccpp_outliers = ccpp_outliers.loc[np.abs(ccpp_outliers.resid_stu)<=2,]

結果顯示，確實存在異常值點，且異常值的數量佔了3.7%。對於異常值的處理，我們可以考慮下面幾種辦法：

      這裡為了簡單起見，我們將3.7%的異常值作刪除處理。

# 重新建模
fit2 = sm.formula.ols('PE~AT+V+AP',data = ccpp_outliers).fit()fit2.summary()

# 計算模型的RMSE值
pred2 = fit2.predict()np.sqrt(mean_squared_error(ccpp_outliers.PE, pred2))

通過對比fit和fit2，將異常值刪除後重新建模的話，效果會更理想一點，具體表現為：信息準則（AIC和BIC）均變小，同時RMSE（誤差均方根）也由原來的4.89降低到4.26。

正態性檢驗

      當模型的殘差服從正態性假設時，才能保證模型偏回歸係數對於的t值和模型的F值是有效的。故模型建好後，要對模型的殘差進行正態性檢驗。關於正態性檢驗由兩類方法，即定性的圖形法（直方圖、PP圖和QQ圖）和定量的非參數法（Shapiro檢驗和K-S檢驗）。

# 殘差的正態性檢驗（直方圖法）
resid = fit2.resid

# 中文和負號的正常顯示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = False

plt.hist(resid, # 繪圖數據        bins = 100, # 指定直方圖條的個數        normed = True, # 設置為頻率直方圖        color = 'steelblue', # 指定填充色        edgecolor = 'k') # 指定直方圖的邊界色

# 設置坐標軸標籤和標題
plt.title('殘差直方圖')plt.ylabel('密度值')

# 生成正態曲線的數據
x1 = np.linspace(resid.min(), resid.max(), 1000)normal = mlab.normpdf(x1, resid.mean(), resid.std())

# 繪製正態分布曲線
plt.plot(x1,normal,'r-', linewidth = 2, label = '正態分布曲線') # 生成核密度曲線的數據
kde = mlab.GaussianKDE(resid)x2 = np.linspace(resid.min(), resid.max(), 1000)

# 繪製核密度曲線
plt.plot(x2,kde(x2),'k-', linewidth = 2, label = '核密度曲線')

# 去除圖形頂部邊界和右邊界的刻度
plt.tick_params(top='off', right='off')

# 顯示圖例
plt.legend(loc='best')
# 顯示圖形
plt.show()

從殘差的直方圖來看，核密度曲線與理論的正態分布曲線的趨勢還是比較吻合的，即使殘差不服從正態分布，也能反映其基本近似於正態分布。

# 殘差的正態性檢驗（PP圖和QQ圖法）
pp_qq_plot = sm.ProbPlot(resid)pp_qq_plot.ppplot(line = '45')plt.title('P-P圖')pp_qq_plot.qqplot(line = 'q')plt.title('Q-Q圖')# 顯示圖形plt.show()

從PP圖和QQ圖來看，有一部分樣本點並沒有落在參考線上，但絕大多數樣本點還是與參考線保持一致的步調。

# 殘差的正態性檢驗（非參數法）
standard_resid = (resid-np.mean(resid))/np.std(resid)stats.kstest(standard_resid, 'norm')

由於shapiro正態性檢驗對樣本量的要求是5000以內；而本次數據集的樣本量由9000多，故選擇K-S來完成正態性檢驗。從K-S檢驗的P值來看，拒絕了殘差服從正態分布的假設，即認為殘差並不滿足正態性假設這個前提。如果殘差不服從正態分布的話，建議對Y變量進行box-cox變換處理。由於fit2模型的殘差並沒有特別明顯的偏態（偏度為0.058，接近於0），故這裡就不對Y變量進行box-cox變換了。如果需要變換的話，可以以下面的代碼為例：

import scipy.stats as stats# 找到box-cox變換的lambda係數lamd = stats.boxcox_normmax(your_data_frame.y, method = 'mle')# 對Y進行變換your_data_frame['trans_y'] = stats.boxcox(your_data_frame.y, lamd)# 建模fit = sm.formula.ols('y~x1+x2+...', data = your_data_frame).fit()fit.summary()

      關於線性回歸模型的異常點識別、線性性假設、無多重共線性假設和殘差的正態性假設在Python中的實現就介紹到這裡。下一期，我們將針對線性回歸模型殘差的方差齊性和獨立性假設進行講解，接下來我們再用R語言對上面的內容復現一遍。

R語言腳本復現

# 加載第三方包
library(readxl)
library(GGally)

# 讀取數據
ccpp <- read_excel(path = file.choose())summary(ccpp)

# 繪製各變量之間的散點圖與相關係數
ggpairs(ccpp)

# 建模
fit <- lm(PE ~ AT + V + AP, data = ccpp)summary(fit)

# 計算模型的RMSE值
RMSE = sqrt(mean(fit$residuals**2))RMSE

# 多重共線性檢驗
vif(fit)

# 異常點檢驗# 高槓桿值點（帽子矩陣）
leverage <- hatvalues(fit)head(leverage)

# dffits值
Dffits <- dffits(fit)head(Dffits)

# 學生化殘差
resid_stu <- Dffits/sqrt(leverage/(1-leverage))head(resid_stu)

# cook距離
cook <- cooks.distance(fit)head(cook)

# covratio值
Covratio <- covratio(fit)head(Covratio)

# 將上面的幾種異常值檢驗統計量與原始數據集合併
ccpp_outliers <- cbind(ccpp, data.frame(leverage, Dffits, resid_stu, cook, Covratio))head(ccpp_outliers)

# 計算異常值數量的比例
outliers_ratio = sum(abs(ccpp_outliers$resid_stu)>2)/nrow(ccpp_outliers)outliers_ratio

# 刪除異常值
ccpp_outliers = ccpp_outliers[abs(ccpp_outliers$resid_stu)<=2,]

# 重新建模
fit2 = lm(PE~AT+V+AP,data = ccpp_outliers)summary(fit2)

# 計算模型的RMSE值
RMSE2 = sqrt(mean(fit2$residuals**2))RMSE2

# 正態性檢驗#繪製直方圖
hist(x = fit2$residuals, freq = FALSE,     breaks = 100, main = 'x的直方圖',     ylab = '核密度值',xlab = NULL, col = 'steelblue')

#添加核密度圖
lines(density(fit2$residuals), col = 'red', lty = 1, lwd = 2)

#添加正態分布圖
x <- fit2$residuals[order(fit2$residuals)]lines(x, dnorm(x, mean(x), sd(x)),      col = 'blue', lty = 2, lwd = 2.5)

#添加圖例
legend('topright',legend = c('核密度曲線','正態分布曲線'),       col = c('red','blue'), lty = c(1,2),       lwd = c(2,2.5), bty = 'n')

# PP圖
real_dist <- ppoints(fit2$residuals)theory_dist <- pnorm(fit2$residuals, mean = mean(fit2$residuals),                     sd = sd(fit2$residuals))

# 繪圖
plot(sort(theory_dist), real_dist, col = 'steelblue',     pch = 20, main = 'PP圖', xlab = '理論正態分布累計概率',     ylab = '實際累計概率')

# 添加對角線作為參考線
abline(a = 0,b = 1, col = 'red', lwd = 2)

# QQ圖
qqnorm(fit2$residuals, col = 'steelblue', pch = 20,       main = 'QQ圖', xlab = '理論分位數',       ylab = '實際分位數')
# 繪製參考線
qqline(fit2$residuals, col = 'red', lwd = 2)

# shapiro正態性檢驗

# shapiro <- shapiro.test(fit2$residuals)
# shapiro

# K-S正態性檢驗
ks <- ks.test(fit2$residuals, 'pnorm',              mean = mean(fit2$residuals),              sd = sd(fit2$residuals))ks

結語

      OK，今天關於線性回歸診斷部分的內容就分享到這裡，限於篇幅，不能一次性寫完，故在下一期將繼續推出診斷的其他內容（殘差的方差齊性和殘差的獨立性檢驗）。希望對數據挖掘或機器學習比較感興趣的朋友，能夠靜下心來好好的復現一遍。如果你有任何問題，歡迎在公眾號的留言區域表達你的疑問。同時，也歡迎各位朋友繼續轉發與分享文中的內容，讓更多的朋友學習和進步。

相關材料下載連結

連結: https://pan.baidu.com/s/1qYNsP0w 密碼: 2g3f

Python愛好者社區歷史文章大合集：

Python愛好者社區歷史文章列表（每周append更新一次）

福利：文末掃碼立刻關注公眾號，「Python愛好者社區」，開始學習Python課程：

關注後在公眾號內回復「課程」即可獲取：

0.小編的Python入門視頻課程！！！

1.崔老師爬蟲實戰案例免費學習視頻。

2.丘老師數據科學入門指導免費學習視頻。

3.陳老師數據分析報告製作免費學習視頻。

4.玩轉大數據分析！Spark2.X+Python 精華實戰課程免費學習視頻。

5.丘老師Python網絡爬蟲實戰免費學習視頻。