Lasso and Ridge
Liu, Tianxin
2019/8/22
话不多说，先看一组数据
head(food, 10)## Y X1 X2 X3 X4 X5
## 1: 98.45 560.20 153.20 6.53 1.23 1.89
## 2: 100.70 603.11 190.00 9.12 1.30 2.03
## 3: 102.80 668.05 240.30 8.10 1.80 2.71
## 4: 133.95 715.47 301.12 10.10 2.09 3.00
## 5: 140.13 724.27 361.00 10.93 2.39 3.29
## 6: 143.11 736.13 420.00 11.85 3.90 5.24
## 7: 146.15 748.91 491.76 12.28 5.13 6.83
## 8: 144.60 760.32 501.00 13.50 5.47 8.36
## 9: 148.94 774.92 529.20 15.29 6.09 10.07
## 10: 158.55 785.30 552.72 18.10 7.97 12.57其中，
\(Y\):市粮食销售量(万吨/年),
\(X_1\):市常住人口数(万人),
\(X_2\):市人均年收入(元/年),
\(X_3\):市肉类销售量(万吨/年),
\(X_4\):市蛋类销售量(万吨/年),
\(X_5\):市鱼虾类销售量(万吨/年)。
我们来直观的看一下什么样的数据能说明变量间存在多重共线性，
首先，看单一自变量\(X\)和\(Y\)之间的关系，
pushViewport(viewport(layout = grid.layout(2,3)))
vplayout <- function(x,y){
viewport(layout.pos.row = x, layout.pos.col = y)
}
print(a, vp=vplayout(1,1))
print(b, vp=vplayout(1,2))
print(c, vp=vplayout(1,3))
print(d, vp=vplayout(2,1))
print(e, vp=vplayout(2,2))
我们可以猜测，每一个自变量\(X\)对\(Y\)都是正相关的。
然后我们看一下，单一自变量\(X\)对\(Y\)的贡献是否显著，
summary(linear_model_X1)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -90.9207421 19.32928786 -4.703781 5.109817e-04
## X1 0.3169248 0.02608088 12.151615 4.204073e-08summary(linear_model_X2)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 99.61349252 6.43124182 15.488998 2.691506e-09
## X2 0.08146954 0.01073788 7.587119 6.440540e-06summary(linear_model_X3)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 74.648244 8.2889893 9.005711 1.096718e-06
## X3 4.892712 0.5635782 8.681514 1.612586e-06summary(linear_model_X4)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.864721 5.9343305 18.344904 3.814392e-10
## X4 5.739752 0.8387557 6.843175 1.789484e-05summary(linear_model_X5)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 113.374731 6.0771328 18.655957 3.137309e-10
## X5 3.080811 0.5122999 6.013688 6.087610e-05
我们可以猜测，每一个自变量\(X\)对\(Y\)都是有显著贡献的。
如果我们将这些自变量\(X\)合在一起，对\(Y\)进行预测，有趣的事情就发生了。
summary(linear_model)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = food)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8657 -1.3698 -0.1279 2.9657 7.2070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.49656 30.00659 -0.117 0.9101
## X1 0.12533 0.05914 2.119 0.0669 .
## X2 0.07367 0.03788 1.945 0.0877 .
## X3 2.67759 1.25729 2.130 0.0658 .
## X4 3.45345 2.45085 1.409 0.1965
## X5 -4.49112 2.21486 -2.028 0.0771 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.72 on 8 degrees of freedom
## Multiple R-squared: 0.9704, Adjusted R-squared: 0.952
## F-statistic: 52.53 on 5 and 8 DF, p-value: 6.645e-06
从结果上看，决定系数\(R^2\)很高，这难道不是一个很优秀的模型吗？
但是，我们可以看出自变量\(X_5\)的回归系数，由正变负，已经违背了业务逻辑和生活常识，这样的预测是错误的。
这些现象都说明了自变量\(X\)之间具有多重共线性。
定义：多重共线性是指线性回归模型中的自变量\(X\)之间由于存在高度相关关系而使模型难以估计准确。
接下来，我们看一下这种自变量\(X\)之间得高度相关有哪几种判断方法
cor(food[,-1])## X1 X2 X3 X4 X5
## X1 1.0000000 0.8665519 0.8822931 0.8524491 0.8213054
## X2 0.8665519 1.0000000 0.9458957 0.9647730 0.9825321
## X3 0.8822931 0.9458957 1.0000000 0.9405058 0.9483613
## X4 0.8524491 0.9647730 0.9405058 1.0000000 0.9819792
## X5 0.8213054 0.9825321 0.9483613 0.9819792 1.0000000
从自变量之间的相关系数看，当\(r\in(-\infty,-0.8]\)或\(r\in[0.8,+\infty)\)的时候，我们可以称为自变量\(X\)之间高度相关。
vif(linear_model)## X1 X2 X3 X4 X5
## 8.715944 48.410630 15.417850 39.281292 105.034818
从方差膨胀因子\(VIF\)看，当\(VIF>10\)的时候，我们可以称为自变量\(X\)之间高度相关。
kappa(food[,-1])## [1] 1172.407
从条件数\(\kappa\)来看，当\(\kappa>1000\)的时候，我们可以称为自变量\(X\)之间高度相关。
当自变量\(X\)出现高度相关的情况时，会造成\(|X^TX|\approx 0\)，从而使最小二乘法不成立。
为了解决这个问题，岭回归应运而生。
定义：岭回归是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于最小二乘法。
也就是说，岭回归=改良后的最小二乘法。
为了解决\(|X^TX|\approx 0\)的问题，岭回归加入了\(\lambda I\)使得\(|X^TX+\lambda I|\neq0\)可逆，最小二乘法成立。
艰难的抉择
对于最小二乘法而言，随着模型复杂度的提升，模型的效果就越好，即预测值与真实值的差距就越小；但是同时模型的回归系数越大，即模型的方差就越大。
对于岭回归而言，\(\lambda\)越小表示模型预测值与真实值的差距越小；\(\lambda\)越大表示\((X^TX+\lambda I)^{-1}\)越小，模型的方差越小。
所以岭回归的关键是找到一个合理的\(\lambda\)值来平衡模型的方差和偏差。
那么岭回归的最小化问题等价于寻找 \[\begin{equation}
\beta^{ridge} = argmin\big[\sum(Y-\beta X)^2+\lambda\sum\beta^2\big]
\end{equation}\]
即找到 \[\begin{equation}
\lambda\sum\beta^2\leq C
\end{equation}\]
其中\(C\)为常数，也就是说所有系数\(\beta\)平方和小于某个阈值\(C\)。
以两个变量为例， \[\begin{equation} Y=\left[ \begin{matrix} y_1\\ y_2\\ ...\\ y_n \end{matrix} \right], \beta=\left[ \begin{matrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{matrix} \right], X=\left[ \begin{matrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ ... & ... & ...\\ 1 & x_{n1} & x_{n2} \end{matrix} \right], \epsilon=\left[ \begin{matrix} \epsilon_1\\ \epsilon_2\\ ...\\ \epsilon_n \end{matrix} \right] \end{equation}\]
在没有岭回归的帮助时，残差平方和\(SSE\)可以表示为\(\beta_1\)和\(\beta_2\)的一个二次函数，数学上可以用一个抛物面表示。
从\(\beta_1\)、\(\beta_2\)平面理解，即为抛物面等高线在水平面的投影和圆的交点，如下图所示
所以，当\(\lambda=0\)时，得到最小二乘解。当参数\(\lambda\)趋向更大时，回归系数\(\beta\)估计趋向于\(0\)。
那么我们应该如何选取\(\lambda\)呢？接下来用代码实现一下。
ridge_model <- lm.ridge(Y~X1+X2+X3+X4+X5, data=food, lambda = seq(0,0.15,0.01))
# 画出图形，并作出lambdaGCV取最小值时的那条竖直线
matplot(ridge_model$lambda, t(ridge_model$coef),
xlab = expression(lambda), ylab = expression(beta),
type = "l", lty = 1, col=colors)
legend('bottomright',inset = 0.01, title = 'coefficients', c('beta_1','beta_2','beta_3','beta_4','beta_5'),lty = 1 ,col=colors)
abline(v = ridge_model$lambda[which.min(ridge_model$GCV)], col = "dark red", lty=3)
这种方法叫岭迹法，这个图叫做岭迹图。
岭迹图应该怎么看呢？
横坐标左边\(\lambda=0\)，则跟最小二乘解完全一样；
当\(\lambda\)略有增大，则各\(\beta\)系数取值迅速减小，即从不稳定趋于稳定。
图中类似喇叭形状的岭迹图，一般都存在多重共线性。
那\(\lambda\)应该如何选择呢？
一般通过观察，选取喇叭口附近（渐缩）的值，此时各\(\beta\)值已趋于稳定，但总的\(SSE\)又不是很大。
副产品：如果岭迹图中有回归系数始终趋于\(0\)值的变量，那么，这个变量就可以删除了。
但说实在的，从图中还是无法具体判断\(\lambda\)的具体值，只能确定在一个范围。
为了确定\(\lambda\)的取值，我们需要广义交叉验证的辅助。
# 下面的语句绘出lambda同GCV之间关系的图形
plot(ridge_model$lambda, ridge_model$GCV, type = "l", xlab = expression(lambda), ylab = "Generalized Cross-Validation: RMSE" )
abline(v = ridge_model$lambda[which.min(ridge_model$GCV)], col = "dark red", lty=3)
交叉验证就是通过把数据集分成若干等份，每一份轮流作测试集，其余数据作训练集进行轮流训练与测试。
广义交叉验证是在交叉验证的基础上，取多次验证结果的均值。
这里的交叉验证是为了找到最小的\(RMSE\)，从而定量地找到岭回归的最佳参数\(\lambda=0.05\)。
ridge_model## X1 X2 X3 X4 X5
## 0.00 -3.496563 0.1253296 0.07366672 2.677589 3.453448 -4.491117
## 0.01 -6.908865 0.1339418 0.06699854 2.561683 3.068801 -4.021525
## 0.02 -9.604300 0.1408163 0.06160782 2.466969 2.760539 -3.641925
## 0.03 -11.772263 0.1464115 0.05715854 2.387893 2.508439 -3.328599
## 0.04 -13.540958 0.1510376 0.05342317 2.320693 2.298816 -3.065490
## 0.05 -15.000261 0.1549120 0.05024216 2.262732 2.122083 -2.841342
## 0.06 -16.215045 0.1581917 0.04750030 2.212104 1.971320 -2.648031
## 0.07 -17.233312 0.1609930 0.04511229 2.167400 1.841413 -2.479545
## 0.08 -18.091338 0.1634037 0.04301366 2.127552 1.728499 -2.331339
## 0.09 -18.817051 0.1654914 0.04115473 2.091738 1.629605 -2.199917
## 0.10 -19.432298 0.1673090 0.03949662 2.059316 1.542408 -2.082540
## 0.11 -19.954415 0.1688985 0.03800845 2.029773 1.465068 -1.977038
## 0.12 -20.397337 0.1702937 0.03666540 2.002698 1.396105 -1.881664
## 0.13 -20.772386 0.1715220 0.03544725 1.977757 1.334320 -1.794999
## 0.14 -21.088863 0.1726061 0.03433741 1.954673 1.278727 -1.715878
## 0.15 -21.354473 0.1735644 0.03332211 1.933218 1.228512 -1.643335
当\(\lambda=0.05\)时， \[\begin{equation}
\beta=\left[
\begin{matrix}
- 15.000261\\
0.1549120\\
0.05024216\\
2.262732\\
2.1220827\\
- 2.841342
\end{matrix}
\right]
\end{equation}\]
让我们回忆一下之前的回归参数， \[\begin{equation}
\beta=\left[
\begin{matrix}
-3.49656\\
0.12533\\
0.07367\\
2.67759\\
3.45345\\
-4.49112
\end{matrix}
\right]
\end{equation}\]
可以发现，通过将\(\beta\)控制在一个合理范围内，缩小了\(\beta\)从而平衡了偏差和方差。
food[,Y_hat:=X %*% beta]
SSE <- sum((food$Y-food$Y_hat)^2)
SSR <- sum((food$Y_hat-mean(food$Y))^2)
SST <- SSE+SSR
R_sqr <- SSR/SST;R_sqr## [1] 0.9680447
此时，\(R^2=0.9680\)，预测图形为
ggplot(data = food, aes(x=Y_hat, y=Y)) +
geom_point() +
geom_abline(intercept=0, slope=1, col='dark red', lty=3)
岭回归总结
1.岭回归可以解决特征数量比样本量多的问题
2.岭回归作为一种缩减算法可以判断哪些特征重要或者不重要，有点类似于降维的效果
3.缩减算法可以看作是对一个模型增加偏差的同时减少方差
通过上一节课程，我们知道不只有岭回归，还有\(Lasso\)回归。
那么，同样作为解决过拟合和多重共线性的回归方式，\(Lasso\)和岭回归的区别是什么呢？
区别
惩罚项不同 \[\begin{equation} \beta^{lasso} = argmin\big[\sum(Y-\beta X)^2+\lambda\sum|\beta|\big] \end{equation}\]
对比一下岭回归， \[\begin{equation} \beta^{ridge} = argmin\big[\sum(Y-\beta X)^2+\lambda\sum\beta^2\big] \end{equation}\]
发现区别在哪儿了么？
这样的区别会导致哪里不同呢？
从几何的角度出发，
\(Lasso\)回归的惩罚函数映射到二维空间的话，就会形成角，一旦角与抛物面相交，就会导致\(\beta\)为\(0\)，这样\(\beta\)对应的变量就是一个可抛弃的变量。
所以\(Lasso\)更易求解，而且对应的很多系数能为0，也就起到了筛选变量的作用。
从数据上看
X <- as.matrix(food[,-c(1,7)])
Y <- as.matrix(food[,1])
lasso_model <- lars(X, Y, type='lasso')
\(Lasso\)的代码需要将\(X\)和\(Y\)分开，且矩阵化才能实现。
lasso_model##
## Call:
## lars(x = X, y = Y, type = "lasso")
## R-squared: 0.97
## Sequence of LASSO moves:
## X1 X3 X2 X5 X4
## Var 1 3 2 5 4
## Step 1 2 3 4 5
由此可见，\(Lasso\)的变量选择依次是\(X_1\)，\(X_3\)，\(X_2\)，\(X_5\)，\(X_4\)
接下来通过\(Lasso\)的图来验证一下
plot(lasso_model, lty=1, col=colors, xlim=c(0,1.5))
legend('bottomright',inset = 0.01, title = 'coefficients', c('beta_1','beta_2','beta_3','beta_4','beta_5'),lty = 1 ,col=colors)
从图中可以看出，\(X_1\)的参数首先下降到\(0\)，然后依次是\(X_3\)，\(X_2\)，\(X_5\)和\(X_4\)。
summary(lasso_model)## LARS/LASSO
## Call: lars(x = X, y = Y, type = "lasso")
## Df Rss Cp
## 0 1 8854.4 258.6543
## 1 2 4786.4 136.3068
## 2 3 1174.5 27.9007
## 3 4 398.4 6.1773
## 4 5 371.7 7.3617
## 5 6 261.7 6.0000
具体我们要保留哪几个变量呢？从上表的\(Cp\)可以看出，当step=4时，\(Cp\)值出现了第一次上升，所以我们将变量锁定在step=3时的那些变量，也就是\(X_1\)，\(X_3\)和\(X_2\)
接下来我们看一下，当step=3时的\(\beta\)值
lasso_model$beta[4,]## X1 X2 X3 X4 X5
## 0.205445969 0.009196321 1.447159599 0.000000000 0.000000000
当step=3时， \[\begin{equation}
\beta=\left[
\begin{matrix}
-33.4\\0.205445969\\0.009196321\\1.447159599\\0\\0
\end{matrix}
\right]
\end{equation}\]
让我们回忆一下之前的回归参数， \[\begin{equation}
\beta=\left[
\begin{matrix}
-3.49656\\
0.12533\\
0.07367\\
2.67759\\
3.45345\\
-4.49112
\end{matrix}
\right]
\end{equation}\]
可以发现，我们直接将两个变量的系数给了\(0\)值，从而达到缩小了\(\beta\)平衡了偏差和方差的作用。
X <- as.matrix(food[,-c(1,7)])
food[,Y_lasso:=predict(lasso_model, X)$fit[,4]]
SSE <- sum((food$Y-food$Y_lasso)^2)
SSR <- sum((food$Y_lasso-mean(food$Y))^2)
SST <- SSE+SSR
R_sqr <- SSR/SST;R_sqr## [1] 0.9543767
此时，\(R^2=0.9544\)，预测图形为
ggplot(data = food, aes(x=Y_lasso, y=Y)) +
geom_point() +
geom_abline(intercept=0, slope=1, col='dark red', lty=3)
\(Lasso\)回归总结
\(Lasso\)回归是在岭回归的基础上发展起来的，如果模型的特征非常多，需要压缩，那么\(Lasso\)回归是很好的选择。