1、数据理解和数据准备

> library(pacman)
> p_load(car, corrplot, leaps, glmnet, caret, dplyr)
> load("./data_set/prostate.RData")
> str(prostate)

## 'data.frame':    97 obs. of  10 variables:
##  $ lcavol : num  -0.58 -0.994 -0.511 -1.204 0.751 ...
##  $ lweight: num  2.77 3.32 2.69 3.28 3.43 ...
##  $ age    : int  50 58 74 58 62 50 64 58 47 63 ...
##  $ lbph   : num  -1.39 -1.39 -1.39 -1.39 -1.39 ...
##  $ svi    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ lcp    : num  -1.39 -1.39 -1.39 -1.39 -1.39 ...
##  $ gleason: int  6 6 7 6 6 6 6 6 6 6 ...
##  $ pgg45  : int  0 0 20 0 0 0 0 0 0 0 ...
##  $ lpsa   : num  -0.431 -0.163 -0.163 -0.163 0.372 ...
##  $ train  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

 lcavol ：肿瘤体积的对数值
 lweight ：前列腺重量的对数值
 age ：患者年龄（以年计）
 lbph ： 良性前列腺增生（BPH）量的对数值，非癌症性质的前列腺增生
 svi ：贮精囊侵入，一个指标变量，表示癌细胞是否已经透过前列腺壁侵入贮精囊（ 1 =是，0 =否）
 lcp ：包膜穿透度的对数值，表示癌细胞扩散到前列腺包膜之外的程度
 gleason ：患者的Gleason评分；由病理学家进行活体检查后给出（2~10），表示癌细胞的变异程度——评分越高，程度越危险
 pgg45 ：Gleason评分为4或5所占的百分比（高等级癌症）
 lpsa ：PSA值的对数值，响应变量
 train ：一个逻辑向量（TRUE 或 FALSE，用来区分训练数据和测试数据）

通过散点图来理解数据：

> plot(prostate)

散点图

可以看出，结果变量lpsa和预测变量lcavol之间确实存在明显的线性关系。
lpsa为结果变量，我们先看看它的数据分布情况：

> table(prostate$gleason)

## 
##  6  7  8  9 
## 35 56  1  5

> ggplot2::ggplot(prostate, aes(x = as.factor(gleason), y = lpsa)) + 
+     stat_boxplot(outlier.color = "red") + 
+     theme_bw()

箱线图

只有1个观测的Gleason评分是8.0，只有5个观测的Gleason评分是9.0。看了上面的图，我认为最好的选择是，将这个特征转换为一个指标变量，0表示评分为6，1表示评分为7或更高，删除特征可能会损失模型的预测能力。缺失值也可能会在我们将要使用的glmnet包中引起问题。

> prostate$gleason <- ifelse(prostate$gleason < 7, 0, 1)
> table(prostate$gleason)

## 
##  0  1 
## 35 62

查看变量之间的相关性：

> prostate %>% cor(.) %>% corrplot.mixed(.)

变量之间的相关性

可以看到，首先，lpsa和肿瘤体积（lcavol）高度相关。其次，多重共线性是个问题。例如，肿瘤体积（lcavol）还与包膜穿透（lcp）相关，而包膜穿透还与贮精囊侵入（svi）相关。

拆分训练集和测试集。数据集中已经有一个特征指明训练集和非训练集，直接提取：

> train <- subset(prostate, train == TRUE) %>% select(-"train")
> test <- subset(prostate, train == FALSE) %>% select(-"train")
> dim(train)

## [1] 67  9

> dim(test)

## [1] 30  9

2、模型构建与评价

2.1 最优子集

建立一个最小子集对象：

> subfit <- regsubsets(lpsa ~ ., data = train)
> b.sum <- summary(subfit)
> # 使用贝叶斯信息准则确定最优子集
> which.min(b.sum$bic)

## [1] 3

结果告诉我们，3特征模型具有最小的BIC值。可以通过一个统计图查看模型性能和子集组合之间的关系：

> par(mfrow = c(1, 2))
> plot(b.sum$bic, type = "line", xlab = "# of Features", ylab = "BIC", main = "BIC score by Features Inclusion")
> plot(subfit, scale = "bic", main = "Best Subset Features")

统计图

上右图告诉我们具有最小BIC值的模型中的3个特征是：lcavol、lweight和gleason。建立线性模型：

> ols <- lm(lpsa ~ lcavol + lweight + gleason, data = train)
> # 拟合值与实际值的对比
> tibble(Predicted = ols$fitted.values, Actual = train$lpsa) %>% ggplot2::ggplot(aes(Predicted, 
+     Actual)) + geom_point() + geom_smooth() + theme_bw()

