數據分析中常見的七種回歸分析以及R語言實現（三）---嶺回歸

02-03

我們在回歸分析的時候，古典模型中有一個基本的假定就是自變數之間是不相關的，但是如果我們在擬合出來的回歸模型出現了自變數之間高度相關的話，可能對結果又產生影響，我們稱這個問題為多重共線性，多重共線性又分為兩種，一種是完全多重共線性，還有一種是不完全多重共線性；

產生的原因有幾個方面

1、變數之間存在內部的聯繫

2、變數之間存在共同的趨勢等

造成的後果分兩部分

完全多重共線性造成的後果

1、當自變數線性相關的時候，參數將無法唯一確定，參數的方差將趨近於無窮大，這時候無法使用最小二乘法

不完全多重共線性造成的後果

1、參數估計量的方差隨著多重共線性的嚴重程度的增加而增加，但是參數是可以估計的

2、進行統計檢驗時容易刪除掉重要解釋變數

因為當多重共線性的時候容易造成自變數對因變數不顯著，從模型中錯誤的剔除，這樣容易刪除重要解釋變數的設定；

3、參數的置信區間明顯擴大

因為由於存在多重共線性。我們的參數估計都有較大的標準差，因此參數真值的置信區間也將增大

那麼我們怎麼去判斷一個模型上存在多重共線性呢？

根據經驗表明，多重共線性存在的一個標誌就是就模型存在較大的標準差，和較小的T統計量，如果一個模型的可決係數R^2很大，F檢驗高度限制，但偏回歸係數的T檢驗幾乎都不顯著，那麼模型很可能是存在多重共線性了。因為通過檢驗，雖然各個解釋變數對因變數的共同影響高度顯著，但每個解釋變數的單獨影響都不顯著，我們無法判斷哪個解釋變數對被解釋變數的影響更大

1、可以利用自變數之間的簡單相關係數檢驗

這個方法是一個簡便的方法，一般而言，如果每兩個解釋變數的簡單相關係數一般較高，則可以認為是存在著嚴重的多重共線性

2、方差膨脹因子

在回歸中我們用VIF表示方差膨脹因子

表達式 VIF=1/(1-R^2)

隨著多重共線性的嚴重程度增強，方差膨脹因子會逐漸的變大，一般的當VIF>=10的時候，我們就可以認為存在嚴重多重共線性；

在R語言中car包中的vif()函數可以幫我們算出這個方差膨脹一找你

這就介紹這兩個了，其實還有好多方法，大家可以可以私底下查，或者和我一起交流；

多重共線性的解決辦法

因為存在多重共線性，我們還是擬合模型的；當然會有解決辦法，這裡我就介紹一下常用的方法嶺回歸；其他的方法也有，這裡就不說了；

這裡就說說大概的思想，具體推導的步驟這裡就不寫，有興趣的可以網上查查；在多重共線性十分嚴重下，兩個共線變數的係數之間的二維聯合分布是一個山嶺曲面，曲面上的每一個點對應一種殘差平方和，點的位置越高，相應的殘差平方和越小。因此山嶺最高點和殘差平方和的最小值相對應，相應的參數值便是參數的最小二乘法估計值，但由於多重共線性的存在最小二乘法估計量已經不適用，一個自然的想法就是應尋找其他的更適合的估計量，這種估計量既要具有較小的方差，又不能使殘差平方和過分偏離其極小值。在參數的聯合分布曲面上，能滿足這種要求的點只能沿著山嶺尋找，這就是嶺回歸法；

這個方法實質是犧牲了無偏性來尋求參數估計的最小方差性；

缺點：通常嶺回歸方程的R平方值會稍低於普通回歸分析，但回歸係數的顯著性往往明顯高於普通回歸，在存在共線性問題和病態數據偏多的研究中有較大的實用價值

R語言建模

這裡使用可能要使用到car和MASS，由於謝佳標老師已經寫了詳細的過程，這裡我就全程照搬了，偷了個懶，寫個代碼過程其實也有些累的；

原鏈接用R建立嶺回歸和lasso回歸 - jiabiao1602的專欄 - 博客頻道 - CSDN.NET

1 分別使用嶺回歸和Lasso解決薛毅書第279頁例6.10的回歸問題

