1 分别使用嶺回歸和Lasso解決薛毅書第279頁例6.10的回歸問題
例6.10的問題如下：
輸入例題中的資料，生成資料集,并做簡單線性回歸，檢視效果
cement <- data.frame(X1 = c(7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10), X2 = c(26,
29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68), X3 = c(6, 15, 8, 8, 6,
9, 17, 22, 18, 4, 23, 9, 8), X4 = c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26,
34, 12, 12), Y = c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1,
115.9, 83.8, 113.3, 109.4))
cement
## X1 X2 X3 X4 Y
## 1 7 26 6 60 78.5
## 2 1 29 15 52 74.3
## 3 11 56 8 20 104.3
## 4 11 31 8 47 87.6
## 5 7 52 6 33 95.9
## 6 11 55 9 22 109.2
## 7 3 71 17 6 102.7
## 8 1 31 22 44 72.5
## 9 2 54 18 22 93.1
## 10 21 47 4 26 115.9
## 11 1 40 23 34 83.8
## 12 11 66 9 12 113.3
## 13 10 68 8 12 109.4
lm.sol <- lm(Y ~ ., data = cement)
summary(lm.sol)
##
## Call:
## lm(formula = Y ~ ., data = cement)
## Residuals:
## Min 1Q Median 3Q Max
## -3.175 -1.671 0.251 1.378 3.925
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.405 70.071 0.89 0.399
## X1 1.551 0.745 2.08 0.071 .
## X2 0.510 0.724 0.70 0.501
## X3 0.102 0.755 0.14 0.896
## X4 -0.144 0.709 -0.20 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 2.45 on 8 degrees of freedom
## Multiple R-squared: 0.982, Adjusted R-squared: 0.974
## F-statistic: 111 on 4 and 8 DF, p-value: 4.76e-07
# 從結果看，截距和自變量的相關系數均不顯著。
# 利用car包中的vif（）函數檢視各自變量間的共線情況
library(car)
vif(lm.sol)
## X1 X2 X3 X4
## 38.50 254.42 46.87 282.51
# 從結果看，各自變量的VIF值都超過10，存在多重共線性，其中，X2與X4的VIF值均超過200.
plot(X2 ~ X4, col = "red", data = cement)
接下來，利用MASS包中的函數lm.ridge()來實作嶺回歸。下面的計算試了151個lambda值，最後選取了使得廣義交叉驗證GCV最小的那個。
library(MASS)
## Attaching package: 'MASS'
## The following object is masked _by_ '.GlobalEnv':
## cement
ridge.sol <- lm.ridge(Y ~ ., lambda = seq(0, 150, length = 151), data = cement,
model = TRUE)
names(ridge.sol) # 變量名字
## [1] "coef" "scales" "Inter" "lambda" "ym" "xm" "GCV" "kHKB"
## [9] "kLW"
ridge.sol$lambda[which.min(ridge.sol$GCV)] ##找到GCV最小時的lambdaGCV
## [1] 1
ridge.sol$coef[which.min(ridge.sol$GCV)] ##找到GCV最小時對應的系數
## [1] 7.627
par(mfrow = c(1, 2))
# 畫出圖形，并作出lambdaGCV取最小值時的那條豎直線
matplot(ridge.sol$lambda, t(ridge.sol$coef), xlab = expression(lamdba), ylab = "Cofficients",
type = "l", lty = 1:20)
abline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])
# 下面的語句繪出lambda同GCV之間關系的圖形
plot(ridge.sol$lambda, ridge.sol$GCV, type = "l", xlab = expression(lambda),
ylab = expression(beta))
par(mfrow = c(1, 1))
# 從上圖看，lambda的選擇并不是那麼重要，隻要不離lambda=0太近就沒有多大差别。
# 下面利用ridge包中的linearRidge()函數進行自動選擇嶺回歸參數
library(ridge)
mod <- linearRidge(Y ~ ., data = cement)
summary(mod)
## linearRidge(formula = Y ~ ., data = cement)
## Estimate Scaled estimate Std. Error (scaled) t value (scaled)
## (Intercept) 83.704 NA NA NA
## X1 1.292 26.332 3.672 7.17
## X2 0.298 16.046 3.988 4.02
## X3 -0.148 -3.279 3.598 0.91
## X4 -0.351 -20.329 3.996 5.09
## Pr(>|t|)
## (Intercept) NA
## X1 7.5e-13 ***
## X2 5.7e-05 ***
## X3 0.36
## X4 3.6e-07 ***
## Ridge parameter: 0.01473, chosen automatically, computed using 2 PCs
## Degrees of freedom: model 3.01 , variance 2.84 , residual 3.18
# 從模型運作結果看，測嶺回歸參數值為0.0147，各自變量的系數顯著想明顯提高（除了X3仍不顯著）
最後，利用Lasso回歸解決共線性問題
library(lars)
## Loaded lars 1.2
x = as.matrix(cement[, 1:4])
y = as.matrix(cement[, 5])
(laa = lars(x, y, type = "lar")) #lars函數值用于矩陣型資料
## lars(x = x, y = y, type = "lar")
## R-squared: 0.982
## Sequence of LAR moves:
## X4 X1 X2 X3
## Var 4 1 2 3
## Step 1 2 3 4
# 由此可見，LASSO的變量選擇依次是X4，X1，X2，X3
plot(laa) #繪出圖