R做線性回歸

本文轉載自「R語言」，己獲授權。

筆者邀請您，先思考：

1 線性回歸是什麼？

2 線性回歸怎麼應用？

本文解釋了如何在R中運行線性回歸。本教程將介紹線性回歸的假設以及如果假設不滿足如何處理。 它還包括擬合模型和計算模型性能指標以檢查線性回歸模型的性能。 線性回歸是最流行的統計技術之一。 它已被使用了三十多年。 它幾乎在每個領域都被廣泛接受，因為它很容易理解線性回歸的輸出。

線性回歸

它是一種發現稱為一個連續的因變數或者目標變數與一個或者多個（連續或者不連續）的自變數之間的關係的方法論。

這是一條直線曲線。 在上圖中，對角紅線是一條回歸線，也稱為最佳擬合直線。 點與回歸線之間的距離是誤差（errors）。線性回歸旨在通過最小化點與回歸線之間的垂直距離的平方和來找到最佳擬合直線。

線性回歸需要因變數是連續的，即數值（無類別或組）。

簡單與元線性回歸

當只有一個自變數時，線性回歸是簡單的線性回歸（一元線性回歸）。 而多元線性回歸會有多個自變數。

回歸方程

解釋：

當所有自變數變數（Xs）等於0且時，b0是因變數（Y）的期望平均值的截距。b1是斜率，它代表因變數（Y）的變化量，如果我們改變X1一個單位保持其他變數不變。

觀察到的（實際）因變數值與從回歸線預測的因變數值之間的差異。

演算法

線性回歸基於最小二乘估計，其表示回歸係數（估計值）應該以使每個觀察到的響應與其擬合值的平方距離的總和最小化的方式來選擇。

最小樣本量

線性回歸分析中需要每個自變數5個案例。

線性回歸分析的假設

1.線性關係：線性回歸需要依賴因變數和自變數之間的線性關係。

2.殘差的正態性：線性回歸需要殘差應該正態分布。

3.同方差性：線性回歸假定所有預測的因變數值的殘差大致相等。換句話說，它意味著誤差方差的不變性。

4.沒有離群值問題

5.多重共線性：這意味著自變數之間有很高的相關性。線性回歸模型不得面臨多重共線性問題。

6.誤差的獨立性 - 無自相關

它指出，與一個觀測相關的誤差與任何其他觀測的誤差都不相關。當您使用時間序列數據時，這是一個問題。假設你從八個不同的地區收集了勞動者的數據。每個地區內的勞動力可能更傾向於來自不同地區的勞動力，也就是說，他們的錯誤不是獨立的。

線性回歸假設目標或因變數是正態分布的。 正態分布與高斯分布相同。 它使用高斯家族的連接函數。

標準化係數

標準化或標準化係數（估計）的概念出現在預測變數（又稱自變數）以不同單位表示時。 假設你有3個獨立變數 -年齡，身高和體重。 變數「年齡」以年表示，身高以厘米為單位，體重以千克為單位。如果我們需要根據非標準化係數對這些預測因子進行排序，那麼它們將不會是一個公平的比較，因為這些變數的單位並不相同。

標準化係數（或估計）主要用於排列預測變數（或自變數或解釋變數），因為它消除了自變數和因變數的測量單位）。 我們可以用標準化係數的絕對值對自變數進行排序。 最重要的變數將具有最大的標準化係數絕對值。

標準化係數的解釋

1.25的標準化係數值表示自變數中一個標準偏差的變化導致因變數增加1.25個標準偏差。

模型性能測量1. R平方

它用模型中的所有自變數來解釋因變數的方差比例。 它假定模型中的每個獨立變數都有助於解釋因變數的變化。 事實上，有些變數不會影響因變數，它們也無助於建立一個好的模型。

上述方程式的分子，yi-hat是預測值。 Y的平均值出現在分母中。

R平方越高，模型越適合您的數據。 在心理調查或研究中，我們通常發現低R平方值低於0.5。 這是因為我們試圖預測人類行為，預測人類並不容易。 在這些情況下，如果您的R平方值很低，但您具有統計學上顯著的獨立變數（又稱預測變數），您仍然可以生成關於預測變數值中的變化如何與響應值變化相關聯的見解。

R平方可否為負？

是的，當水平線比您的模型更好地解釋數據時。 它主要發生在你不包括截距的情況下。 沒有截距，在預測目標變數方面，回歸可能會比樣本均值差。 這不僅是因為沒有截距。 即使包含截距，它也可能是負的。

在數學上，當模型的誤差平方大於水平線上的總平方和時，這是可能的。

R-squared = 1 - [(Sum of Square Error)/(Total Sum of Square)]

2.調整R平方

它衡量的只是那些真正影響因變數的自變數所解釋的方差比例。 它懲罰你添加不影響因變數的自變數。

調整後的R-Squared比R-squared更重要

