R與生物專題 | 第四十八講 R-逐步回歸

2021-02-19 投必得學術

在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。以從淺入深，層層遞進的形式在投必得醫學公眾號更新。



在我們學習了線性回歸和邏輯回歸等回歸模型對連續型變量或分類變量進行預測模型建模（可點擊上面「R與生物統計專題」進一步查找詳情）。同時，我們也在之前的講解中，多次提到，如果預測變量（自變量）與結果變量（因變量）間不相關，即該預測模型的係數b與0無統計學差異，我們在最終的預測模型中，需要將該預測變量刪除。


那麼，假如我們有很多個預測變量，想要一一確定這些預測變量是否真對結果變量起到了預測功能，我們有沒有辦法讓這個「刪除」變量的步驟可以自動實現呢？於是我們引入了逐步回歸。



逐步回歸的基本思想是將變量逐個引入模型，每引入一個預測變量（解釋變量）後都要進行F檢驗（第十九講 F檢驗：兩樣本方差比較），並對已經選入的預測變量逐個進行t檢驗（第十講 R-兩獨立樣本t檢驗），當原來引入的預測變量由於後面預測變量的引入變得不再顯著時，則將其刪除。以確保每次引入新的變量之前回歸方程中只包含顯著性變量。這是一個反覆的過程，直到既沒有顯著的預測變量選入回歸方程，也沒有不顯著的預測變量從回歸方程中剔除為止。以保證最後所得到的預測變量集是最優的。


依據上述思想，可利用逐步回歸篩選並剔除引起多重共線性的變量（第三十七講 R-多元線性回歸中的多重共線性和方差膨脹因子），其具體步驟如下：先用結果變量對每一個所考慮的預測變量做簡單回歸，然後以對結果變量貢獻最大的預測變量所對應的回歸方程為基礎，再逐步引入其餘預測變量。經過逐步回歸，使得最後保留在模型中的解釋變量既是重要的，又沒有嚴重多重共線性。


逐步回歸有三種實現策略，我們一般指的逐步回歸為第三種：


1.正向（Forward）選擇，從模型中沒有預測因素開始，反覆添加最有幫助的預測因素，直到沒有顯著的預測變量選入回歸方程 。

2.向後（Backward）選擇（也稱向後消除），從完整模型（即包含所有可能預測變量的模型）中的所有預測變量開始，以迭代方式刪除貢獻最小的預測變量，直到沒有不顯著的預測變量從回歸方程刪除。

3.逐步（stepwise）選擇（也稱順序替換），這是向前和向後選擇的組合。您從沒有預測變量開始，然後順序添加最有貢獻的預測變量（例如正向選擇）。添加每個新變量後，刪除所有不再改善模型擬合的變量（例如向後選擇），直到既沒有顯著的預測變量選入回歸方程，也沒有不顯著的預測變量從回歸方程中剔除為止。



·正向選擇和逐步選擇可以應用在高維數據情況下，即樣本數量n小於預測變量p的數量，例如在基因組領域，全基因組數量遠大於測試人群。


·向後選擇要求樣本數n大於變量p數，以便可以擬合整個模型。




·tidyverse 便於數據操作和可視化

·caret 簡化工作流程

·leaps 用於計算逐步回歸

library(tidyverse)library(caret)library(leaps)


我們將使用R的內置數據集swiss。

library(datasets)data("swiss")


Fertility Agriculture Examination Education Catholic Infant.Mortality1      77.3        89.7           5         2   100.00             18.32      76.1        35.3           9         7    90.57             26.63      83.1        45.1           6         9    84.84             22.2

研究問題：根據社會經濟相關的多個指標（Agriculture，Examination，Education，Catholic，Infant.Mortality）預測生育力得分（Fertility）。



有許多函數和R包可用於計算逐步回歸。


·stepAIC()[MASS包]，它根據AIC值選擇最佳模型。它具有一個選項direction，該選項可以選擇以下值：

「both」（用於逐步回歸 ）
「backward」（用於後向選擇）
「forward」（用於前向選擇）

它會返回最佳的最終模型。
備註：
AIC：（Akaike information criterion，簡稱AIC）是評估統計模型的複雜度和衡量統計模型「擬合」優度的指標之一。


 建立完全模型
full.model<-lm(Fertility~.,data=swiss)summary(full.model)

 輸出結果
