獨家 | 手把手教你用R語言做回歸後的殘差分析(附代碼)

作者：Abhijit Telang

翻譯：張睿毅

校對：丁楠雅

本文約2600字，建議閱讀10分鐘。

本文介紹了做殘差分析的方法及其重要性，以及利用R語言實現殘差分析。

在這篇文章中，我們通過探索殘差分析和用R可視化結果，深入研究了R語言。

殘差本質上是當一個給定的模型（在文中是線性回歸）不完全符合給定的觀測值時留下的gap。

醫學中的病理學發現的殘留分析是一個形象的比喻。人們通常用代謝殘留水平來作為衡量藥物吸收的指標。

殘差是用於建模的原始值與作為模型結果的對於原始值的估計之間的差異。

殘差=y-y-hat，其中y是初始值，y-hat是計算值。

期望這個錯誤儘可能接近於零，並且通過殘差找到任何異常值。

找到異常值的一個快速方法是使用標準化殘差。第一種方法是簡單地求出殘差與其標準差的比值，因此，任何超過3個標準差的情況都可以被視為異常值。

## 標準化殘差-相對於其標準偏差的比例殘差

residueStandard<-rstandard(lmfit)

df[residueStandard>3,]

以下是得到的結果：

days.instant days.atemp days.hum days.windspeed days.casual

442          442   0.505046 0.755833       0.110704        3155

456          456   0.421708 0.738333       0.250617        2301

463          463   0.426129 0.254167       0.274871        3252

470          470   0.487996 0.502917       0.190917        2795

471          471   0.573875 0.507917       0.225129        2846

505          505   0.566908 0.456250       0.083975        3410

512          512   0.642696 0.732500       0.198992        2855

513          513   0.641425 0.697083       0.215171        3283

533          533   0.594708 0.504167       0.166667        2963

624          624   0.585867 0.501667       0.247521        3160

645          645   0.538521 0.664167       0.268025        3031

659          659   0.472842 0.572917       0.117537        2806

當然，我希望我的模型是無偏的，至少我想這樣。因此回歸線兩邊的任何殘差，如果沒有在這條線上，都是隨機的，也就是說，沒有任何特定的模式。

也就是說，我希望我的剩餘誤差分布遵循一個普通的正態分布。

使用R語言，只需兩行代碼就可以優雅地完成這項工作。

hist(lmfit$residuals)

qqnorm(lmfit$residuals);qqline(lmfit$residuals)

於是，我們知道這個圖偏離了正常值（正常值用直線表示）。

但這種非黑即白的信息一般是不夠的。因此，我們應該檢查偏態和峰度，以了解分布的分散性。

首先，我們將計算偏態；我們將使用一個簡單的高爾頓偏態（Galton’s skewness）公式。

## 分布對稱性檢驗：偏態

summary<-summary(lmfit$residuals)

Q1<-summary[[2]]

Q2<-summary[[3]]

Q3<-summary[[5]]

skewness<-((Q1-Q2)+(Q3-Q2))/(Q3-Q1)

skewness

[1] 0.2665384

這量化了分布的方向（分布於右側或左側，或者是否完全對稱或均勻）。

另一個度量角度是峰度，它顯示分布是朝向中心（+ve值）還是遠離中心（-ve值）。峰度是一個量化離群值可能性的指標。正值表示離群的存在。

kurtosis(lmfit$residuals)

[1] 2.454145

接下來，除了可視化殘差的分布，我還想找出殘差是否相互關聯。

對於一個模型來說，為了解釋觀測值的所有變化，殘差必須隨機發生，並且彼此不相關。

Durbin-Watson測試允許檢驗殘差彼此之間的獨立性。

dwt(lm(df$days.windspeed~df$days.atemp+df$days.hum),method='resample', alternative='two.sided')

 lag Autocorrelation D-W Statistic p-value

   1       0.2616792      1.475681       0

 Alternative hypothesis: rho != 0

對於給定的自由度和觀測次數，需要將統計值與臨界值表確定的下限和上限進行比較。文中案例的值域是[1.55,1.67]。

由於計算的D-W統計值低於該範圍的較低值，我們拒絕了殘差不相關的零假設。因此取而代之的假設是殘差之間可能存在相關性。

