本文對應《R語言實戰》第8章：回歸
回歸是一個廣義的概念，通指那些用一個或多個預測變量（也稱自變量或解釋變量）來預測響應變量（也稱因變量、效標變量或結果變量）的方法。通常，回歸分析可以用來挑選與相應變量相關的解釋變量，可以描述兩者的關系，也可以生成一個等式，通過解釋變量來預測響應變量。
回歸分析的各種變體
| 回歸類型 |
用途 |
| 簡單線性 |
用一個量化的解釋變量預測一個量化的響應變量 |
| 多項式 |
用一個量化的解釋變量預測一個量化的響應變量，模型的關系是n階多項式 |
| 多元線性 |
用兩個或多個量化的解釋變量預測一個量化的響應變量 |
| 多變量 |
用一個或多個解釋變量預測多個響應變量 |
| Logistic |
用一個或多個解釋變量預測一個類別型響應變量 |
| 泊松 |
用一個或多個解釋變量預測一個代表頻數的響應變量 |
| Cox比例風險 |
用一個或多個解釋變量預測一個事件（死亡、失敗或舊病復發）發生的時間 |
| 時間序列 |
對誤差項相關的時間序列數據建模 |
| 非線性 |
用一個或多個量化的解釋變量預測一個量化的響應變量，不過模型是非線性的 |
| 非參數 |
用一個或多個量化的解釋變量預測一個量化的響應變量，模型的形式源自數據形式，不事先設定 |
| 穩健 |
用一個或多個量化的解釋變量預測一個量化的響應變量，能抵御強影響點的干擾 |
本章主要介紹OLS（最小二乘回歸法）
統計假設：
正態性：對於固定自變量值，因變量值成正態分布
獨立性：因變量值之間相互獨立
線性：因變量與自變量之間線性相關
同方差性：因變量的方差不隨自變量的水平不同而變化
基本格式：
myfit <- lm(formula, data)
表達式中常用符號
| ~ |
分隔符號，左邊為響應變量，右邊為解釋變量 |
| + |
分隔預測變量 |
| : |
表示預測變量的交互項 e.g. x:z |
| * |
表示所有可能交互項的簡潔方式 e.g. x * z = x + z + x:z |
| ^ |
表示交互項達到某個次數 |
| . |
表示除因變量外的所有變量 |
| - |
表示從formula中移除某變量 |
| -1 |
刪除截距項（強制回歸直線通過原點） |
| I() |
從算術的角度來解釋括號中的元素（即在I中的表達式符號為算術符號） |
| function |
可以在表達式中用的數學函數 |
擬合模型的其他常用函數
| summary() |
展示擬合模型的詳細結果 |
| coefficients() |
列出擬合模型的模型參數 |
| confint() |
提供模型參數的置信區間（默認95%） |
| fitted() |
列出擬合模型的預測值 |
| residuals() |
列出擬合模型的殘差值 |
| anova() |
生成一個擬合模型的方差分析表，或者比較兩個或更多擬合模型的方差分析表 |
| vcov() |
列出模型參數的協方差矩陣 |
| AIC() |
輸出赤池信息統計量 |
| plot() |
生成評價擬合模型的診斷圖 |
| predict() |
用擬合模型對新的數據集預測響應變量值 |
代碼示例
#簡單線性回歸 fit <- lm(weight ~ height, data = women) summary(fit) #多項式回歸 fit2 <- lm(weight ~ height + I(height ^ 2), data = women) summary(fit2) fit3 <- lm(weight ~ height + I(height ^ 2) + I(height ^ 3), data = women) summary(fit3) #多元線性回歸 fit4 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states) summary(fit4) #有交互項的多元線性回歸 fit5 <- lm(mpg ~ hp + wt + hp:wt) summary(fit5)
回歸診斷：
標准方法：plot(fit)
四張圖：Residuals vs Fitted(殘差圖與擬合圖), Normal Q-Q(正態Q-Q圖), Scale-Location(位置尺度圖), Residuals vs Leverage(殘差與杠桿圖)
診斷統計假設方式：
正態性：Normal Q-Q圖上的點在45度直線上
獨立性：無法從圖中讀出，但可以收集的數據中驗證
線性：Residuals vs Fitted隨機分布（若有明顯趨勢，則暗示回歸模型有相應的趨勢分布）
同方差性：Scale-Location水平線周圍的點隨機分布
Residuals vs Leverage圖可以鑒別離群點、高杠桿值點和強影響點（由於該圖可讀性差，通過其他方式完成這個工作）
改進的方法：car包
| qqplot() |
分位數比較圖 |
| durbinWatsonTest() |
對誤差自相關性做Durbin-Watson檢驗 |
| crPlots() |
成分與殘差圖 |
| ncvTest() |
對非恆定的誤差方差做得分檢驗 |
| spreadLevelPlot() |
分散水平檢驗 |
| outlierTest() |
Bonferroni離群點檢驗 |
| avPlots() |
添加的變量圖形 |
| inluencePlot() |
回歸影響圖 |
| scatterplot() |
增強的散點圖 |
| scatterplotMatrix() |
增強的散點圖矩陣 |
| vif() |
方差膨脹因子 |
正態性：
qqPlot()函數提供了更為精確的正態假設檢驗方法：畫出了n-p-1個自由度的t分布下的學生化殘差（studentized residual, 也稱學生化刪除殘差或折疊化殘差）圖形，其中n是樣本大小，p是回歸參數的數目（包括截距項）
#id.method = “identify”選項能夠交互式繪圖，鼠標單擊圖形內的點，將會標注函數中label選項的設定值 #simulate = TRUE時，95%的置信區間將會用參數自助法 library(car) fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states) qqPlot(fit, labels = row.names(states), id.method = “identify”, simulate = TRUE, main = “Q-Q Plot”)
誤差的獨立性：
最好的方法是依據收集數據方式的先驗知識。car包中Durbin-Watson檢驗函數也可以檢測誤差的序列相關性（p值不顯著說明無自相關性）。
線性：
成分殘差圖（component plus residual plot，或稱偏殘差圖，partial residual plot）
若圖形存在非線性，則說明你可能對預測變量的函數形式建模不夠充分，那么就需要添加一些曲線成分，或對一個或多個變量進行變換，或用其他回歸變體形式而不是線性回歸。
同方差性：
ncvTest()函數生成一個計分檢驗，零假設為誤差方差不變，備擇假設為誤差方差隨着擬合值水平的變化而變化。（若檢驗顯著，則說明存在異方差性）
spreadLevelPlot()函數創建一個添加了最佳擬合曲線的散點圖，展示標准化殘差絕對值與擬合值的關系。另外會輸出一個建議冪次變換（suggested power transformation），含義是：經過p次冪變換，非恆定的方差將會平穩。
線性模型假設的綜合驗證：
使用gvlma包中的函數gvlma()
library(gvlma) gvmodel <- gvlma(fit) summary(gvmodel)
輸出項可以看到滿足OLS回歸模型的所有統計假設，若p<0.05可認為違反假設條件
多重共線性：
兩個預測變量之間存在高度的相關性，造成模型參數的置信區間過大，使單個系數解釋起來很困難。（簡單來說，回歸系數測量的是當其他預測變量不變時，某個預測變量與響應變量之間的關系，而有兩個預測變量之間相關時，就變成了假定某變量不變，測量響應變量與該變量之間的關系）
多重共線性可用統計量VIF（variance inflation factor, 方差膨脹因子）進行檢測。VIF的平方根表示變量回歸參數的置信區間能膨脹為與模型無關的預測變量的程度。car包中的vif()函數可以計算。一般原則下，就表明存在多重共線性問題。
library(car) vif(fit) sqrt(vif(fit)) > 2 #返回各變量及布爾值，TRUE代表存在問題，FALSE代表不存在問題
異常觀測值：離群點、高杠桿值、強影響點
離群點：指那些模型預測效果不佳的觀測點（很大的殘差）。鑒別：Q-Q圖中落在置信區間外的點，或標准化殘差值大於2或小於-2的點。
car包中的函數outlierTest()可以檢測到單個離群點（p值顯著），要檢驗其他離群點，需要先將第一步的離群點刪除重新檢測。
高杠桿值點：與其他預測變量有關的離群點，即許多異常的預測變量值組合起來的，與響應變量值沒有關系。
鑒別：通過帽子統計量（hat statistic）判斷。對於一個給定的數據集，帽子均值為p/n，其中，p是模型估計的參數數目（包含截距項），n為樣本量。一般來說，觀測點的帽子值大於帽子均值的2倍或3倍，即可認定為高杠桿值點。
高杠桿值可能是強影響點，也可能不是，這要看它們是否是離群點
#自定義畫出帽子值分布，並給出閾值線
hat.plot <- function (fit) {
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit), main = “Index Plot of Hat Values”)
abline(h = c(2, 3)*p/n, col = “red”, lty = 2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)
強影響點：即對模型參數估計值影響有些比例失衡的點。
檢測：Cook距離（或稱D統計量）或者變量添加圖（added variable plot）。一般來說，Cook’s D值大於4/(n-k-1),則表明它是強影響點，其中n為樣本量大小，k是預測變量數目。（注：如果以D=1為判別標准可能更具有一般性）
#繪制Cook’s D圖形 cutoff <- 4/(nrow(states) – length(fit$coefficients) - 2) plot(fit, which = 4, cook.levels = cutoff) abline(h = cutoff, lty = 2, col = “red”)
Cook’s D圖有助於鑒別強影響點，但是並不能提供關於這些點如何影響模型的信息。變量添加圖彌補了這個缺陷。
library(car) avPlots(fit, ask = FALSE, onepage = TRUE, id.method = “identify”)
將離群點、高杠桿值和強影響點信息整合到一幅圖形中：
library(car) influencePlot(fit, id.method = “identify”, main = “Influence Plot”, sub = “Circle size is proportional to Cook’s distance”)
注：圖中，縱坐標超過+2或小於-2的可被認為是離群點，橫坐標超過0.2或0.3的有高杠桿值，圓圈大小與影響成比例，圓圈很大的點可能是對模型參數的估計造成的不成比例影響的強影響點。
改進措施：針對違背回歸假設的問題
刪除觀測點；變量變換；添加或刪除變量；使用其他回歸方法
刪除觀測點：刪除離群點或強影響點。謹慎。若是數據記錄錯誤、沒有遵守規程等的離群點可以刪除，否則發掘為何該觀測點不同於其他點也是研究的方向之一。
變量變換：模型不符合正態性、線性或者同方差性假設時，一個或多個變量變換可以改善效果。模型違反正態假設時，通常對響應變量嘗試變換；違反線性假設時，對預測變量變換有效，同時也可以改善異方差性。需要說明的是，變換后的變量需要有意義，否則很難解釋。
增刪變量：刪除變量可以處理多重共線性。但是如果僅是預測並不需要處理，處理多重共線性最重要的是對預測變量進行解釋。
選擇“最佳”的回歸模型
#檢驗嵌套模型的擬合優度 fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states) fit2 <- lm(Murder ~ Population + Illiteracy, data = states) anova(fit2, fit1)
如果檢驗不顯著（p > 0.05），則使用較少變量的線性模型。
使用AIC（Akaike Information Criterion，赤池信息准則）比較模型：AIC值越小的模型越優先選擇。
變量選擇：逐步回歸法（stepwise method）；全子集回歸（all-subsets regression）
逐步回歸法：向前逐步回歸（forward stepwise）向后逐步回歸（backward stepwise）向前向后逐步回歸（stepwise stepwise）
library(MASS) fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states) stepAIC(fit1, direction = “backward”)
逐步回歸法其實存在爭議，雖然它可能會找到一個好模型，但是不能保證模型就是最佳模型，因為不是每一個可能的模型都被評價了。為了克服這個限制，便有了全子集回歸法。
全子集回歸法：所有可能的模型都會被檢驗，可用leaps包中的regsubsets()函數實現。一般來說，變量的自動選擇應該被看作是對模型選擇的一種輔助方法，而不是直接方法。擬合效果佳而沒有意義的模型對你毫無幫助，主題背景知識的理解才能最終指引你獲得理想的模型。
交叉驗證：將一定比例的數據挑選出來作為訓練樣本，另外的樣本作保留樣本，先在訓練樣本上獲得回歸方程，然后在保留樣本上做預測。
在k折交叉驗證中，樣本被分為k個子樣本，輪流將k-1個子樣本組合作為訓練集，另外一個子樣本作為保留集，獲得k個預測方程，記錄k個保留樣本的預測表現結果，求均值。
相對重要性：比較標准化的回歸系數，表示當其他預測變量不變時，該預測變量一個標准差的變化可引起的響應變量的預期變化（以標准差單位度量）。在進入回歸分析之前，可用scale()函數將數據標准化為均值為0、標准差為1的數據集，這樣用R回歸即可獲得標准化的回歸系數。（注意scale()函數返回的是一個矩陣，而lm()函數要求一個數據框）
zstates <- as.data.frame(scale(states)) zfit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data = zstates) coef(zfit)
另外可以根據相對權重判斷相對重要性。比較各個回歸系數，最大的就是最重要的。