对连续性变量进行LASSO回归分析 cox

对连续性变量进行LASSO回归

jmzeng@163.com

6/19/2017

• 我的博客
• 我们的论坛
• 捐赠我

　安装并加载必须的packages

如果你还没有安装，就运行下面的代码安装：

install.packages('lars')
install.packages('glmnet')

如果你安装好了，就直接加载它们即可

library(lars)
# https://cran.r-project.org/web/packages/lars/lars.pdf
library(glmnet)
data(diabetes)

首先查看测试数据的基本信息

这是一个糖尿病相关的数据集，包含两个x矩阵和一个y结果向量，分别是：

• x a matrix with 10 columns (自变量)
• y a numeric vector (因变量)
• x2 a matrix with 64 columns

attach(diabetes)
summary(x)

##       age                 sex                bmi           
##  Min.   :-0.107226   Min.   :-0.04464   Min.   :-0.090275  
##  1st Qu.:-0.037299   1st Qu.:-0.04464   1st Qu.:-0.034229  
##  Median : 0.005383   Median :-0.04464   Median :-0.007284  
##  Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.000000  
##  3rd Qu.: 0.038076   3rd Qu.: 0.05068   3rd Qu.: 0.031248  
##  Max.   : 0.110727   Max.   : 0.05068   Max.   : 0.170555  
##       map                  tc                 ldl           
##  Min.   :-0.112400   Min.   :-0.126781   Min.   :-0.115613  
##  1st Qu.:-0.036656   1st Qu.:-0.034248   1st Qu.:-0.030358  
##  Median :-0.005671   Median :-0.004321   Median :-0.003819  
##  Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.000000  
##  3rd Qu.: 0.035644   3rd Qu.: 0.028358   3rd Qu.: 0.029844  
##  Max.   : 0.132044   Max.   : 0.153914   Max.   : 0.198788  
##       hdl                 tch                 ltg           
##  Min.   :-0.102307   Min.   :-0.076395   Min.   :-0.126097  
##  1st Qu.:-0.035117   1st Qu.:-0.039493   1st Qu.:-0.033249  
##  Median :-0.006584   Median :-0.002592   Median :-0.001948  
##  Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.000000  
##  3rd Qu.: 0.029312   3rd Qu.: 0.034309   3rd Qu.: 0.032433  
##  Max.   : 0.181179   Max.   : 0.185234   Max.   : 0.133599  
##       glu           
##  Min.   :-0.137767  
##  1st Qu.:-0.033179  
##  Median :-0.001078  
##  Mean   : 0.000000  
##  3rd Qu.: 0.027917  
##  Max.   : 0.135612

summary(y)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    25.0    87.0   140.5   152.1   211.5   346.0

boxplot(y)

#summary(x2)

其中x矩阵含有10个变量，分别是：“age” “sex” “bmi” “map” “tc” “ldl” “hdl” “tch” “ltg” “glu” 它们都在一定程度上或多或少的会影响个体糖尿病状态。

数据的详细介绍见 Efron, Hastie, Johnstone and Tibshirani (2003) “Least Angle Regression” (with discussion) Annals of Statistics;

先看看因变量和每一个自变量的关系

oldpar=par()
par(mfrow=c(2,5))
for(i in 1:10){
  plot(x[,i], y)
  abline(lm(y~x[,i]))
}

par(oldpar)

再对数据进行简单的线性回归

model_ols <- lm(y ~ x)
summary(model_ols)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -155.829  -38.534   -0.227   37.806  151.355 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  152.133      2.576  59.061  < 2e-16 ***
## xage         -10.012     59.749  -0.168 0.867000    
## xsex        -239.819     61.222  -3.917 0.000104 ***
## xbmi         519.840     66.534   7.813 4.30e-14 ***
## xmap         324.390     65.422   4.958 1.02e-06 ***
## xtc         -792.184    416.684  -1.901 0.057947 .  
## xldl         476.746    339.035   1.406 0.160389    
## xhdl         101.045    212.533   0.475 0.634721    
## xtch         177.064    161.476   1.097 0.273456    
## xltg         751.279    171.902   4.370 1.56e-05 ***
## xglu          67.625     65.984   1.025 0.305998    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.15 on 431 degrees of freedom
## Multiple R-squared:  0.5177, Adjusted R-squared:  0.5066 
## F-statistic: 46.27 on 10 and 431 DF,  p-value: < 2.2e-16