拟合值与实际值

从图中可以看出，在训练集上线性拟合表现得很好，也不存在异方差性。然后看看模型在测试集上的表现：

> pred.subfit <- predict(ols, newdata = test)
> tibble(Predicted = pred.subfit, Actual = test$lpsa) %>% ggplot2::ggplot(aes(Predicted, 
+     Actual)) + geom_point() + geom_smooth() + theme_bw()

测试集表现

这个图还不是太难看。总体来说，图中呈现一种线性关系，只不过当PSA值比较高时，有两个离群点。
计算均方误差：

> resid.subfit <- test$lpsa - pred.subfit
> mse <- mean(resid.subfit^2)
> print(mse)

## [1] 0.5084126

2.2 岭回归

glmnet程序包要求输入特征存储在矩阵中，而不是在数据框中。岭回归的命令形式为 glmnet(x=输入矩阵, y=响应变量, family=分布函数, alpha=0) 。这里的alpha为 0 时，表示进行岭回归；alpha为 1 时，表示进行LASSO回归。

> x <- as.matrix(train[, 1:8])
> y <- train[, 9]
> ridge <- glmnet(x, y, family = "gaussian", alpha = 0)
> 
> # 展示非0系数的数量，解释偏差百分比以及相应的λ值。程序包中算法默认的计算次数是100，
> # 但如果偏差百分比在两个λ值之间的提高不是很显著的话，算法会在100次计算之前停止。
> print(ridge)

## 
## Call:  glmnet(x = x, y = y, family = "gaussian", alpha = 0) 
## 
##     Df    %Dev Lambda
## 1    8 0.00000 878.90
## 2    8 0.00559 800.80
## 3    8 0.00613 729.70
## 4    8 0.00672 664.80
## 5    8 0.00737 605.80
##  ------------------------------
## 95   8 0.69230   0.14
## 96   8 0.69350   0.13
## 97   8 0.69460   0.12
## 98   8 0.69550   0.11
## 99   8 0.69640   0.10
## 100  8 0.69710   0.09

以第100行为例。可以看出非0系数，即模型中包含的特征的数量为8。请记住，在岭回归中，这个数是不变的。还可以看出解释偏差百分比为0.6971，以及这一行的调优系数λ的值为0.09。
此处即可决定选择在测试集上使用哪个λ。这个λ值应该是0.09，但是为了简单起见，在测试集上可以试一下0.10。

> # xvar='lambda'看系数值如何随着λ的变化而变化
> # xvar='dev'看系数值如何随解释偏差百分比变化 
> plot(ridge, xvar = "lambda", label = TRUE)

系数随lambda变化

指定参数 s=0.1 ，告诉glmnet在拟合模型时使用具体的λ值，而不是从λ值的两侧选值插入。如下所示：

> ridge.coef <- coef(ridge, s = 0.1)
> ridge.coef

## 9 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)  0.130475478
## lcavol       0.457279371
## lweight      0.645792042
## age         -0.017356156
## lbph         0.122497573
## svi          0.636779442
## lcp         -0.104712451
## gleason      0.346022979
## pgg45        0.004287179

需要特别注意的是， age、lcp 和 pgg45 的系数非常接近0，但还不是0。

> plot(ridge, xvar = "dev", label = TRUE)

系数值随解释偏差百分比变化

测试集上测试：

> newx <- as.matrix(test[, 1:8])
> ridge.y <- predict(ridge, newx = newx, type = "response", s = 0.1)
> tibble(Predicted = ridge.y, Actual = test$lpsa) %>% ggplot2::ggplot(aes(Predicted, 
+     Actual)) + geom_point() + geom_smooth() + labs(title = "Ridge Regression") + 
+     theme_bw()

测试集上的表现

岭回归中预测值和实际值关系的统计图看上去与最优子集的非常相似，同样地，在PSA测量结果比较大的一端有两个有趣的离群点。在实际情况下，我建议对离群点进行更深入的研究，搞清楚是它们真的与众不同，还是我们忽略了什么。这就是领域专家的用武之地。与MSE基准的比较可能告诉我们一些不同的事。先算出残差，然后算出残差平方的平均值：

> ridge.resid <- ridge.y - test$lpsa
> mse.ridge <- mean(ridge.resid^2)
> print(mse.ridge)

## [1] 0.4783559

岭回归给出的MSE稍好一点。

2.3 LASSO回归

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

