R语言自适应LASSO 多项式回归、二元逻辑回归和岭回归应用分析(上)

as_data_frame(y)

## # A tibble: 100 x 1
##            V1
##         <dbl>
##  1 -1.2748860
##  2  1.8434251
##  3  0.4592363
##  4  0.5640407
##  5  1.8729633
##  6  0.5275317
##  7  2.4346589
##  8 -0.8945961
##  9 -0.2059384
## 10  3.1101188
## # ... with 90 more rows

## 执行岭回归
glmnet(x , y 
                 ## “alpha=1”是套索惩罚, “alpha=0”是岭惩罚。
                 alpha = 0)

## 用10折CV进行岭回归
cv.glmnet(
                       ## 类型.测量：用于交叉验证的丢失。
                       type.measure = "mse",
                       ## K = 10 是默认值。
                       nfold = 10,
                       ##“alpha=1”是套索惩罚，“alpha=0”是岭惩罚。
                       alpha = 0)
## 惩罚vs CV MSE图

## 在误差最小λ处提取系数
cv$lambda.min

## [1] 0.1789759

## s：需要进行预测的惩罚参数“lambda”的值。默认值是用于创建模型的整个序列。
coef( s = lambda.min)

## 21 x 1 sparse Matrix of class "dgCMatrix"
##                        1
## (Intercept)  0.149041059
## V1           1.302684272
## V2           0.035835380
## V3           0.719936146
## V4           0.036473087
## V5          -0.863490158
## V6           0.605750873
## V7           0.123446432
## V8           0.376890626
## V9          -0.040012847
## V10          0.105999328
## V11          0.240967604
## V12         -0.066363634
## V13         -0.042048935
## V14         -1.092107794
## V15         -0.119566353
## V16         -0.035859663
## V17         -0.038827463
## V18          0.061785988
## V19         -0.001409608
## V20         -1.079879797

## 截距估计应该剔除。
(coef(cv, s = lambda.min))[-1]