直觀地看，這個假設可以通過研究模型在試圖捕獲原始Y值的增加值時失敗的原因來了解。當捕獲增加值時，隨著y的增加，殘差與y成正比。

將其與繪製擬合y-hat值與y值進行比較。當y-hat值趨於落後時，殘差似乎與y共同增長，故此，因為過去的殘值似乎繼續沿著固定的坡度值運行，過去的殘值似乎是當前值的更好預測因子。同時，在達爾文-沃森檢驗（Darwin-Watson tests）中在殘差與先前值之間的差的平方和，與所有觀測的給定殘差之和的比較和對比中，發現了相關性。

因此，殘餘誤差之間的相關性就像心率測試失靈一樣。如果你的模型不能跟上Y的快速變化，它會越來越和一個檢驗玩命的跑步者的情況類似，沒記錄上的步數似乎比跑步者實際的步數更相關。

重構你的線性回歸模型中的下意識影響。

有一點很重要：當對任何給定的觀測集進行線性回歸時，因變量（符號上表示為y）的計算估計量（符號上表示為y-hat）的每個值不僅依賴於當前值（例如，觀察值），還依賴於每次觀測。

畢竟，你的多元線性回歸模型只在最小化因變量y的實際值和所有觀測值y（y-hat）的計算估計值之間的誤差後，才計算出每個影響因素的係數。

在數學上，這可以用簡單的加權平均模型表示，如下所示。

Yhat(j)= w1j*Y(1)+w2j*Y(2)+w3j*Y(3)+w4j*Y(4)++wnj*Y(n);

where n is number of observations and j is one particular estimated value.

如果你想用矩陣表示法來考慮這個問題，矩陣w的維數是[n，n]，其中n是觀察到的數目。同樣，以矩陣形式表示的y具有尺寸[n，1]，因此y-hat[j]表示為矩陣，因為它是兩個矩陣（w和y）的乘積，並且具有與y相同的尺寸[n，1]。

因此，計算出的估計量的每一個值都可以表示為所有初始y值的加權平均值，表明給定的計算估計量在多大程度上受到每個初始y值的影響。

更簡單地說，如果我看到一堆分散的觀察結果，我必須以儘可能小的偏差，畫一條直線穿過它們，也就是說，如果我必須遺漏一些觀察結果，我必須在這條回歸線的兩側平等地遺漏。

但是在這樣做的同時，我（或選擇的回歸程序）不管來源何處，下意識地估計了每次觀察對我回歸線的斜率的影響程度。

我現在所做的就是重新計算這些影響。為了改變我的回歸線，將它們包括在內或排除在外，每一個觀察值的權重是多少？

權重的統計術語是hat values，因為它們連接了計算的估計量和它們的原始對應的估計量。

以下是用R語言計算的方法：

modelmatrix<-model.matrix(lmfit)

hatvalues<-hat(modelmatrix)

首先，我們得到一個矩陣形式的模型。然後我們計算hat值。

或者，可以使用以下函數獲得類似的結果。

hatvalues<-lm.influence(lmfit)$hat

讓我們考慮一下可以施加在每個權重上的限制。

顯然，權重的最小可能值等於所有原始Y值貢獻相等的可能性（因為它們必須為線性回歸程序貢獻一些東西，通過對所有觀測進行工作和優化來估計係數）。

在這種情況下，其值域的下限為1/n，其中n是觀測總數。在任何情況下，因為n總是+ve，所以權重總是<1。

現在試著將hat值加和，你會看到一個有趣的結果…

sum(hatvalues)

[1] 4

它們等於線性回歸模型為計算考慮的影響因素數量+1。

例如，在示例數據集中，我們有三個因素，即溫度、溼度和風速。

接下來，我們如何找到最重要或最有影響的觀察結果？

一種優雅的方式是：

將hat值切分為四分位數。

應用95%標準過濾最異常值。

將該過濾標準應用於觀察結果。

R語言允許你一步完成這些操作！

df[hatvalues>quantile(hatvalues,0.95),]

我們可以得到滿足這一標準的具體觀測結果。

 days.instant days.atemp days.hum days.windspeed days.casual

9              9  0.1161750 0.434167       0.361950          54