可以看到,如果用普通的线性回归分析,有4个变量(bmi,sex,map,ltg),还有截距都显著影响着糖尿病.

接下来进行LASSO回归

很简单，一个命令即可。目前最好用的拟合广义线性模型的R package是glmnet，由LASSO回归的发明人，斯坦福统计学家Trevor Hastie领衔开发。它的特点是对一系列不同λ值进行拟合，每次拟合都用到上一个λ值拟合的结果，从而大大提高了运算效率。

model_lasso <- glmnet(x, y, family="gaussian", nlambda=50, alpha=1)

参数family规定了回归模型的类型:

• family=“gaussian” 适用于一维连续因变量(univariate)
• family=“mgaussian” 适用于多维连续因变量(multivariate)
• family=“poisson” 适用于非负次数因变量(count)
• family=“binomial” 适用于二元离散因变量(binary)
• family=“multinomial” 适用于多元离散因变量(category)

参数nlambda=50让算法自动挑选50个不同的λ值，拟合出50个系数不同的模型。

参数alpha=1输入α值，1是它的默认值。参数α来控制应对高相关性(highly correlated)数据时模型的性状.

LASSO回归的结果解析

print(model_lasso)

## 
## Call:  glmnet(x = x, y = y, family = "gaussian", alpha = 1, nlambda = 50) 
## 
##       Df   %Dev    Lambda
##  [1,]  0 0.0000 45.160000
##  [2,]  2 0.1329 37.420000
##  [3,]  2 0.2352 31.010000
##  [4,]  2 0.3055 25.700000
##  [5,]  3 0.3542 21.290000
##  [6,]  3 0.3937 17.640000
##  [7,]  4 0.4213 14.620000
##  [8,]  4 0.4433 12.120000
##  [9,]  4 0.4584 10.040000
## [10,]  4 0.4688  8.319000
## [11,]  4 0.4759  6.893000
## [12,]  5 0.4833  5.712000
## [13,]  5 0.4913  4.733000
## [14,]  6 0.4968  3.922000
## [15,]  7 0.5009  3.250000
## [16,]  7 0.5052  2.693000
## [17,]  7 0.5081  2.232000
## [18,]  7 0.5101  1.849000
## [19,]  7 0.5115  1.532000
## [20,]  7 0.5125  1.270000
## [21,]  7 0.5131  1.052000
## [22,]  8 0.5137  0.871900
## [23,]  8 0.5143  0.722500
## [24,]  8 0.5147  0.598700
## [25,]  8 0.5149  0.496100
## [26,]  8 0.5151  0.411100
## [27,]  8 0.5152  0.340600
## [28,]  8 0.5153  0.282300
## [29,] 10 0.5156  0.233900
## [30,] 10 0.5163  0.193800
## [31,] 10 0.5167  0.160600
## [32,] 10 0.5170  0.133100
## [33,] 10 0.5172  0.110300
## [34,] 10 0.5174  0.091390
## [35,]  9 0.5174  0.075730
## [36,]  9 0.5174  0.062750
## [37,] 10 0.5175  0.052000
## [38,] 10 0.5175  0.043090
## [39,] 10 0.5176  0.035700
## [40,] 10 0.5176  0.029590
## [41,] 10 0.5177  0.024520
## [42,] 10 0.5177  0.020310
## [43,] 10 0.5177  0.016830
## [44,] 10 0.5177  0.013950
## [45,] 10 0.5177  0.011560
## [46,] 10 0.5177  0.009578
## [47,] 10 0.5177  0.007937

每一行代表了一个模型。

列Df是自由度，代表了非零的线性模型拟合系数的个数。

列%Dev代表了由模型解释的残差的比例，对于线性模型来说就是模型拟合的R^2(R-squred)。它在0和1之间，越接近1说明模型的表现越好，如果是0，说明模型的预测结果还不如直接把因变量的均值作为预测值来的有效。

列Lambda当然就是每个模型对应的λ值。

我们可以看到，随着λ的变小，越来越多的自变量被模型接纳进来，%Dev也越来越大。

第29行时，模型包含了所有20个自变量，%Dev也在0.5以上。

其实我们本应该得到50个不同的模型，但是连续几个%Dev变化很小时glmnet()会自动停止。