vif <- function(x, lambda) {
    ZtZ <- cor(x)
    diag(solve(ZtZ + lambdaI  %*% ZtZ %*% solve(ZtZ + lambdaI)
 
##
    ggplot(mapping = aes(x = lambda, y = value, group = key, color = key)) +
    geom_line() +

## 执行自适应LASSO
glmnet(x =  y =
                  ## 类型。度量:用于交叉验证的损失。
                  ##“alpha=1”是套索惩罚，“alpha=0”是岭惩罚。
                  alpha = 1,
                  ##
                  ## 惩罚系数：可以对每个系数应用单独的惩罚因子。这是一个乘以“lambda”以允许差异收缩的数字。对于某些变量可以是0， 这意味着没有收缩，而且这个变量总是包含在模型中。对于所有变量，默认值为1（对于“exclude”中列出的变量，默认值为无限大）。注意：惩罚因子在内部被重新调整为与nvars相加，lambda序列将反映这种变化。

## 使用10折CV执行自适应套索
                        ## 类型。度量:用于交叉验证的损失。
类型。测量= " mse ",
                        ## K = 10 是默认值。
                        nfold = 10,
                        ## ‘alpha = 1’ 是套索惩罚，'alpha=0'是岭惩罚。
                        ##
                        ## 惩罚系数：可以对每个系数应用单独的惩罚因子。这是一个乘以“lambda”以允许差异收缩的数字。对于某些变量可以为0，这意味着没有收缩，并且该变量始终包含在模型中。对于所有变量，默认值为1（对于“exclude”中列出的变量，默认值为无限大）。注意：惩罚因子在内部被重新调整为与nvars相加，lambda序列将反映这种变化。
## 惩罚vs CV MSE图

## 在误差最小λ处提取系数
lambda.min

## [1] 0.7193664

## s：需要进行预测的惩罚参数“lambda”的值。默认值是用于创建模型的整个序列。
best_alasso_coef1

## 21 x 1 sparse Matrix of class "dgCMatrix"
##                      1
## (Intercept)  0.1269539
## V1           1.3863728
## V2           .        
## V3           0.7573538
## V4           .        
## V5          -0.8937983
## V6           0.5718800
## V7           .        
## V8           0.3654255
## V9           .        
## V10          .        
## V11          0.1824140
## V12          .        
## V13          .        
## V14         -1.1150736
## V15          .        
## V16          .        
## V17          .        
## V18          .        
## V19          .        
## V20         -1.1268794

##  R^2函数
## https://en.wikipedia.org/wiki/Coefficient_of_determination
    ##  总SS
    ss_tot <- sum((y - ybar)^2)
    ## 剩余 SS
    ss_res <- sum((y - yhat)^2)
    ## R^2 = 1 - ss_res/ ss_tot
## 调整R^2函数
## n个样本，p个参数
## 获取 R^2
 r_sq(as.vector(y_cont), as.vector(predict(alasso1, newx =

## [1] 0.906806

##获得调整R ^ 2
adj_r_sq(r_squared_alasso1, n = nrow(y_cont),

## [1] 0.9007934

## 交叉验证测试集R^2
## alasso1_cv$cvm[1] 是截距模型的交叉验证测试集均方误差。
1 - cvm[lambda == lambda.min] / cvm[1]

## [1] 0.8854662

lapply(unique(  foldid), function(id) {
    ## 拟合排除测试集 (foldid == id)
 glmnet(x = x_cont[alasso1_cv$foldid != id,],
                  y = y_cont[alasso1_cv$foldid != id], 
    ## 使用模型拟合最佳lambda测试集Yïhat
 predict(fit, newx = x_cont[alasso1_cv$foldid == id,],  
    ## 测试组 R^2
    1 - sum((y - y_pred)^2) / sum((y - mean(y))^2)
}) %>%

##  [1] 0.8197796 0.9090972 0.9499495 0.8019303 0.8637534 0.7184797 0.8579943 0.9250376 0.8300891
## [10] 0.9188004

## [1] 0.8594911

## # A tibble: 500 x 30
##            V1         V2          V3         V4         V5         V6          V7         V8
##         <dbl>      <dbl>       <dbl>      <dbl>      <dbl>      <dbl>       <dbl>      <dbl>
##  1  0.8212500  1.2155090 -0.64860899 -0.7001262 -1.9640742  1.1692107  0.28598652 -0.1664266
##  2  0.9264925 -1.1855031 -1.18297879  0.9828354  1.0693610 -0.2302219  0.57772625 -0.8738714
##  3 -1.5719712  0.8568961 -0.02208733  1.7445962 -0.4148403 -2.0289054 -1.31228181 -1.2441528
##  4  0.7419447 -0.9452052 -1.61821790  1.0015587 -0.4589488  0.5154490  0.29189973  0.1114092
##  5 -0.1333660  0.5085678  0.04739909 -0.4486953 -0.2616950 -0.1554108 -1.24834832 -1.0498054
##  6 -0.5672062  0.6020396 -2.10300909  0.3119233  0.3272173 -0.8671885  0.97512759 -0.7216256
##  7  1.9683411  2.5162198  1.61109738  1.0047913 -0.5194647  1.0738680 -0.16176095 -0.4267418
##  8  0.2857727 -1.7017703  1.41062569 -0.5823727 -1.3330908  1.7929250  0.06396841 -0.6818909
##  9 -0.5339434  0.1725089  0.93504676 -1.9956942 -0.9021089 -0.2624043  0.97406411  0.5166823
## 10  0.8081052 -0.9662501  0.54666915 -0.8388913  0.9665053  1.4039598  0.63502500  0.3429640
## # ... with 490 more rows, and 22 more variables: V9 <dbl>, V10 <dbl>, V11 <dbl>, V12 <dbl>,
## #   V13 <dbl>, V14 <dbl>, V15 <dbl>, V16 <dbl>, V17 <dbl>, V18 <dbl>, V19 <dbl>, V20 <dbl>,
## #   V21 <dbl>, V22 <dbl>, V23 <dbl>, V24 <dbl>, V25 <dbl>, V26 <dbl>, V27 <dbl>, V28 <dbl>,
## #   V29 <dbl>, V30 <dbl>

as_data_frame(y)

## # A tibble: 500 x 1
##    value
##    <dbl>
##  1     3
##  2     2
##  3     2
##  4     2
##  5     3
##  6     3
##  7     3
##  8     1
##  9     1
## 10     1
## # ... with 490 more rows

plot(ridge2, xvar = "lambda")