岭回归

岭回归(英文名：ridge regression, Tikhonov regularization)是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于最小二乘法。

LASSO

LASSO是由1996年Robert Tibshirani首次提出，全称Least absolute shrinkage and selection operator。该方法是一种压缩估计。它通过构造一个惩罚函数得到一个较为精炼的模型，使得它压缩一些回归系数，即强制系数绝对值之和小于某个固定值；同时设定一些回归系数为零。因此保留了子集收缩的优点，是一种处理具有复共线性数据的有偏估计。

弹性网络

1987年由德宾（Durbin）和威尔萧（Willshaw）在「自然」月刊所提出弹性网络（Elastic Network）。这篇只有2页的文章（郝伯非尔德和坦克论文足足有11页），给人的第一印象是直觉，简单，而有力。弹性网络是一种社交联系。你所有的联系人都显示在缩略图列表中，排序依据是你与此人之间联系的强弱。每当Color侦测到你和其他用户在物理上接近时（即跟你在一起），这个应用程序就会调整你们两人的关系的强度。所以当你开启这个应用程序，并打开你的联系人名单时，你会看到亲密的朋友和家人排在最前面。

数据理解和准备

本文使用的数据在ElemStatLearn包里面，不过目前这个包在r里面不能直接下载了，这里为大家提供一种方法
下载方法

以ElemStatLearn包为例