分析模型输出我们可以看到当Df大于9的时候，%Dev就达到了0.5，而且继续缩小λ，即增加更多的自变量到模型中，也不能显著提高%Dev。

所以我们可以认为当λ接近0.23时，得到的包含10个自变量的模型，已经是最好的描述这组数据的模型了。

我们也可以通过指定λ值，抓取出某一个模型的系数:

coef(model_lasso, s=c(model_lasso$lambda[29],0.23))

## 11 x 2 sparse Matrix of class "dgCMatrix"
##                        1            2
## (Intercept)  152.1334842  152.1334842
## age           -0.1744351   -0.3396002
## sex         -227.4412952 -227.6460806
## bmi          526.5411179  526.4325058
## map          315.1207791  315.2735154
## tc          -227.5703945 -236.5240061
## ldl           26.5812439   33.7330891
## hdl         -139.7688154 -135.9670209
## tch          109.4083069  110.4666063
## ltg          541.8941364  545.2162329
## glu           64.7035129   64.7539943

作图观察LASSO回归模型的系数是如何变化的

plot.glmnet(model_lasso, xvar = "norm", label = TRUE)

plot(model_lasso, xvar="lambda", label=TRUE)

图中的每一条曲线代表了每一个自变量系数的变化轨迹，纵坐标是系数的值，下横坐标是loga(λ)，上横坐标是此时模型中非零系数的个数。我们可以看到，黑线代表的自变量1在λ值很大时就有非零的系数，然后随着λ值变小不断变大。我们还可以尝试用xvar=“norm”和xvar=“dev”切换下横坐标。

用自带glmnet()函数模型进行预测

pre <- predict(model_lasso, newx=x , s=c(model_lasso$lambda[29],0.23))
head(cbind(y,pre))

##        y         1         2
## [1,] 151 204.54545 204.56986
## [2,]  75  70.58607  70.54513
## [3,] 141 175.79128 175.80730
## [4,] 206 163.01751 163.08081
## [5,] 135 127.76397 127.77541
## [6,]  97 105.66422 105.67519

用交叉验证来确定最优lambda值

公式中的lambda是重要的设置参数，它控制了惩罚的严厉程度，如果设置得过大，那么最后的模型参数均将趋于0，形成拟合不足。如果设置得过小，又会形成拟合过度。所以lambda的取值一般需要通过交叉检验来确定。这时候需要用lars包的glmnet()函数啦！

cv_fit <- cv.glmnet(x=x, y=y, alpha = 1, nlambda = 1000)

参数nlambda=1000让算法自动挑选1000个不同的λ值，拟合出1000个系数不同的模型。

其实cv.glmnet()函数的参数值得注意的还有两个，其中type.measure是用来指定交叉验证选取模型时希望最小化的目标参量。

还可以用nfolds指定fold数，或者用foldid指定每个fold的内容。

可视化交叉验证结果

plot.cv.glmnet(cv_fit)

因为交叉验证，对于每一个λ值，在红点所示目标参量的均值左右，我们可以得到一个目标参量的置信区间。两条虚线分别指示了两个特殊的λ值:

c(cv_fit$lambda.min,cv_fit$lambda.1se)

## [1] 1.140661 8.913526

lambda.min是指在所有的λ值中，得到最小目标参量均值的那一个。

而lambda.1se是指在lambda.min一个方差范围内得到最简单模型的那一个λ值。

因为λ值到达一定大小之后，继续增加模型自变量个数即缩小λ值，并不能很显著的提高模型性能，lambda.1se给出的就是一个具备优良性能但是自变量个数最少的模型。

用cv.glmnet()函数模型进行预测

pre <-  predict(cv_fit, newx=x[1:5,], type="response", s="lambda.1se") 
head(cbind(y,pre))

## Warning in cbind(y, pre): number of rows of result is not a multiple of
## vector length (arg 1)

##        y         1
## [1,] 151 197.19236
## [2,]  75  88.54846
## [3,] 141 176.50576
## [4,] 206 153.14809
## [5,] 135 123.88797

提取模型的系数

cv_fit$lambda.min

## [1] 1.140661

fit <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
fit$beta