## 
## Call:  glmnet(x = x, y = y, family = "gaussian", alpha = 1) 
## 
##    Df    %Dev  Lambda
## 1   0 0.00000 0.87890
## 2   1 0.09126 0.80080
## 3   1 0.16700 0.72970
## 4   1 0.22990 0.66480
## 5   1 0.28220 0.60580
## --------------------------------
## 7   1 0.36150 0.50290
## 8   1 0.39140 0.45820
##---------------------------------
## 32  7 0.67460 0.04914
## 33  7 0.67650 0.04477
## 34  8 0.67970 0.04079
## --------------------------------
## 65  8 0.70180 0.00228
## 66  8 0.70180 0.00208
## 67  8 0.70180 0.00189
## 68  8 0.70180 0.00172
## 69  8 0.70180 0.00157

模型构建过程在69步之后停止了，因为解释偏差不再随着λ值的增加而减小。还要注意，Df列现在也随着λ变化。初看上去，当λ值为0.00157时，所有8个特征都应该包括在模型中。然而，出于测试的目的，我们先用更少特征的模型进行测试，比如7特征模型。从结果行中可以看到，λ值大约为0.045时，模型从7个特征变为8个特征。因此，使用测试集评价模型时要使用这个λ值。

> plot(lasso, xvar = "lambda", label = TRUE)

系数随lambda变化

请注意标号为8、3和6的曲线的表现，这几条曲线分别对应特征 pgg45 、 age 和 lcp 。看上去 lcp 一直接近于0，直到作为最后一个特征被加入模型。

> # 查看7特征模型的系数值
> lasso.coef <- coef(lasso, s = 0.045)
> print(lasso.coef)

## 9 x 1 sparse Matrix of class "dgCMatrix"
##                         1
## (Intercept) -0.1305900670
## lcavol       0.4479592050
## lweight      0.5910476764
## age         -0.0073162861
## lbph         0.0974103575
## svi          0.4746790830
## lcp          .           
## gleason      0.2968768129
## pgg45        0.0009788059

LASSO算法在λ值为0.045时，将 lcp 的系数归零。
LASSO模型在测试集上的表现：

> lasso.y <- predict(lasso, newx = newx, type = "response", s = 0.045)
> tibble(Predicted = lasso.y, Actual = test$lpsa) %>% ggplot2::ggplot(aes(Predicted, 
+     Actual)) + geom_point() + geom_smooth() + labs(title = "LASSO") + theme_bw()

LASSO模型测试集表现

计算MSE：

> lasso.resid <- lasso.y - test$lpsa
> mse.lasso <- mean(lasso.resid^2)
> tibble(mse = mse, mse.ridge = mse.ridge, mse.lasso = mse.lasso) %>% print()

## # A tibble: 1 x 3
##     mse mse.ridge mse.lasso
##   <dbl>     <dbl>     <dbl>
## 1 0.508     0.478     0.444

看起来我们的统计图和前面一样，只是MSE值有了一点点改进。

2.4 弹性网络

回忆一下，α = 0表示岭回归惩罚，α = 1表示LASSO惩罚，弹性网络参数为0≤α≤1。弹性网络要做的调整是不但要解出λ值，还要解出弹性网络参数α。

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

##       .lambda
## .alpha 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
##    0   1    1    1    1    1   1    1    1    1    1   1
##    0.2 1    1    1    1    1   1    1    1    1    1   1
##    0.4 1    1    1    1    1   1    1    1    1    1   1
##    0.6 1    1    1    1    1   1    1    1    1    1   1
##    0.8 1    1    1    1    1   1    1    1    1    1   1
##    1   1    1    1    1    1   1    1    1    1    1   1

可以确认这就是我们想要的结果——α值在0和1之间，λ值在0和0.2之间。

> # LOOCV重取样方法
> control <- trainControl(method = "LOOCV")
> enet.train <- train(lpsa ~ ., data = train, method = "glmnet", trControl = control, 
+     tuneGrid = grid)
> enet.train

## glmnet 
## 
## 67 samples
##  8 predictor
## 
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 66, 66, 66, 66, 66, 66, ... 
## Resampling results across tuning parameters:
## 
##   alpha  lambda  RMSE       Rsquared   MAE      
##   0.0    0.00    0.7498311  0.6091434  0.5684548
##   0.0    0.02    0.7498311  0.6091434  0.5684548
##   0.0    0.04    0.7498311  0.6091434  0.5684548
##   0.0    0.06    0.7498311  0.6091434  0.5684548
##   0.0    0.08    0.7498311  0.6091434  0.5684548
##   1.0    0.08    0.7819175  0.5755415  0.6082954
##   ----------------------------------------------------------------
##   1.0    0.12    0.7921572  0.5692525  0.6146159
##   1.0    0.14    0.7999326  0.5642789  0.6198758
##   1.0    0.16    0.8089248  0.5583637  0.6265797
##   1.0    0.18    0.8185327  0.5521348  0.6343174
##   1.0    0.20    0.8259445  0.5494268  0.6411104
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0 and lambda = 0.08.