cement <- data.frame(X1 = c(7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10), X2 = c(26, n 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68), X3 = c(6, 15, 8, 8, 6, n 9, 17, 22, 18, 4, 23, 9, 8), X4 = c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, n 34, 12, 12), Y = c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, n 115.9, 83.8, 113.3, 109.4))ncementn## X1 X2 X3 X4 Yn## 1 7 26 6 60 78.5n## 2 1 29 15 52 74.3n## 3 11 56 8 20 104.3n## 4 11 31 8 47 87.6n## 5 7 52 6 33 95.9n## 6 11 55 9 22 109.2n## 7 3 71 17 6 102.7n## 8 1 31 22 44 72.5n## 9 2 54 18 22 93.1n## 10 21 47 4 26 115.9n## 11 1 40 23 34 83.8n## 12 11 66 9 12 113.3n## 13 10 68 8 12 109.4nlm.sol <- lm(Y ~ ., data = cement)nsummary(lm.sol)n## n## Call:n## lm(formula = Y ~ ., data = cement)n## n## Residuals:n## Min 1Q Median 3Q Max n## -3.175 -1.671 0.251 1.378 3.925 n## n## Coefficients:n## Estimate Std. Error t value Pr(>|t|) n## (Intercept) 62.405 70.071 0.89 0.399 n## X1 1.551 0.745 2.08 0.071 .n## X2 0.510 0.724 0.70 0.501 n## X3 0.102 0.755 0.14 0.896 n## X4 -0.144 0.709 -0.20 0.844 n## ---n## Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1n## n## Residual standard error: 2.45 on 8 degrees of freedomn## Multiple R-squared: 0.982, Adjusted R-squared: 0.974 n## F-statistic: 111 on 4 and 8 DF, p-value: 4.76e-07n# 從結果看，截距和自變數的相關係數均不顯著。n# 利用car包中的vif（）函數查看各自變數間的共線情況nlibrary(car)nvif(lm.sol)n## X1 X2 X3 X4 n## 38.50 254.42 46.87 282.51n# 從結果看，各自變數的VIF值都超過10，存在多重共線性，其中，X2與X4的VIF值均超過200.nplot(X2 ~ X4, col = "red", data = cement)n

接下來，利用MASS包中的函數lm.ridge()來實現嶺回歸。下面的計算試了151個lambda值，最後選取了使得廣義交叉驗證GCV最小的那個。

library(MASS)n## n## Attaching package: MASSn## n## The following object is masked _by_ .GlobalEnv:n## n## cementnridge.sol <- lm.ridge(Y ~ ., lambda = seq(0, 150, length = 151), data = cement, n model = TRUE)nnames(ridge.sol) # 變數名字n## [1] "coef" "scales" "Inter" "lambda" "ym" "xm" "GCV" "kHKB" n## [9] "kLW"nridge.sol$lambda[which.min(ridge.sol$GCV)] ##找到GCV最小時的lambdaGCVn## [1] 1nridge.sol$coef[which.min(ridge.sol$GCV)] ##找到GCV最小時對應的係數n## [1] 7.627npar(mfrow = c(1, 2))n# 畫出圖形，並作出lambdaGCV取最小值時的那條豎直線nmatplot(ridge.sol$lambda, t(ridge.sol$coef), xlab = expression(lamdba), ylab = "Cofficients", n type = "l", lty = 1:20)nabline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])n# 下面的語句繪出lambda同GCV之間關係的圖形nplot(ridge.sol$lambda, ridge.sol$GCV, type = "l", xlab = expression(lambda), n ylab = expression(beta))nabline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])n

par(mfrow = c(1, 1))n# 從上圖看，lambda的選擇並不是那麼重要，只要不離lambda=0太近就沒有多大差別。n# 下面利用ridge包中的linearRidge()函數進行自動選擇嶺回歸參數nlibrary(ridge)nmod <- linearRidge(Y ~ ., data = cement)nsummary(mod)n## n## Call:n## linearRidge(formula = Y ~ ., data = cement)n## n## n## Coefficients:n## Estimate Scaled estimate Std. Error (scaled) t value (scaled)n## (Intercept) 83.704 NA NA NAn## X1 1.292 26.332 3.672 7.17n## X2 0.298 16.046 3.988 4.02n## X3 -0.148 -3.279 3.598 0.91n## X4 -0.351 -20.329 3.996 5.09n## Pr(>|t|) n## (Intercept) NA n## X1 7.5e-13 ***n## X2 5.7e-05 ***n## X3 0.36 n## X4 3.6e-07 ***n## ---n## Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1n## n## Ridge parameter: 0.01473, chosen automatically, computed using 2 PCsn## n## Degrees of freedom: model 3.01 , variance 2.84 , residual 3.18n# 從模型運行結果看，測嶺回歸參數值為0.0147，各自變數的係數顯著想明顯提高（除了X3仍不顯著）n最後，利用Lasso回歸解決共線性問題nlibrary(lars)n## Loaded lars 1.2nx = as.matrix(cement[, 1:4])ny = as.matrix(cement[, 5])n(laa = lars(x, y, type = "lar")) #lars函數值用於矩陣型數據n## n## Call:n## lars(x = x, y = y, type = "lar")n## R-squared: 0.982 n## Sequence of LAR moves:n## X4 X1 X2 X3n## Var 4 1 2 3n## Step 1 2 3 4n# 由此可見，LASSO的變數選擇依次是X4，X1，X2，X3nplot(laa) #繪出圖n

summary(laa) #給出Cp值n## LARS/LARn## Call: lars(x = x, y = y, type = "lar")n## Df Rss Cpn## 0 1 2716 442.92n## 1 2 2219 361.95n## 2 3 1918 313.50n## 3 4 48 3.02n## 4 5 48 5.00n# 根據課上對Cp含義的解釋（衡量多重共線性，其值越小越好），我們取到第3步，使得Cp值最小，也就是選擇X4，X1，X2這三個變數n