每次向模型添加自變數時，即使自變數不顯著，R平方也會增加。 它永不衰落。 而調整R平方僅在自變數顯著並且影響因變數時增加。

3. RMSE（均方根誤差）

它解釋了實際數據點與模型預測值的接近程度。 它測量殘差的標準差。

在上面的公式中，yi是因變數的實際值，yi-hat是預測值，n是樣本大小。

RMSE與因變數具有相同的單位。 例如，因變數是以美元計量的銷售額。 假設這種銷售模式的RMSE出來21.我們可以說它是21美元。

什麼是好的RMSE分數？

沒有關於好或壞RMSE評分的拇指規則。 這是因為它依賴於你的因變數。如果您的目標變數介於0到100之間。RMS的0.5可以認為是好的，但同樣的0.5 RMSE可以被認為是差分，如果因變數的範圍從0到10.因此，通過簡單地查看 在價值。

較低的RMSE值表示更合適。 RMSE是衡量模型如何準確預測響應的一個很好的衡量標準，如果模型的主要目的是預測，那麼它就是擬合最重要的標準。

如果您已經建立了良好的模型，那麼您的訓練和測試集的RMSE應該非常相似。 如果測試集的RMSE遠高於訓練集的RMSE，那麼很可能是您的數據嚴重過擬合，即您創建了一個在樣本中運行良好的模型，但在樣本外的測試集上幾乎沒有預測價值。

RMSE與MAE

與平均絕對誤差（MAE）相比，RMSE嚴格懲罰大的誤差。

R平方是個比例並且沒有與目標變數相關聯的單位，而RMSE具有與目標變數相關的單位。 因此，R平方是擬合的相對度量，RMSE是擬合的絕對度量。

下面的代碼涵蓋了對模型性能的假設測試和評估：

1. 數據準備

2. 多重共線性檢驗

3. 多重共線性的處理

4. 檢查自相關

5. 檢查異常值

6. 檢查異方差

7. 殘差的正態性檢驗

8. 前進，後退和逐步選擇

9. 計算RMSE

10. 因變數的Box Cox變換

11. 手動計算R平方和調整R平方

12. 計算殘差和預測值

13. 計算標準化係數

R腳本：線性回歸

理論部分結束了。 讓我們用R -

載入所需的包

1. library(ggplot2)

2. library(car)

3. library(caret)

4. library(corrplot)

確保上面列出的軟體包已經安裝並載入到R中。如果它們尚未安裝，則需要使用命令install.packages（「package-name」）進行安裝。

我們將使用cars包中的mtcars數據集。 這些數據摘自Motor Trend US雜誌，包括32種汽車的燃料消耗和10個方面的汽車設計和性能。

1. data(mtcars)

2. str(mtcars)

以上變數的變數描述將在下面列出它們各自的變數名稱。

數據摘要

在這個數據集中，mpg是一個目標變數。 使用head（）函數查看前6行數據。

1. head(mtcars)

要查看變數的分布情況，請使用summary（）函數。

1. summary(mtcars)

數據準備

確保分類變數存儲為因子。 在下面的程序中，我們將變數轉換為因子。

1. mtcars$am <> as.factor(mtcars$am)

2. mtcars$cyl <> as.factor(mtcars$cyl)

3. mtcars$vs <> as.factor(mtcars$vs)

4. mtcars$gear <> as.factor(mtcars$gear)

識別和修正共線性

在這一步中，我們正在識別彼此高度相關的自變數。 由於mpg是一個因變數，我們將在下面的代碼中刪除它。

1. mtcars_a <> subset(mtcars, select = -c(mpg))

2. numericData <> mtcars_a[sapply(mtcars_a, is.numeric)]

3. descrCor <> cor(numericData)

4. print(descrCor)

1. corrplot(

2. descrCor,

3. order = "FPC",

4. method = "color",

5. type = "lower",

6. tl.cex = 0.7,

7. tl.col = rgb(0,0,0)

8. )

1. highlyCorrelated <> findCorrelation(descrCor, cutoff=0.7)

2. highlyCorCol <> colnames(numericData)[highlyCorrelated]

3. dat3 <> mtcars[, -which(colnames(mtcars) %in% highlyCorCol)]

4. dim(dat3)

有三個變數「hp」「disp」「wt」被發現是高度相關的。 我們已經刪除它們以避免共線。 現在，我們有7個獨立變數和1個因變數。

在這一步，我們正在建立多元線性回歸模型。

1. fit <> lm(mpg ~ ., data=dat3)

2. summary(fit)

3. summary(fit)$coeff

4. anova(fit)

5. par(mfrow=c(2,2))

6. plot(fit)

查看線性回歸模型的係數和ANOVA表

性回歸模型測試估計等於零的零假設。 具有小於0.05的p值的獨立變數意味著我們拒絕5％顯著性水平的零假設。這意味著該變數的係數不等於0.大p值意味著變數對預測目標變數沒有意義。