Call:lm(formula = Fertility ~ ., data = swiss) Residuals:          Min       1Q   Median       3Q      Max  -15.2743  -5.2617   0.5032   4.1198  15.3213       Coefficients:                       Estimate Std. Error t value Pr(>|t|)  (Intercept)      66.91518   10.70604   6.250 1.91e-07 *** Agriculture      -0.17211    0.07030  -2.448  0.01873 *  Examination      -0.25801    0.25388  -1.016  0.31546    Education        -0.87094    0.18303  -4.758 2.43e-05 *** Catholic          0.10412    0.03526   2.953  0.00519 ** Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** --- Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1  Residual standard error: 7.165 on 41 degrees of freedom Multiple R-squared:  0.7067,   Adjusted R-squared:  0.671 F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

 逐步回歸
step.model <- stepAIC(full.model,direction = "both",                      trace =FALSE) summary(step.model)

Call:lm(formula = Fertility ~ Agriculture + Education + Catholic +    Infant.Mortality, data = swiss)  Residuals:              Min       1Q   Median       3Q      Max         -14.6765  -6.0522   0.7514   3.1664  16.1422 Coefficients:                  Estimate Std. Error t value Pr(>|t|)   (Intercept)      62.10131    9.60489   6.466 8.49e-08 ***Agriculture      -0.15462    0.06819  -2.267  0.02857 * Education        -0.98026    0.14814  -6.617 5.14e-08 ***Catholic          0.12467    0.02889   4.315 9.50e-05 ***Infant.Mortality  1.07844    0.38187   2.824  0.00722 **---Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1 Residual standard error: 7.168 on 42 degrees of freedomMultiple R-squared:  0.6993,   Adjusted R-squared:  0.6707F-statistic: 24.42 on 4 and 42 DF,  p-value: 1.717e-10


train()函數[caret包]提供了一個簡單的工作流程，可以使用leaps包和MASS包進行逐步選擇。它有一個名為的選項method，可以選擇以下值：

"leapBackward"，以使線性回歸與向後選擇相適應
"leapForward"，以使線性回歸與正向選擇相適應
"leapSeq"，以逐步選擇擬合線性回歸

您還需要指定調整參數nvmax，該參數對應於要納入模型的最大預測變量數。
例如，您可以nvmax在1到5之間變化。在這種情況下，該函數從搜索大小不同的最佳模型開始，一直到最佳有5個預測變量數的變量模型。也就是說，它搜索最佳的1變量模型，最佳2變量模型，…，最佳5變量模型。

以下示例method = "leapBackward"使用swiss數據集執行後向選擇，根據社會經濟指標確定預測生育力的最佳模型。

由於數據集僅包含5個預測變量，因此我們將nvmax在1到5 個變量之間變化，以識別出具有不同大小的5種最佳模型：最佳1變量模型，最佳2變量模型，…，最佳5變量模型。

我們將使用10倍交叉驗證來估計5個模型中每個模型的平均預測誤差（RMSE）（第四十二講 R-回歸預測模型的交叉驗證）。RMSE統計量用於比較5個模型並自動選擇最佳模型，其中最佳定義為RMSE最小的模型。

install.packages('leaps')library(leaps)library(MASS)

 設立隨機數種子

設置重複的k倍交叉驗證
train.control <- trainControl(method = "cv",number = 10)

 訓練模型
step.model <- train(Fertility ~.,data =swiss,                   method = "leapBackward",                   tuneGrid =data.frame(nvmax = 1:5),                                  trControl =train.control                                    )step.model$results
 輸出結果
nvmax RMSE Rsquared  MAE RMSESD RsquaredSD MAESD1     1 9.30    0.408 7.91   1.53      0.390  1.652     2 9.08    0.515 7.75   1.66      0.247  1.403     3 8.07    0.659 6.55   1.84      0.216  1.574     4 7.27    0.732 5.93   2.14      0.236  1.675     5 7.38    0.751 6.03   2.23      0.239  1.64

上面的輸出顯示了不同的指標及其標準誤，用於比較5種最佳模型的準確性。各列內容是：


在我們的示例中，可以看到具有4個變量（nvmax = 4）的模型是具有最低RMSE的模型。您可以顯示由train()函數自動選擇的最佳調整值（nvmax），如下所示：

這表明最好的模型是nvmax = 4個變量的模型。該函數summary()報告每種模型大小的最佳變量集，直到最佳的4變量模型。

