广义线性回归模型估计：所有线性回归的大仓库
The following article comes from 数据挖掘与R语言 Author 2万计量学者
1 引言
说起glmnet包就厉害了，它是斯坦福大学的几位统计大牛写的包，解决了带惩罚项的广义线性模型的参数估计问题。具体的模型涉及很广，包括线性回归、逻辑回归、泊松计数模型、cox回归模型、多分类逻辑回归以及多响应线性回归等等。
1.1 问题描述
在介绍glmnet的用法之前，我们首先用数学语言来描述glmnet包解决了哪类问题，如下：
该表达式前半部分即损失函数，后半部分为惩罚项。其中
glmnet有一些令人兴奋的特性，比如支持稀疏矩阵作为输入，可以限制待估系数的范围等，这些在我们的后续章节会进行详细的介绍。
它内部采用循环坐标下降算法来优化目标函数。而且它的底层是fortran写的，所以运行效率比不必担心，它执行起来非常快。
1.2 glmnet包安装
glmnet安装起来非常简单，因为已经发布到CRAN上了，我们可以直接用以下命令来安装：
install.packages("glmnet", repos = "https://cran.us.r-project.org")安装好以后便可以开启我们的glmnet包之旅了。
2 快速入门
这节介绍glmnet包中的主要函数以及它们的一般用法，对常用函数的输入参数以及输出结果做了简要的说明。
2.1 加载glmnet包
首先，我们载入glmnet包:
library(glmnet)这节中我们就以线性回归为例，来说明glmnet包的用法。
2.2 准备数据
接下来，我们载入包中自带的数据：
data(QuickStartExample)这时你会发现R的工作环境中新增了两个对象，x即输入矩阵，y即响应变量。
2.3 拟合模型
数据有了，我们就可以调用包中与之同名的glmnet函数来做线性回归了：
fit <- glmnet(x, y)该函数默认调用的就是线性回归，这里生成的结果“fit”是一个“glmnet”类的对象，它包含了模型拟合结果。本质上它是一个list对象，可以按照提取list元素的方式来提取，但不推荐。因为在R中，有一系列泛型函数可以来提取和使用模型对象中的数据， 比如plot, print, coef和predict等等。接下来，我们就用这些泛型函数来分析我们得到的模型。
2.4 模型对象的提取和可视化
如果想看到glmnet计算的每一步路径，可以采用print函数：
print(fit)##
## Call: glmnet(x = x, y = y)
##
## Df %Dev Lambda
## [1,] 0 0.00000 1.631000
## [2,] 2 0.05528 1.486000
## [3,] 2 0.14590 1.354000
## [4,] 2 0.22110 1.234000
## [5,] 2 0.28360 1.124000
## [6,] 2 0.33540 1.024000
## [7,] 4 0.39040 0.933200
## [8,] 5 0.45600 0.850300
## [9,] 5 0.51540 0.774700
## [10,] 6 0.57350 0.705900
## [11,] 6 0.62550 0.643200
## [12,] 6 0.66870 0.586100
## [13,] 6 0.70460 0.534000
## [14,] 6 0.73440 0.486600
## [15,] 7 0.76210 0.443300
## [16,] 7 0.78570 0.404000
## [17,] 7 0.80530 0.368100
## [18,] 7 0.82150 0.335400
## [19,] 7 0.83500 0.305600
## [20,] 7 0.84620 0.278400
## [21,] 7 0.85550 0.253700
## [22,] 7 0.86330 0.231200
## [23,] 8 0.87060 0.210600
## [24,] 8 0.87690 0.191900
## [25,] 8 0.88210 0.174900
## [26,] 8 0.88650 0.159300
## [27,] 8 0.89010 0.145200
## [28,] 8 0.89310 0.132300
## [29,] 8 0.89560 0.120500
## [30,] 8 0.89760 0.109800
## [31,] 9 0.89940 0.100100
## [32,] 9 0.90100 0.091170
## [33,] 9 0.90230 0.083070
## [34,] 9 0.90340 0.075690
## [35,] 10 0.90430 0.068970
## [36,] 11 0.90530 0.062840
## [37,] 11 0.90620 0.057260
## [38,] 12 0.90700 0.052170
## [39,] 15 0.90780 0.047540
## [40,] 16 0.90860 0.043310
## [41,] 16 0.90930 0.039470
## [42,] 16 0.90980 0.035960
## [43,] 17 0.91030 0.032770
## [44,] 17 0.91070 0.029850
## [45,] 18 0.91110 0.027200
## [46,] 18 0.91140 0.024790
## [47,] 19 0.91170 0.022580
## [48,] 19 0.91200 0.020580
## [49,] 19 0.91220 0.018750
## [50,] 19 0.91240 0.017080
## [51,] 19 0.91250 0.015570
## [52,] 19 0.91260 0.014180
## [53,] 19 0.91270 0.012920
## [54,] 19 0.91280 0.011780
## [55,] 19 0.91290 0.010730
## [56,] 19 0.91290 0.009776
## [57,] 19 0.91300 0.008908
## [58,] 19 0.91300 0.008116
## [59,] 19 0.91310 0.007395
## [60,] 19 0.91310 0.006738
## [61,] 19 0.91310 0.006140
## [62,] 20 0.91310 0.005594
## [63,] 20 0.91310 0.005097
## [64,] 20 0.91310 0.004644
## [65,] 20 0.91320 0.004232
## [66,] 20 0.91320 0.003856
## [67,] 20 0.91320 0.003513可以看到结有三列，第一列Df是非零系数的个数，第三列和第二列分别是λ以及该λ对应的解释偏差百分比%dev。
通过coef来提取模型的系数：
coef(fit,s=0.1)## 21 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.150928072
## V1 1.320597195
## V2 .
## V3 0.675110234
## V4 .
## V5 -0.817411518
## V6 0.521436671
## V7 0.004829335
## V8 0.319415917
## V9 .
## V10 .
## V11 0.142498519
## V12 .
## V13 .
## V14 -1.059978702
## V15 .
## V16 .
## V17 .
## V18 .
## V19 .
## V20 -1.021873704采用plot函数对拟合出的模型系数进行可视化：
plot(fit)上图中，每一条曲线代表一个变量的系数。该图展示了系数随L1-范数变化的图，图形顶部的横轴表示当前选定的λ下非零系数的个数。
2.5 预测
预测采用predict函数，参数newx用来设置输入数据，s用来设置λ的值：
nx <- matrix(rnorm(10*20),10,20)
predict(fit,newx=nx,s=c(0.1,0.05))## 1 2
## [1,] -1.17525458 -1.26294557
## [2,] -0.28860312 -0.46563599
## [3,] 1.15571622 1.18391784
## [4,] 3.48520956 3.63599739
## [5,] -0.94155250 -0.94327820
## [6,] -0.61023849 -0.50550565
## [7,] 0.07087784 0.08926638
## [8,] 0.13483451 -0.05207817
## [9,] -0.95596685 -0.89136569
## [10,] 3.79826343 3.980801122.6 交叉验证
glmnet提供了一系列的模型可供选择(随着λ的不同，估计出的模型也不尽相同)，而在大多数情况下我们需要从中挑选出一个最合适的来用就可以了。这时就需要通过交叉验证的方法来筛选最优的λ值了，cv.glmnet函数实现了这一功能。
继续沿用之前的样本数据，调用cv.glmnet函数：
cvfit <- cv.glmnet(x, y)可以看到，cv.glmnet返回的结果是一个cv.glmnet类的对象，该对象的类型和glmnet函数返回的结果一样，它们本质上都是R中的list。
我们用可视化的图形来展示cv.glmnet的结果：
plot(cvfit)该图中，红色的散点为交叉验证的散点图，横轴为logλ，纵轴为均方误差，每个点的标准偏差上界和下界也画出来了。图的顶部字数表示非零系数的个数，两条垂直的虚线即为交叉验证后选定的λ。
最优的λ 值可以直接采用如下命令来提取：
cvfit$lambda.min## [1] 0.07569327cvfit$lambda.1se## [1] 0.1322761其中，lambda.min是指交叉验证中使得均方误差最小的那个值λ，lambda.1se为离最小均方误差一倍标准差的λ值。
用coef函数来提取回归模型的系数：
coef(cvfit, s = "lambda.min")## 21 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.14867414
## V1 1.33377821
## V2 .
## V3 0.69787701
## V4 .
## V5 -0.83726751
## V6 0.54334327
## V7 0.02668633
## V8 0.33741131
## V9 .
## V10 .
## V11 0.17105029
## V12 .
## V13 .
## V14 -1.07552680
## V15 .
## V16 .
## V17 .
## V18 .
## V19 .
## V20 -1.05278699可以看到回归模型的系数是采用稀疏矩阵的形式来存储的。由于计算出的模型系数经常是稀疏的，这时采用稀疏矩阵的方式来存储和计算更有效率。如果你不习惯稀疏矩阵的输出形式，可以用as.matrix()将其转化为传统的矩阵形式。
预测同glmnet，直接采用predict泛型函数：
predict(cvfit, newx = x[1:5,], s = "lambda.min")## 1
## [1,] -1.3638848
## [2,] 2.5713428
## [3,] 0.5729785
## [4,] 1.9881422
## [5,] 1.5179882自此，glmnet的入门介绍完了，你可以用来他做一些基本的回归模型了。
接下来，我们对glmnet包进行更为深入的介绍。
3 线性回归
glmnet中的线性回归主要包含两类。一定是高斯簇gaussian，还有一类是多响应高斯簇mgaussian。我们依次介绍：
3.1 高斯簇
gaussian是glmnet函数中的默认函数簇，它本质上是带正则项的多元线性回归的估计问题。
3.1.1 优化目标
假定观测值
其中λ ≥ 0是模型复杂度参数；0 ≤ α ≤ 1,当α = 0时为岭回归，当α = 1时为lasso回归，0 ＜ α ＜ 1则为两者的折中。
3.1.2 glmnet参数设置
glmnet提供了很多参数可以供我们选择。下面介绍一些常用的参数设置：
alpha之前介绍过，它是弹性网的参数，取值范围是[0, 1]，当α = 1时是lasso回归，当α = 0时是岭回归。默认α = 1。weights配置观测的权重。默认每个观测的权重取值均为1。nlambdaλ取值的个数。默认值是100。lambda一般是程序自动构建，也可以自己定义。standardize表示在拟合模型前，x变量是否需要标准化。默认standardize=TRUE。
更多参数设置参考帮助文档help(glmnet)。
我们用下面的例子来看看这些参数的用法： 还是用原来的样本数据，不同的是取α = 0.2（接近岭回归的正则项），设置观测的权重以及减少λ序列的数量：
fit <- glmnet(x, y, alpha = 0.2, weights = c(rep(1,50),rep(2,50)), nlambda = 20)用print函数打印结果：
print(fit)##
## Call: glmnet(x = x, y = y, weights = c(rep(1, 50), rep(2, 50)), alpha = 0.2, nlambda = 20)
##
## Df %Dev Lambda
## [1,] 0 0.0000 7.939000
## [2,] 4 0.1789 4.889000
## [3,] 7 0.4445 3.011000
## [4,] 7 0.6567 1.854000
## [5,] 8 0.7850 1.142000
## [6,] 9 0.8539 0.703300
## [7,] 10 0.8867 0.433100
## [8,] 11 0.9025 0.266700
## [9,] 14 0.9101 0.164300
## [10,] 17 0.9138 0.101200
## [11,] 17 0.9154 0.062300
## [12,] 17 0.9160 0.038370
## [13,] 19 0.9163 0.023630
## [14,] 20 0.9164 0.014550
## [15,] 20 0.9164 0.008962
## [16,] 20 0.9165 0.005519
## [17,] 20 0.9165 0.003399打印结果之前已做过说明，这里不再赘述。可以看到这里λ并没有达到预设的20。这是因为在偏差解释率达到0.999或者其变化小于10e-5时计算就会终止。而这些预设的计算终止条件可以通过glmnet.control来设置，详见help(glmnet.control)。
快速入门中我们简单介绍了plot函数的用法，这里我们再详细介绍下plot函数中其他参数的设置：
3.1.3 plot参数设置
plot函数可以用xvar来定义X轴的度量，有三种选项：
“norm” 表示系数的L1-范数
“lambda” 表示对数lambda值
“dev” 表示偏差解释率
在plot函数中添加参数label = TRUE可以显示变量的标签。以下是选取不同参数设置下的plot：
plot(fit, xvar = "lambda", label = TRUE)plot(fit, xvar = "dev", label = TRUE)3.1.4 coef参数设置
coef函数中最常用的两个参数为:
s指定λ值exact表示是否需要根据提供的λ值进行精确地提取。如果exact = TRUE，那么就回去返回的结果中取精确地匹配λ，如果没有找到匹配项，那么模型会根据新的λ重新拟合模型；如果exact=FALSE(默认)，在没有找到精确匹配的λ 时，它会直接用线性插值来获取，不会重新去拟合模型。
以下是一个示例：
any(fit$lambda == 0.5)## [1] FALSEcoef.exact <- coef(fit, s = 0.5, exact = TRUE)
coef.apprx <- coef(fit, s = 0.5, exact = FALSE)
cbind2(coef.exact, coef.apprx)## 21 x 2 sparse Matrix of class "dgCMatrix"
## 1 1
## (Intercept) 0.19657091 0.199098747
## V1 1.17495906 1.174650452
## V2 . .
## V3 0.52934375 0.531934651
## V4 . .
## V5 -0.76126245 -0.760959480
## V6 0.46627405 0.468209413
## V7 0.06148079 0.061926756
## V8 0.38048793 0.380301491
## V9 . .
## V10 . .
## V11 0.14213614 0.143260991
## V12 . .
## V13 . .
## V14 -0.91090226 -0.911207368
## V15 . .
## V16 . .
## V17 . .
## V18 . 0.009196628
## V19 . .
## V20 -0.86099392 -0.863117051可以看出，当exact选取不同的参数时，提取的系数也存在一定程度的差异，但差距不大。没有特别要求的话，使用线性插值得到的结果已经够用了。
3.1.5 predict参数设置
predict函数与coef函数相比多了一些参数的设置：
newx是待预测的输入数据集。
type有多个选项可供选择：
“link” 给出预测值
“response” 对于gaussian簇，同“link”
“coefficients” 计算给定
s下的系数矩阵“nonzero” list对象，存储每个
s下非0系数对应的下标
predict(fit, newx = x[1:5,], type = "response", s = 0.05)## 1
## [1,] -0.9802591
## [2,] 2.2992453
## [3,] 0.6010886
## [4,] 2.3572668
## [5,] 1.7520421上述命令表示在λ = 0.05时计算x头5条观测的预测值。这里的s可以是一个向量，当s是一个多数值向量时，预测值则为一个矩阵。
3.1.6 交叉验证
这小节对cv.glmnet函数的参数做简要说明：
nfolds– 交叉验证数据集划分的份数foldid– 自定义划分数据type.measure– 定义交叉验证的损失函数，“deviance”和“mse”用的是平方损失，“mae”用的是平均绝对损失
# 做20重交叉验证，采用平均绝对损失
cvfit <- cv.glmnet(x, y, type.measure = "mse", nfolds = 20)除了可以设置nfolds来寻找合适的模型外，我们还可以通过foldid设置相同的数据划分来选择最优的α值。
foldid <- sample(1:10,size=length(y),replace=TRUE)
cv1 <- cv.glmnet(x,y,foldid=foldid,alpha=1)
cv.5 <- cv.glmnet(x,y,foldid=foldid,alpha=.5)
cv0 <- cv.glmnet(x,y,foldid=foldid,alpha=0)par(mfrow=c(2,2))
plot(cv1);plot(cv.5);plot(cv0)
plot(log(cv1$lambda),cv1$cvm,pch=19,col="red",xlab="log(Lambda)",ylab=cv1$name)
points(log(cv.5$lambda),cv.5$cvm,pch=19,col="grey")
points(log(cv0$lambda),cv0$cvm,pch=19,col="blue")
legend("topleft",legend=c("alpha= 1","alpha= .5","alpha0"),pch=19,col=c("red","grey","blue"))我们可以看到这里选择lasso(alpha=1)时，模型的均方误差最小。
cv.glmnet函数内部支持并行计算，但目前还不支持windows系统。如果你感兴趣，可以在除windows平台以外的系统上尝试，代码如下：
require(doMC)
registerDoMC(cores=2)
X = matrix(rnorm(1e4 * 200), 1e4, 200)
Y = rnorm(1e4)system.time(cv.glmnet(X, Y))## user system elapsed
## 2.440 0.080 2.518system.time(cv.glmnet(X, Y, parallel = TRUE))## user system elapsed
## 2.450 0.157 1.567可以看出并行计算对计算效率的提升比较明显。
另外，最新版本的glmnet中还新增了两个有意思的参数：
设定系数的上下界
比如，我们想要在拟合模型的同时，希望估计出的系数在-0.7到0.5之间，那么可以用upper.limits和lower.limits这两个参数：
tfit <- glmnet(x,y,lower.limits=-.7,upper.limits=.5)
plot(tfit)这里就会有很好的应用场景了，比如就就想看看模型中对其响应变量有正影响的那些变量进行建模，那么设置lower.limit=0就可以了。
设置参数值的时候也不是完全没有限制的，upper.limits的值不能小于0，lower.limits的值不能大于0。另外，如果想对每一个变量的系数对不同的限定，需要将这里的单点值即标量改为向量的形式就可以了。
设定惩罚因子
这个参数可以给每一个系数提供一个单独的惩罚因子。该惩罚因子默认是1，它也支持自定义。如果将惩罚因子全部设置成为0的话，相当于就没有惩罚项了。
看看以下公式就一目了然了：
这个参数设置选项很有用，假如我们知道了一些先验信息，知道了其中一些变量很重要，需要在建模正则化的同时一直保留这些变量，那么可以把这些变量对应系数的惩罚因子设置为0。
同样用之前的数据，我们把第5、10、15个变量对应的惩罚因子设置为0：
p.fac <- rep(1, 20)
p.fac[c(5, 10, 15)] <- 0
pfit <- glmnet(x, y, penalty.factor = p.fac)
plot(pfit, label = TRUE)从上图中可以看到，变量5、10、15对应的系数一直都在模型中。
还有一些其它的有用的参数，比如，exclude参数可以用来限制指定的变量入选模型；intercept参数可以用来设定模型是否含有截距项等等。更多设置参考帮助文档help(cv.glmnet).
3.2 多响应高斯簇
多响应高斯簇模型的估计需要在glmnet函数中设置family = "mgaussian"。与以上单变量响应模型类似，它只是响应变量增多了，我们通常称之为“多任务学习”问题。虽然响应变量增多了，但是建模时所选择的自变量是完全一样的，只是待估系数不同而已。
很显然，模型的因变量不再是一个向量形式，而是一个二维矩阵，这时估计出的系数也会是一个矩阵，先看看多响应高斯簇模型解决的问题：
这里的βj是p×K维的系数矩阵β的第j行。
3.2.1 载入演示数据
data(MultiGaussianExample)3.2.2 拟合模型
mfit <- glmnet(x, y, family = "mgaussian")多响应高斯模型与单响应高斯模型的大部分参数设置相同，如alpha, weights, nlambda, standardize。但是mgaussian簇有一个额外的参数standardize.response，它可以用来给响应变量做标准化，默认为FALSE。
3.2.3 查看拟合效果
用plot函数查看系数的变化：
plot(mfit, xvar = "lambda", label = TRUE, type.coef = "2norm")其中，xvar和label两个参数的设置同单响应的高斯模型一样。这里的type.coef = "2norm"表示每个变量的系数以二范数的形式展现。默认设置为type.coef = "coef"，这时每个响应变量会展示一张系数变化的图。
3.2.4 预测
同上，predict函数用来预测：
predict(mfit, newx = x[1:5,], s = c(0.1, 0.01))## , , 1
##
## y1 y2 y3 y4
## [1,] -4.7106263 -1.1634574 0.6027634 3.740989
## [2,] 4.1301735 -3.0507968 -1.2122630 4.970141
## [3,] 3.1595229 -0.5759621 0.2607981 2.053976
## [4,] 0.6459242 2.1205605 -0.2252050 3.146286
## [5,] -1.1791890 0.1056262 -7.3352965 3.248370
##
## , , 2
##
## y1 y2 y3 y4
## [1,] -4.6415158 -1.2290282 0.6118289 3.779521
## [2,] 4.4712843 -3.2529658 -1.2572583 5.266039
## [3,] 3.4735228 -0.6929231 0.4684037 2.055574
## [4,] 0.7353311 2.2965083 -0.2190297 2.989371
## [5,] -1.2759930 0.2892536 -7.8259206 3.2052113.2.5 交叉验证
同样的交叉验证用cv.glmnet函数：
cvmfit <- cv.glmnet(x, y, family = "mgaussian")画出交叉验证的结果：
plot(cvmfit)想要查看最优的λ,采用如下命令:
cvmfit$lambda.min## [1] 0.05193158cvmfit$lambda.1se## [1] 0.1316655《计量经济圈的圈子社群》
写在后面：各位圈友，一个等待数日的好消息，是计量经济圈应圈友提议，09月04日创建了“计量经济圈的圈子”知识分享社群，如果你对计量感兴趣，并且考虑加入咱们这个计量圈子来受益彼此，那看看这篇介绍文章和操作步骤哦（戳这里）。进去之后一定要看“群公告”，不然接收不了群信息。若需要获得计量经济学视频资料，那可以（戳这里）。