線性回歸模型的關鍵圖如下所示 :-

計算模型性能度量

1. summary(fit)$r.squared

2. summary(fit)$adj.r.squared

3. AIC(fit)

4. BIC(fit)

更高的R平方和調整的R平方值，更好的模型。 然而，更低的AIC和BIC得分，更好的模型。

理解AIC和BIC

AIC和BIC是擬合度的衡量標準。 他們懲罰複雜的模型。 換句話說，它會懲罰更多的估計參數。 它相信一個概念，即一個具有較少參數的模型將被一個具有更多參數的模型要好。 一般來說，BIC比AIC更為免費參數懲罰模型。 兩個標準都取決於估計模型的似然函數L的最大值。

AIC值大致等於參數的數量減去整個模型的似然性。 假設你有兩個模型，AIC和BIC分數較低的模型更好。

變數選擇方法

有三種變數選擇方法 - 向前，向後，逐步。1.以單個變數開始，然後基於AIC（「前進」）一次添加一個變數2.從所有變數開始，基於AIC（』後退』）迭代地去除那些重要性低的變數3.雙向運行（』逐步』）

1. library(MASS)

2. step <> stepAIC(fit, direction="both")

3. summary(step)

4. step <> stepAIC(fit, direction="forward")

5. summary(step)

6. n <> dim(dat3)[1]

7. stepBIC <> stepAIC(fit,k=log(n))

8. summary(stepBIC)

在基於BIC執行逐步選擇之後，查看以上估計值。 變數已經減少，但調整後的R-Squared保持不變（稍微改善）。 AIC和BIC分數也下降，這表明一個更好的模型。

1. AIC(stepBIC)

2. BIC(stepBIC)

計算標準化係數

標準化係數有助於根據標準化估計值的絕對值排列預測值。 值越高，變數越重要。

1. #使用QuantPsyc包的lm.beta函數計算

2. library(QuantPsyc)

3. lm.beta(stepBIC)

5. #自定義函數計算

6. stdz.coff <> function (regmodel)

7. { b <> summary(regmodel)$coef[-1,1]

8. sx <> sapply(regmodel$model[-1], sd)

9. sy <> sapply(regmodel$model[1], sd)

10. beta <> b * sx / sy

11. return(beta)

12. }

14. std.Coeff <> data.frame(Standardized.Coeff = stdz.coff(stepBIC))

15. std.Coeff <> cbind(Variable = row.names(std.Coeff), std.Coeff)

16. row.names(std.Coeff) <> NULL

計算方差膨脹因子（VIF）

與獨立變數高度相關的情況相比，差異膨脹因子衡量的是係數的變化幅度。 它應該小於5。

1. vif(stepBIC)

測試其它假設

1. #Autocorrelation Test

2. durbinWatsonTest(stepBIC)

4. #Normality Of Residuals (Should be > 0.05)

5. res=residuals(stepBIC,type="pearson")

6. shapiro.test(res)

8. #Testing for heteroscedasticity (Should be > 0.05)

9. ncvTest(stepBIC)

11. #Outliers – Bonferonni test

12. outlierTest(stepBIC)

14. #See Residuals

15. resid = residuals(stepBIC)

17. #Relative Importance

18. install.packages("relaimpo")

19. library(relaimpo)

20. calc.relimp(stepBIC)

查看實際值和預測值

1. #See Predicted Value

2. pred <> predict(stepBIC,dat3)

3. #See Actual vs. Predicted Value

4. finaldata <> cbind(mtcars,pred)

5. print(head(subset(finaldata, select = c(mpg,pred))))

其它有用的函數

1. #Calculating RMSE

2. rmse = sqrt(mean((dat3$mpg - pred)^2))

3. print(rmse)

5. #Calculating Rsquared manually

6. y = dat3[,c("mpg")]

7. R.squared = 1 - sum((y-pred)^2)/sum((y-mean(y))^2)

8. print(R.squared)

10. #Calculating Adj. Rsquared manually

11. n = dim(dat3)[1]

12. p = dim(summary(stepBIC)$coeff)[1] - 1

13. adj.r.squared = 1 - (1 - R.squared) * ((n - 1)/(n-p-1))

14. print(adj.r.squared)

16. #Box Cox Transformation

17. library(lmSupport)

18. modelBoxCox(stepBIC)

K-fold交叉驗證

在下面的程序中，我們正在進行5倍交叉驗證。 在5倍交叉驗證中，數據被隨機分成5個相同大小的樣本。 在5個樣本中，隨機20％數據的單個樣本保留為驗證數據，其餘80％用作訓練數據。 然後這個過程重複5次，5個樣本中的每一個都只用作驗證數據一次。 稍後我們將結果平均。

1. library(DAAG)

2. kfold = cv.lm(data=dat3, stepBIC, m=5)

作者：Deepanshu原文鏈接：https://www.listendata.com/2015/09/linear-regression-with-r.html譯者：陸勤