21            21  0.1578330 0.457083       0.353242          75

22            22  0.0790696 0.400000       0.171970          93

26            26  0.2036000 0.862500       0.293850          34

32            32  0.2345300 0.829565       0.053213          47

36            36  0.2430580 0.929167       0.161079         100

45            45  0.3983500 0.375833       0.417908         208

50            50  0.3914040 0.187917       0.507463         532

65            65  0.3662520 0.948261       0.343287         114

69            69  0.3856680 0.000000       0.261877          46

90            90  0.2575750 0.918333       0.217646         179

94            94  0.5429290 0.426250       0.385571         734

95            95  0.3983500 0.642083       0.388067         167

106          106  0.4254920 0.888333       0.340808         121

153          153  0.6439420 0.305000       0.292287         736

239          239  0.6355560 0.850000       0.375617         226

249          249  0.5152000 0.886957       0.343943         204

293          293  0.4665250 0.636250       0.422275         471

302          302  0.2279130 0.882500       0.351371          57

341          341  0.4002460 0.970417       0.266175          50

368          368  0.1262750 0.441250       0.365671          89

383          383  0.2752540 0.443333       0.415429         109

388          388  0.2430580 0.911250       0.110708         145

408          408  0.1016580 0.464583       0.409212          73

421          421  0.2556750 0.395833       0.421642         317

433          433  0.5246040 0.567500       0.441563         486

434          434  0.3970830 0.407083       0.414800         447

463          463  0.4261290 0.254167       0.274871        3252

465          465  0.4766380 0.317500       0.358196         905

478          478  0.3895040 0.835417       0.344546         120

543          543  0.5947040 0.373333       0.347642        1077

627          627  0.5650670 0.872500       0.357587         371

667          667  0.4677710 0.694583       0.398008         998

668          668  0.4394000 0.880000       0.358200           2

694          694  0.2487420 0.404583       0.376871         532

722          722  0.2361130 0.441250       0.407346         205

726          726  0.2203330 0.823333       0.316546           9

最終，你可以用顏色編碼使他們具象化。

##具象化

plot(df$days.casual,lmfit$fitted.values,col=ifelse(hatvalues>quantile(hatvalues,0.95),'red','blue'))

就是這樣！這是進行殘差分析和其重要性的簡單概述。

原文標題：

Doing Residual Analysis Post Regression in R

原文連結：

https://dzone.com/articles/doing-residual-analysis-post-regression-in-r

編輯：王菁

校對：龔力

張睿毅，北京郵電大學大二物聯網在讀。我是一個愛自由的人。在郵電大學讀第一年書我就四處跑去蹭課，折騰整一年驚覺，與其在當下焦慮，不如在前輩中沉澱。於是在大二以來，堅持讀書，不敢稍歇。資本主義國家的科學觀不斷刷新我的認知框架，同時因為出國考試很早出分，也更早地感受到自己才是那個一直被束縛著的人。太多真英雄在社會上各自閃耀著光芒。這才開始，立志終身向遇到的每一個人學習。做一個純粹的計算機科學裡面的小學生。

工作內容：需要一顆細緻的心，將選取好的外文文章翻譯成流暢的中文。如果你是數據科學/統計學/計算機類的留學生，或在海外從事相關工作，或對自己外語水平有信心的朋友歡迎加入翻譯小組。

你能得到：定期的翻譯培訓提高志願者的翻譯水平，提高對於數據科學前沿的認知，海外的朋友可以和國內技術應用發展保持聯繫，THU數據派產學研的背景為志願者帶來好的發展機遇。

其他福利：來自於名企的數據科學工作者，北大清華以及海外等名校學生他們都將成為你在翻譯小組的夥伴。

點擊文末「閱讀原文」加入數據派團隊~

轉載須知

如需轉載，請在開篇顯著位置註明作者和出處（轉自：數據派ID：datapi），並在文章結尾放置數據派醒目二維碼。有原創標識文章，請發送【文章名稱-待授權公眾號名稱及ID】至聯繫郵箱，申請白名單授權並按要求編輯。

發布後請將連結反饋至聯繫郵箱（見下方）。未經許可的轉載以及改編者，我們將依法追究其法律責任。

點擊「閱讀原文」擁抱組織