我们选择最优模型的原则是RMSE值最小，模型最后选定的最优参数组合是α = 0，λ = 0.08。
在测试集上验证模型：

> enet <- glmnet(x, y, family = "gaussian", alpha = 0, lambda = 0.08)
> enet.coef <- coef(enet, s = 0.08)
> print(enet.coef)

## 9 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)  0.137811097
## lcavol       0.470960525
## lweight      0.652088157
## age         -0.018257308
## lbph         0.123608113
## svi          0.648209192
## lcp         -0.118214386
## gleason      0.345480799
## pgg45        0.004478267

> enet.y <- predict(enet, newx = newx, type = "response", s = 0.08)
> tibble(Predicted = enet.y, Actual = test$lpsa) %>% ggplot2::ggplot(aes(Predicted, 
+     Actual)) + geom_point() + geom_smooth() + labs(title = "LASSO") + theme_bw()

弹性网络模型测试集表现

计算MSE：

> enet.resid <- enet.y - test$lpsa
> mse.enet <- mean(enet.resid^2)
> tibble(mse = mse, mse.ridge = mse.ridge, mse.lasso = mse.lasso, mse.enet = mse.enet) %>% 
+     print()

## # A tibble: 1 x 4
##     mse mse.ridge mse.lasso mse.enet
##   <dbl>     <dbl>     <dbl>    <dbl>
## 1 0.508     0.478     0.444    0.480

这个模型的误差与岭回归很接近。在测试集上，LASSO模型在误差方面表现最好。模型可能过拟合了！

2.5 使用glmnet进行交叉验证

在K折交叉验证中，数据被划分成k个相同的子集（折），每次使用k-1个子集拟合模型，然后使用剩下的那个子集做测试集，最后将k次拟合的结果综合起来（一般取平均数），确定最后的参数。

> set.seed(111)
> # 数据量比较少，使用3折交叉验证
> lasso.cv <- cv.glmnet(x, y, nfolds = 3)
> plot(lasso.cv)

交叉验证

图中两条垂直的虚线表示取得MSE最小值的logλ（左侧虚线）和距离最小值一个标准误差的logλ。
如果有过拟合问题，那么距离最小值一个标准误的位置是非常好的解决问题的起点。你还可以得到这两个λ的具体值:

> # 最小值
> lasso.cv$lambda.min

## [1] 0.06495561

> # 一个标准误距离
> lasso.cv$lambda.1se

## [1] 0.2877937

使用lambda.1se查看系数并在测试集上进行模型验证：

> coef(lasso.cv, s = "lambda.1se")

## 9 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept) 0.883375319
## lcavol      0.425544199
## lweight     0.270837436
## age         .          
## lbph        .          
## svi         0.105666285
## lcp         .          
## gleason     0.006825567
## pgg45       .

> lasso.y.cv <- predict(lasso.cv, newx = newx, type = "response", s = "lambda.1se")
> lasso.cv.resid <- lasso.y.cv - test$lpsa
> mse.lasso.cv <- mean(lasso.cv.resid^2)
> print(mse.lasso.cv)

## [1] 0.5292918

这个模型的误差为0.5292918，只有4个特征，排除了age、lcp、lbph和pgg45。看上去，lambda.1se还不是最优的选择。我们看看使用lambda.min选择的模型是否可以再次改善样本预测结果：

> lasso.y.cv.min <- predict(lasso.cv, newx = newx, type = "response", s = "lambda.min")
> lasso.min.resid <- lasso.y.cv.min - test$lpsa
> mse.lasso.min <- mean(lasso.min.resid^2)
> print(mse.lasso.min)

## [1] 0.4468133

3、模型选择

> tibble(mse = mse, mse.ridge = mse.ridge, mse.lasso = mse.lasso, mse.enet = mse.enet, 
+     mse.lasso.cv = mse.lasso.cv, mse.lasso.min = mse.lasso.min) %>% print()

## # A tibble: 1 x 6
##     mse mse.ridge mse.lasso mse.enet mse.lasso.cv mse.lasso.min
##   <dbl>     <dbl>     <dbl>    <dbl>        <dbl>         <dbl>
## 1 0.508     0.478     0.444    0.480        0.529         0.447

仅看误差的话，7特征LASSO模型表现最好。
在本例的样本规模之下，仅改变随机数种子或重新划分训练集和测试集都可能使结果发生大的改变。到头来，这些结果非但不能提供答案，还可能引起更多问题。所以我们需要更多数据。
同时，正则化是一项强大的技术，与其他建模技术相比，既可以提高计算效率，还可以提取更有意义的特征。