## 10 x 1 sparse Matrix of class "dgCMatrix"
##             s0
## age    .      
## sex -190.88636
## bmi  521.26169
## map  293.55728
## tc   -95.35910
## ldl    .      
## hdl -221.79409
## tch    .      
## ltg  509.63138
## glu   51.29918

coef(cv_fit$glmnet.fit, s =  min(cv_fit$lambda))

## 11 x 1 sparse Matrix of class "dgCMatrix"
##                      1
## (Intercept)  152.13348
## age            .      
## sex         -222.99096
## bmi          526.27917
## map          312.40981
## tc          -183.52085
## ldl            .      
## hdl         -162.22910
## tch           92.59096
## ltg          528.00455
## glu           63.40113

cv_fit$lambda.1se

## [1] 8.913526

fit <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)
fit$beta

## 10 x 1 sparse Matrix of class "dgCMatrix"
##           s0
## age   .     
## sex   .     
## bmi 483.8524
## map 156.7522
## tc    .     
## ldl   .     
## hdl -78.9380
## tch   .     
## ltg 419.7262
## glu   .

coef(cv_fit$glmnet.fit, s = cv_fit$lambda.1se)

## 11 x 1 sparse Matrix of class "dgCMatrix"
##                     1
## (Intercept) 152.13348
## age           .      
## sex           .      
## bmi         483.78833
## map         156.74421
## tc            .      
## ldl           .      
## hdl         -78.93562
## tch           .      
## ltg         419.75889
## glu           .

因为每个fold间的计算是独立的。 我们还可以考虑运用并行计算来提高运算效率，使用parallel=TRUE可以开启这个功能。但是我们需要先装载package doParallel。

可以用并行加快计算速度

并行的示例代码是：

library(doParallel)
# Windows System
cl <- makeCluster(6)
registerDoParallel(cl)
cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class", parallel=TRUE)
stopCluster(cl)
# Linux System
registerDoParallel(cores=8)
cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class", parallel=TRUE)
stopImplicitCluster()

用model.matrix函数把自变量的离散变量变成矩阵输入

LASSO回归的特点是在拟合广义线性模型的同时进行变量筛选(Variable Selection)和复杂度调整(Regularization)。因此，不论目标因变量(dependent/response varaible)是连续的(continuous)，还是二元或者多元离散的(discrete)， 都可以用LASSO回归建模然后预测。

但是，glmnet只能接受数值矩阵作为模型输入，如果自变量中有离散变量的话，需要把这一列离散变量转化为几列只含有0和1的向量，这个过程叫做One Hot Encoding。

下面是一个例子

df=data.frame(Factor=factor(1:5), Character=c("a","a","b","b","c"),
              Logical=c(T,F,T,T,T), Numeric=c(2.1,2.3,2.5,4.1,1.1))
model.matrix(~., df)

##   (Intercept) Factor2 Factor3 Factor4 Factor5 Characterb Characterc
## 1           1       0       0       0       0          0          0
## 2           1       1       0       0       0          0          0
## 3           1       0       1       0       0          1          0
## 4           1       0       0       1       0          1          0
## 5           1       0       0       0       1          0          1
##   LogicalTRUE Numeric
## 1           1     2.1
## 2           0     2.3
## 3           1     2.5
## 4           1     4.1
## 5           1     1.1
## attr(,"assign")
## [1] 0 1 1 1 1 2 2 3 4
## attr(,"contrasts")
## attr(,"contrasts")$Factor
## [1] "contr.treatment"
## 
## attr(,"contrasts")$Character
## [1] "contr.treatment"
## 
## attr(,"contrasts")$Logical
## [1] "contr.treatment"

除此之外，如果我们想让模型的变量系数都在同一个数量级上，就需要在拟合前对数据的每一列进行标准化(standardize)，即对每个列元素减去这一列的均值然后除以这一列的标准差。这一过程可以通过在glmnet()函数中添加参数standardize=TRUE来实现。

模型其它技巧

指定模型系数的上限与下限

使用upper.limits和lower.limits,上限与下限可以是一个值，也可以是一个向量，向量的每一个值作为对应自变量的参数上下限。

调整惩罚参数

有时，在建模之前我们就想凸显某几个自变量的作用，此时我们可以调整惩罚参数,即设置penalty.factor。每个自变量的默认惩罚参数是1，把其中的某几个量设为0将使得相应的自变量不遭受任何惩罚