在r studio4.0.3版本中输入install.packages(“ElemStatLearn”）无法下载ElemStatLearn包

提示CRAN下载地址为https://cran.r-project.org/web/packages/ElemStatLearn/index.html

选择最新版本下载

下载后在r studio中选择Packages——install——install from改为Packages Archive
File(.zip;.tar.gz)——Package archive选择刚刚下载的程序包——点击“Install”

安装完成即可library该程序包

library(ElemStatLearn)
prostate
1
2

此数据收集了97位男性的十个变量的数据，分别为：

• lcavol：肿瘤体积的对数值
• lweight：前列腺重量的对数值
• age：患者年龄
• lbph：良性前列腺增生（BPH）量的对数值，非癌症性质的前列腺增生
• svi：贮精囊是否侵入（1代表是，0代表否）
• lcp：包膜穿透值的对数值
• gleason：患者的Gleson评分，评分越高越危险
• pgg45：Gleson评分为4或5的百分比（高等级癌症）
• lpsa：PSA值的对数值，响应变量
• train：一个逻辑向量（TRUE or False，用来区分训练数据和测试数据）

一. 数据预处理

由于我们的gleson变量是一个特征变量，我们可以把其转化为指标变量，0代表评分为6,1表示评分为7或者更高，只需要这么一行代码即可

prostate$gleason <- ifelse(prostate$gleason==6,0,1)
1

相关系数矩阵：

library(corrplot)
corrplot(cor(prostate))
1
2

二.训练集和测试集的划分

train <- subset(prostate,train=T)[,1:9]
test <- subset(prostate,train=F)[,1:9]
1
2

prostate数据集的最后一列变量名为train，若train=true则为训练集否则测试集
我们只需要从中提取出训练集和测试集就可以了

三.模型构建与评价

1.最优子集

library(leaps)
a <- regsubsets(lpsa~.,data=train)
b <- summary(a) 
which.min(b$bic)
1
2
3
4

通过BIC信息准则来看，我们应该选择的3个最优子集

plot(a,scale = "bic",main="Best Sunset Features")
1

于是上图告诉我们具有最小BIC值的模型中的3个特征是：lcavol，lweight和svi
下面进行拟合：

ols <- lm(lpsa~lcavol+svi+lweight,data=train)
plot(ols$fitted.values,train$lpsa,xlab="Predicted"
     ,ylab="Acuual",main="Predicted vs Actual")
1
2
3

从图中可以看出，在训练集上线性拟合的很好，也不存在异方差。然后看看模型在测试集上的表现，使用predicte函数指定newdata=test，如下所示：

pred.test <- predict(ols,newdata=test) 
  plot(pred.test,test$lpsa,
       xlab="Predicted"
      ,ylab="Acuual",main="Predicted vs Actual")
1
2
3
4

计算残差的平方：

resid <- test$lpsa-pred.test
mean(resid^2)
1
2

MSE值为0.488以此为基准继续下面的内容

2.岭回归

在岭回归中，我们的模型会包括全部8个特征，使用岭回归的包为glmnet。岭回归需要的数据是矩阵而不是数据框。岭回归的命令格式为glment（x=矩阵，y=响应变量，famliy=分布函数，alpha=0）当alpha=0是表示使用岭回归，当alpha=1时表示使用LASSO回归
首先把数据转换为矩阵形式

x <- as.matrix(train[,1:8])
y <- train[,9]
ridge <- glmnet(x,y,famliy="gaussian",alpha=0) 
 print(ridge)
1
2
3
4

以第一行为列可以看出自由度df=7，即模型包含的特征的数量为7.请记住在岭回归中
，这个数量是不变的。还可以看出解释偏差（%DEV）为0.00，以及这一行的调优系数lambda为878，此处即可决定选择在那个测试集上使用那个lambda。我们首先来看一下，基本的统计图，使用label=T可以给曲线加上注释

plot(ridge,label=TRUE)
1

可以看到横坐标是L1范数，我们还可以将其更改为lambda或者dev

plot(ridge,label=TRUE,xvar = "dev")
plot(ridge,label=TRUE,xvar="lambda")
1
2

这张图非常有价值，它表明当lambda值减少时，压缩参数随之减少，而系数绝对值随之增大。要想看lambda为一个特定值是，可以使用coef命令。现在看一下lambda为0.1时，系数值是多少。指定参数lambda=0.1，如下所示：

ridge.coef <- coef(ridge,s=0.1)
ridge.coef
1
2

需要注意的是pgg45，lcp，age的系数值接近0，但还不是0。接下来看一下在测试集上的效果：

newx <- as.matrix(test[,1:8])
length(newx)
ridge.y <- predict(ridge,newx = newx,type="response",s=0.1)
plot(ridge.y,test$lpsa,
     xlab="Predicted"
    ,ylab="Acuual",main="Ridge Regression")
1
2
3
4
5
6

计算残差平方的平均值：

ridge.resid <- ridge.y-test$lpsa
mean(ridge.resid^2)
1
2

相比较于之前的0.48，这里的44稍微减少的一点，这时候看一下LASSO回归的效果

3.LASSO回归

lasso <- glmnet(x,y,family = "gaussian",alpha=1)
print(lasso)

1
2
3

可以看出第一列是自由度由7变为8时，lambda大约为0.023，因此我们应该使用这个lambda值

plot(lasso,label=T,xvar="lambda")
1

lasso.y <- predict(lasso,newx = newx,type="response",s=0.023)
lasso.resid <- lasso.y-test$lpsa
mean(lasso.resid^2)
1
2
3

输出结果为0.44和岭回归差不多

4.弹性网络

caret包旨在解决分类问题训练回归模型。首先我们需要找到集中于KaTeX parse error: Undefined control sequence: \弹 at position 9: \lambda和\̲弹̲性网络混合参数alpha的最优组合。可以通过下面3个简单的步骤完成：

• 使用R基础包里的expand.grid（）函数，建立一个向量存储我们要研究的$\alpha 和\lambda$的所有的组合
• 使用caret包中的trainControl（）函数确定重取样方法
• P在caret包中的train（）函数使用glmnet（）训练模型来选择$\alpha 和\lambda$
  可以按照下面两个规则来选择超参数
• $\alpha$从0到1，每次增加0.2；请记住$\alpha$被绑定在0到1之间
• $\lambda$从0.到0.2，使lasso回归和岭回归的$\lambda$位于这个$\lambda$之间
  具体操作如下：

grid <- expand.grid(.alpha=seq(0,1,0.2),.lambda=seq(0.,0.2,0.02))
table(grid)
1
2

对于重取样方法，我们要在代码中将method参数指定为LOOCV。然后通过train函数确定最优的弹性网络了。这个函数和lm函数很相似，只需在函数语法中加上method=“glmnet”,
trControl=control,tuneGrid=grid。将结果存储在enent.train的对象中

library(caret)
control <- trainControl(method="LOOCV")
enet.train <- train(lpsa~.,data=train,method="glmnet",
                 trControl=control,tuneGrid=grid)
enet.train
1
2
3
4
5

通过结果可以看出最优的$\alpha =0,\lambda =0.08$
验证模型效果：

enet <- glmnet(x,y,family = "gaussian",alpha=0,lambda=0.08)
coef(enet,s=0.08,exact=T)
enet.y <- predict(enet,newx=newx,type="response",s=0.08)
plot(enet.y,test$lpsa,
xlab="Predicted"
,ylab="Acuual",main="Elastic Net")
mean((enet.y-test$lpsa)^2)
1
2
3
4
5
6
7

可以看出误差大约为0.43相比较前面两个已经得到了改进，因此我们应该选择弹性网路的这种方法
因此得到回归方程：
$$\begin{gathered} lpsa=0.34+0.47lcavol+0.61lweight-0.01age+0.07lbph+ \\ 0.67svi-0.04lcp+0.29gleson+0.002pgg45 \end{gathered}$$