summary(step.model$finalModel)
 Subset selection object 5 Variables  (and intercept)                  Forced in Forced out Agriculture          FALSE      FALSE Examination          FALSE      FALSE  Education            FALSE      FALSE  Catholic             FALSE      FALSE Infant.Mortality     FALSE      FALSE1 subsets of each size up to 4 Selection Algorithm: backward         Agriculture Examination Education Catholic Infant.Mortality 1  ( 1 ) " "         " "         "*"       " "      " "             2  ( 1 ) " "         " "         "*"       "*"      " "             3  ( 1 ) " "         " "         "*"       "*"      "*"             4  ( 1 ) "*"         " "         "*"       "*"      "*"


星號指定給定變量是否包含在相應的模型中。例如，可以看出最好的4變量模型包含Agriculture，Education，Catholic，和Infant.Mortality。Fertility ~ Agriculture + Education + Catholic + Infant.Mortality。

最終模型的回歸係數（id = 4）可以按以下方式訪問：
coef(step.model$finalModel, 4)
輸出結果 
(Intercept)      Agriculture        Education         Catholic Infant.Mortality      62.1013116       -0.1546175       -0.9802638        0.1246664        1.0784422

或者，僅使用選定的預測變量來計算線性模型：
lm(Fertility ~ Agriculture + Education + Catholic + Infant.Mortality,   data = swiss)
輸出結果
Call:lm(formula = Fertility ~ Agriculture + Education + Catholic +    Infant.Mortality, data = swiss)  Coefficients:     (Intercept)       Agriculture         Education          Catholic          62.101            -0.155            -0.980             0.125  Infant.Mortality            1.078

參考內容：
1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R
2. 百度百科：逐步回歸

在之後的更新中，我們會進一步為您介紹R的入門，以及常用生物統計方法和R實現。歡迎關注，投必得學術手把手帶您走入R和生物統計的世界。提前預告一下，下一講我們將學習 R-分類預測模型及診斷性試驗性能指標。
當然啦，R語言的掌握是在長期訓練中慢慢積累的。一個人學習太累，不妨加入「R與統計交流群」，和數百位碩博一起學習。
快掃二維碼撩客服，
帶你進入投必得學術交流群，
讓我們共同進步！
↓↓

來源：投必得醫學R與生物統計
聲明：本文僅做學術分享，版權歸原作者所有，並不代表本平臺的觀點，如有侵權，請先聯繫topedit2021刪除，萬分感謝！
發表SCI論文很迷茫？
來找「投必得」幫忙↓↓↓

