R語言從入門到精通:Day12--R語言統計--回歸分析

　　關注

　　從許多方面來看，回歸分析都是統計學的核心。它其實是一個廣義的概念，通指那些用一個或多個預測變量(也稱自變量)來預測響應變量(也稱因變量) 的方法。通常，回歸分析可以用來挑選與響應變量相關的預測變量，可以描述兩者的關係，也可以生成一個等式，通過預測變量來預測響應變量。

　　回歸分析在現代統計學中非常重要，本次教程內容安排如下：

　　首先：看一看如何擬合和解釋回歸模型，然後回顧一系列鑑別模型潛在問題的方法，並學習如何解決它們；

　　其次：我們將探究變量選擇問題（對於所有可用的預測變量，如何確定哪些變量包含在最終的模型中?）；

　　然後：我們將討論一般性問題（模型在現實世界中的表現到底如何?）；

　　最後：我們再看看相對重要性問題（模型所有的預測變量中，哪個最重要，哪個第二重要，哪個最無關緊要?）。

　　所以這次的內容很多，估計會成為這個系列中最長最複雜的一篇教程，不過堅持下來應該也是收穫最多的一篇教程。

　　

　　溫馨提示

　　1、本節內容重點內容較多，

　　務必緊跟紅色標記。

　　2、測試數據及代碼

　　見文末客服小姐姐二維碼。

　　回歸作為一個廣義的概念，涵蓋了許多變種，R語言中也為其提供了強大而豐富的函數和選項（但顯然選項越多，對初學者越不友好），早在2005年，R中就有200多種關於回歸分析的函數 (https://cran.r-project.org/doc/contrib/Ricci-refcard-regression.pdf，這個文檔提供了部分回歸分析函數列表，供大家參考)。

　　這些R函數對應了回歸分析的各種變體（如Logistic回歸，泊松回歸等等），而這次的內容主要關於OLS（普通最小二乘）回歸法，包括了簡單線性回歸、多項式回歸和多元線性回歸，下次再介紹其它常用的回歸分析。

　　OLS回歸是通過預測變量的加權和來預測量化的因變量，其中權重是通過數據估計而得的參數，一些典型的可以用OLS回歸解決的問題包括了：

　　

　　如果你讀到這裡發現自己還不知道什麼是OLS（普通最小二乘）回歸，建議去補習一下統計學知識了，否則作為未來的統計學家卻不知道基礎的統計學知識，也太不像話了，被說認識我科研貓~~

　　

　　
所以這裡默認大家都了解OLS回歸的數學原理和基礎術語，直接進入R函數的介紹了。

　　1、線性擬合的常用函數

　　在R中，擬合線性模型最基本的函數就是函數lm()，格式為:

　　myfit <- lm(formula, data)

　　回歸分析裡的參數formula對應著要擬合的模型形式，data是一個數據框，包含了用於擬合模型的數據。閱讀過本系列的上一篇關於方差分析的教程的同學對參數formula肯定很熟悉了，這裡同樣提供了一個參數formula中常用符號的表格：

　　

　　表1:參數formula中的常用符號

　　除了函數lm()，表2還列出了其他一些對做簡單或多元回歸分析有用的函數。擬合模型後，將這些函數應用於函數lm()返回的對象，可以得到更多額外的模型信息。

　　

　　表2: 對擬合線性模型非常有用的其他函數

　　2、回歸模型中的變量

　　當回歸模型包含一個因變量和一個自變量時，我們稱為簡單線性回歸。當只有一個預測變量， 但同時包含變量的冪(比如，X、X2、X3)時，我們稱為多項式回歸。當有不止一個預測變量時，則稱為多元線性回歸。首先從一個簡單的線性回歸例子開始，然後逐步展示多項式回歸和多元線性回歸，最後介紹一個包含交互項的多元線性回歸的例子。

　　基礎安裝中的數據集women提供了15個年齡在30~39歲間女性的身高和體重信息，我們通過身高來預測體重，獲得一個等式可以幫助我們分辨出那些過重或過輕的個體。圖1展示了擬合結果，通過代碼的輸出結果，可以得到預測等式：weight` = -87.52+3.45*height。輸出結果中的F 統計量檢驗所有的預測變量預測響應變量是否都在某個機率水平之上。由於簡單回歸只有一個預測變量，此處F檢驗等同於身高回歸係數的t檢驗。

　　

　　圖1:簡單線性回歸圖示

　　圖1 中可以看出我們其實可以用一個彎曲的曲線（多項式回歸）來提高預測的精度，多項式回歸允許你用一個解釋變量預測一個響應變量，它們關係的形式即n次多項式。圖2展示了擬合含二次項等式的結果，可以看出曲線確實擬合得較好。大家可以嘗試自己從代碼的輸出結果中得出預測的結果等式。

　　

　　圖2:多項式回歸

　　當然你也可以用三次甚至更高次的多項式來完成這次回歸分析，但我發現使用比三次更高的項幾乎沒有必要。car包中的函數scatterplot()可以很容易、方便地繪製二元關係圖，大家可以參考後臺代碼學習。

　　當預測變量不止一個時，簡單線性回歸就變成了多元線性回歸，分析也稍微複雜些。從技術上來說，多項式回歸可以算是多元線性回歸的特例:二次回歸有兩個預測變量(X和X2)，三次回歸有三個預測變量(X、X2和X3)。以基礎包中的state.x77數據集為例，探究一個州的犯罪率和其他因素的關係，包括人口、文盲率、平均收入和結霜天數(溫度在冰點以下的平均天數)。多元回歸分析中，第一步檢查一下變量間的相關性（如圖3）。

　　

　　圖3:變量相關性檢查

　　從圖中可以看到，謀殺率是雙峰的曲線，每個預測變量都一定程度上出現了偏斜。謀殺率隨著人口和文盲率的增加而增加，隨著收入水平和結霜天數增加而下降。同時，越冷的州府文盲率越低，收入水平越高。多元回歸的結果顯示文盲率的回歸係數為4.14，表示控制人口、收入和溫度不變時，文盲率上升1%，謀殺率將會上升4.14%，它的係數在p<0.001的水平下顯著不為0。相反，Frost的係數沒有顯著不為0(p=0.954)。

　　以上分析中，並沒有考慮預測變量的交互項。接下來，我們將考慮一個包含此因素的例子。以mtcars數據框中的汽車數據為例，把汽車重量和馬力作為預測變量，並包含交互項來擬合回歸模型。通過effects包中的函數effect()，可以用圖形展示交互項的結果。

　　

　　圖4:交互項圖形

　　回歸分析的結果告訴我們，馬力與車重的交互項是顯著的（p=0.00081），說明每加侖汽油行駛英裡數與汽車馬力的關係依車重不同而不同。若兩個預測變量的交互項顯著，說明響應變量與其中一個預測變量的關係依賴於另外一個預測變量的水平。從圖4中可以很清晰地看出，隨著車重的增加，馬力與每加侖汽油行駛英裡數的關係減弱了。當wt=4.2時，直線幾乎是水平的，表明隨著hp的增加，mpg不會發生改變。

　　3、模型的評估

　　討論完以上內容中，我們使用lm()函數來擬合OLS回歸模型，通過summary()函數獲取模型參數和相關統計量。但是，沒有任何輸出告訴我們模型是否合適，對模型參數推斷的信心依賴於它在多大程度上滿足OLS模型統計假設（這將決定回歸分析得出的模型應用到真實世界中時的預測效果）。下面我們要對模型進行診斷（回歸診斷）。

　　R基礎安裝中提供了大量檢驗回歸分析中統計假設的方法。最常見的方法就是對函數lm()返回的對象使用函數 plot()，可以生成評價模型擬合情況的四幅圖形。

　　

　　圖5:簡單回歸分析的診斷圖

　　理解上面這些圖形需要一些回歸分析的基礎知識，這可能需要你的數學老師花一個上午來講解，在這裡我只能簡單解釋四幅圖的含義：

　　

　　圖5中左上圖中可以清楚地看到一個曲線關係，這暗示著你可能需要對回歸模型加上一個二次項（這裡檢查了回歸分析統計假設中的「線性」）；

　　

　　右上圖檢查正態性，若滿足正態假設，那麼圖上的點應該落在呈45度角的直線上；

　　

　　左下圖檢查同方差性，滿足的條件下水平線周圍的點應該隨機分布；

　　

　　右下圖提供了你可能關注的單個觀測點的信息，從圖形可以鑑別出離群點、高槓桿值點和強影響點。

　　

　　圖6則是二次擬合的診斷圖。圖6表明多項式回歸擬合效果比較理想，基本符合了線性假設、殘差正態性(除了觀測點13)和同方差性(殘差方差不變)。觀測點15看起來像是強影響點(根據是它有較大的 Cook距離值)，刪除它將會影響參數的估計。事實上，刪除觀測點13和15，模型會擬合得會更好。但是對於刪除數據，要非常小心，因為本應是模型去匹配數據，而不是反過來。

　　

　　圖6:二次擬合的診斷

　　最後，用這個方法去診斷多元回歸分析的結果。

　　

　　圖7:多元回歸的診斷

　　這些R中的基礎函數的診斷結果對初學者並不友好，相信你們已經體會到了這一點，不過我們還有更好的工具可以選擇。car包中就提供很多診斷函數，見表3。

　　

　　表3:car包中的回歸診斷函數

　　利用這些函數依次檢查回歸分析的統計假設，函數qqPlot()提供了更為精確的正態假設檢驗方法；函數durbinWatsonTest()檢查誤差的獨立性；函數crPlots()檢查因變量與自變量之間是否呈非線性關係；函數ncvTest()，函數spreadLevelPlot()檢查同方差性；函數vif()檢查多重共線性。

　　

　　圖8:函數qqPlot()的結果

　　除了Nevada，所有的點都離直線很近，並都落在置信區間內，這表明正態性假設符合得很好。代碼中還提供了一個自定義的生成學生化殘差柱狀圖(即直方圖)，感興趣的同學可以試一下。函數durbinWatsonTest()的結果不顯著（p= 0.282)說明無自相關性，誤差項之間獨立。

　　

　　圖9:函數crPlots()的結果

　　圖9說明成分殘差圖證實了你的線性假設，線性模型形式對該數據集看似是合適的（如果不合適，就需要添加一些曲線成分，比如多項式項，或對一個或多個變量進行變換(如用log(X)代 替X)，或用其他回歸變體形式而不是線性回歸）。

　　

　　圖10:同方差性檢查結果

　　可以看到，函數ncvTest()的結果不顯著(p=0.19)，說明滿足方差不變假設，還可以通過分布水平圖 (圖10)看到這一點。

　　函數vif()的結果則表明預測變量不存在多重共線性問題。

　　最後，gvlma包中的函數gvlma()能對線性模型假設進行綜合驗證，同時還能做偏斜度、峰度和異方差性的評價。換句話說，它給模型假設提供了一個單獨的綜合檢驗(通過/不通過)。代碼中已提供示例。

　　4、異常值的處理

　　前面的回歸分析中出現了一些不符合模型的點，當時的建議是刪除這些「不聽話「的點，但這並不是一個嚴謹的辦法。一個全面的回歸分析要覆蓋對異常值的分析，包括離群點、高槓桿值點和強影響點。這些數據點需要更深入的研究，因為它們在一定程度上與其他觀測點不同，可能對結果產生較大的負面影響。

　　離群點是指那些模型預測效果不佳的觀測點。Q-Q圖不失為一種識別離群點的方法，car包也提供了一種離群點的統計檢驗方法，函數outlierTest()。

　　高槓桿值觀測點，即與其他預測變量有關的離群點。代碼中提供了一個自定義的函數來檢查這些點，結果如圖11。

　　

　　圖11:高槓桿值的檢查

　　強影響點，即對模型參數估計值影響有些比例失衡的點。有兩種方法可以檢測強影響點:Cook距離，或稱D統計量，以及變量添加圖。代碼中繪製了一個Cook距離的示例圖，圖12。圖中可以看到三個強影響點。

　　

　　圖12:Cook距離圖

　　Car包中函數avPlots()也有對應的功能，代碼中已提供例子。

　　當然，利用car包中的函數influencePlot()，你還可以將離群點、槓桿值和強影響點的信息整合到一幅圖形中(代碼中已提供例子)。

　　發現了這些異常點之後，一般有四種辦法來處理：刪除、變量變換、變量增刪、使用其他回歸方法。這四種方法中的變量變換在car包中有函數boxTidwell()和函數powerTransform() 幫助我們確定確定該如何進行變換（代碼已提供例子）。

　　變量增刪換一個說法，其實就是選擇合適的自變量，有以下兩種流行的方法:逐步回歸法(stepwise method)和全子集回歸(all-subsets regression)。MASS包中的函數stepAIC()可以實現前者（代碼例子中給出了向後回歸（逐步回歸中的一種））。全子集回歸則可用leaps包中的regsubsets()函數實現（代碼例子中提供兩個繪圖方法，供大家自己選擇）。變量的選擇在有大量變量的情況似乎會帶來很大的困擾，這時就需要背景知識來幫助你做出選擇，不要把時間浪費在毫無實際意義的變量上。

　　而第四種方法則帶來了新的問題：怎麼判斷哪種回歸模型是最適合數據的呢？當只需要在兩個模型之間選擇時，函數anova()和函數AIC()可以解決這個問題（代碼中已提供例子）。如果有100個甚至更多模型呢，交叉驗證就不失為一個好方法了。所謂交叉驗證，即將一定比例的數據挑選出來作為訓練樣本，另外的樣本作保留樣本，先在 訓練樣本上獲取回歸方程，然後在保留樣本上做預測。bootstrap包中的函數crossval()可以實現交叉驗證，在此基礎上可以自定義一個函數來對模型的R平方統計量做了k重交叉驗證（函數及例子見代碼）。

　　做完回歸分析，對模型進行診斷，解決了模型，異常點的問題，最後一個問題是哪些變量對模型最重要呢？或者說根據相對重要性對預測變量進行排序。最簡單的莫過於比較標準化的回歸係數，它表示當其他預測變量不變時，該預測變量一個標準差的變化可引起的響應變量的預期變化（在此之前，需要用函數scale()對數據進行標準化處理，例子見代碼）。同時，relaimpo包涵蓋了一些相對重要性的評價方法。最後，相對權重(relative weight)是一種比較有前景的新方法，它是對所有可能子模型添加一個預測變量引起的R平方平均增加量的一個近似值（代碼中提供了一個生成相對權重的函數）。圖13表明Illiteracy有最大的相對重要性，這和回歸係數的結果是一致的。

　　

　　圖13：相對權重

　　小結

　　總的來說，在統計中，回歸分析是許多方法的一個總稱，它是一個交互性很強的方法，包括擬合模型、檢驗統計假設、修正數據和模型，以及為達到最終結果的再擬合等過程，所以本次教程內容很長。下一次將會討論更加複雜的回歸模型，但總體思路是類似的，學習起來也不會很複雜。

　　加油

　　·

　　R語言統計-回歸分析

　　資源獲取方法新規

　　2

　　0

　　2

　　1

　　新年快樂

　　科研繪圖神器—hiplot，是2020年7月19日openbiox聯合科研貓鄭重推出的全網首個開源繪圖平臺，目前提供基於R語言的70餘種基礎可視化和50餘種進階繪圖的功能，同時還部署了多個 openbiox社區項目（如bget下載文獻附錄、UCSCXenaShiny 等）。

　　截止目前，網站的總訪問量大約37萬餘次，日均訪問量五千餘次，註冊用戶上萬人。

　　https://hiplot.com.cn

　　點擊圖片進入Hiplot平臺介紹

　　所有功能免費

　　更多科研新鮮資訊、文獻精讀和生物信息技能

特別聲明：以上內容(如有圖片或視頻亦包括在內)為自媒體平臺「網易號」用戶上傳並發布，本平臺僅提供信息存儲服務。

Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.