相關焦點

  • R與生物專題 | 第五十四講 R-樣本量及實驗效能計算
    t檢驗（兩個n不相等的樣本）7.1 對t檢驗（兩組樣本相等）進行效能計算的函數pwr.t2n.test7.7 對平衡型方差分析進行樣本量／效能計算的函數pwr.anova.test()第十四講 R-單因素方差分析1第十五講 R-單因素方差分析2第十六講
  • R與生物專題 | 第三十二講 R-回歸分析概述
    回歸分析回歸分析（或回歸模型）是機器學習中監督性學習的方法之一，它是由一個或多個預測變量／自變量（x）的值預測連續性數據類型的結果變量／因變量（y）。線性回歸線性回歸（linear-regression）是預測連續變量的最簡單，最流行的技術。它的前提假設是：結果與預測變量之間存在線性關係。
  • R與生物專題 | 第三十五講 R-多元線性回歸
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。
  • R與生物專題 | 第四十四講 R-非線性回歸
    當自變量（預測變量）與因變量（結果變量）的真實關係為非線性時，之前學過的線性回歸模型將不再適用。這可以使用mgcv包來完成：model<-gam(medv~s(lstat), data=train.data)predictions<-model%>%predict(test.data)data.frame( RMSE=RMSE(predictions, test.data$medv), R2=R2(predictions
  • R與生物專題 | 第三十四講 R-簡單線性回歸模型(2)
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。
  • R與生物專題 | 第四十三講 R-回歸預測模型的自舉重採樣驗證(boostrap-resampling)
    ,data=swiss,method= "lm", trControl=train.control)#輸出結果Linear Regression47 samples5 predictorNo pre-processingResampling: Bootstrapped (100 reps)Summary of sample
  • 第四十二講 R-逐步回歸
    於是我們引入了逐步回歸。逐步回歸的基本思想是將變量逐個引入模型，每引入一個預測變量（解釋變量）後都要進行F檢驗，並對已經選入的預測變量逐個進行t檢驗，當原來引入的預測變量由於後面預測變量的引入變得不再顯著時，則將其刪除。以確保每次引入新的變量之前回歸方程中只包含顯著性變量。
  • R與生物專題 | 第五十二講 統計方法選擇全解析2
    （2）如果要做兩變量之間的相關性，可採用Pearson相關分析。（1）b＋c>40，則用McNemar配對檢驗。R函數：mcnemar.test()（1）配對比較：用McNemar配對檢驗。R函數：mcnemar.test()（2）一致性檢驗，用Kappa檢驗。
  • 第三十一講 R-機器學習與回歸概述
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。
  • R與生物專題 | 第六講 R-數據正態分布檢驗
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。
  • R與生物專題 | 第五十講 R-ROC曲線
    我們就需要用到ROC曲線來幫助大家找到最好的診斷切點，用以定義診斷試驗／預測模型的最佳狀態。 曲線下面積（Area under the Curve of ROC (AUCROC)）總結了預測模型的整體性能，在比較不同的分類模型時，可以將每個模型的ROC曲線都畫出來，比較曲線下面積做為模型優劣的指標。它表示預測模型區分事件與非事件的能力。
  • R與生物專題 | 第七講 R-相關性分析及作圖
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。以從淺入深，層層遞進的形式在投必得學術公眾號更新。在上一講中，我們介紹了第六講 R-數據正態分布檢驗（點我可跳轉）。這兩天大家學習積極性很高，在群裡紛紛交作業，群裡討論也很活躍，想入群的拉倒文末撩小編哦。
  • R語言實戰:回歸
    多層多元線性多變量Logistic泊松Cox比例風險時間序列非線性非參數穩健普通最小二乘 (OLS) 回歸法，包括簡單線性回歸，多項式回歸和多元線性回歸。 擬合回歸模型簡單線性回歸多項式回歸多元線性回歸簡單線性回歸數據集head(women) height weight1 58 1152 59 1173 60 1204 61 1235 62 126
  • R與生物專題 | 第十講 R-兩獨立樣本t檢驗
    單樣本的T檢驗和單樣本Wilcoxon檢驗用於檢驗一組樣本和理論值是否有差異。當我們比較兩組樣本間是否存在差異時，就需要用到今天的統計學方法：兩獨立樣本t檢驗。如果要檢驗男性的體重是否小於女性的體重，請輸入以下內容：t.test(weight ~ group, data = my_data,var.equal = TRUE, alternative = "less")
  • R與生物專題 | 第十四講 R-單因素方差分析1
    在「R與生物統計專題」中，我們會從介紹R的基本知識展開到生物統計原理及其在R中的實現。
  • R與生物專題 | 第十一講 R-兩獨立樣本Wilcoxon檢驗
    Woman 28.510 Man 37.811 Man 80.012 Man 33.413 Man 36.014 Man 89.415 Man 83.316 Man 97.317 Man 81.318 Man 92.4按性別計算統計信息（中位數和四分位數範圍
  • R語言和Python實現回歸分析
    變量之間關係可以分為兩類：函數關係：反映了事務之間某種確定性關係相關關係：兩個變量之間存在某種依存關係，但二者並不是一一對應的；反映了事務間不完全確定關係；相關係數(r)r的取值範圍是[-1,1]，r=1表示完全正相關！r=-1表示完全負相關！r=0表示完全不相關。為什麼要對相關係數進行顯著性檢驗？
  • 第四十講 R-線性回歸:預測模型及可信區間
    在「R與生物統計專題」中，我們會從介紹
  • R與數學專題 第三節: R求解線性代數
    結束語：通過這周的學習，大家對線性方程組的求解有了一定的認識，這周的課程只是帶領大家入門，例如超定方程組可以應用於多元線性回歸方程的求解。關於線性方程的求解方法還有很多，比如高斯消元法，迭代法，這些方法對於求解大矩陣的方程組更有效。
  • R邏輯回歸逐步選擇變量的實現
    本文是Python邏輯回歸逐步選擇變量的實現的續文-用R實現Stepwise。      R中可以使用step()函數進行變量逐步選擇，其通過AIC值來選擇模型和變量。在Python邏輯回歸逐步選擇變量的實現一文中，比較詳細的說明了逐步選擇變量的方法以及SAS的實現結果，並接著用python實現類似SAS的結果，本文採用R來實現類似的結果。