广义线性回归模型估计：所有线性回归的大仓库（2）
The following article is from 数据挖掘与R语言 Author 2万计量学者
4 逻辑回归
逻辑回归是分类问题中最常用的模型之一。如果是一个二分类问题，一般假定响应变量服从二项分布，如果是多分类问题，则假定服从多项式分布。
4.1 二项分布逻辑回归
假定响应变量的取值为
还可以改写为如下形式：
这个带惩罚逻辑回归的目标函数的对数似然如下：
4.1.1 载入示例数据集
data(BinomialExample)这里的输入x与其他分布簇相同，二分类逻辑回归的响应变量y是包含两个水平的因子对象。
4.1.2 拟合模型
做二分类逻辑回归还是用到glmnet函数，仅需要把函数簇改为family = "binomial"即可：
fit <- glmnet(x, y, family = "binomial")4.1.3 查看拟合效果
同样用plot查看结果：
plot(fit, xvar = "dev", label = TRUE)逻辑回归的预测同高斯簇函数的用法有点不同，主要体现在参数type的设置上，详细概括如下：
“link” 线性拟合值
“response” 拟合的概率值
“class” 给出计算出的最大概率对应的类的标签
“coefficients” 计算给定
s下的系数的估计值“nonzero” 返回一个list对象，该list包含每一个
s对应非零系数的索引
predict(fit, newx = x[1:5,], type = "class", s = c(0.05, 0.01))## 1 2
## [1,] "0" "0"
## [2,] "1" "1"
## [3,] "1" "1"
## [4,] "0" "0"
## [5,] "1" "1"4.1.5 交叉验证
逻辑回归的cv.glmnet的用法同高斯簇函数，nfolds, weights, lambda,parallel的设置一样，区别主要在type.measure：
“mse” 用平方损失
“deviance” 用真实偏差
“mae” 用平均绝对误差
“class” 用误分类率
“auc” ROC曲线的下面积(这个选项仅针对两分类逻辑回归)
例如，用误分类率误差为标准做十折交叉验证，代码如下：
cvfit <- cv.glmnet(x, y, family = "binomial", type.measure = "class")用plot查看cv.glmnet生成的结果： .
plot(cvfit)查看最优的λ值：
cvfit$lambda.min## [1] 0.0147554cvfit$lambda.1se## [1] 0.02349477coef和predict与高斯簇类似：
coef(cvfit, s = "lambda.min")## 31 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.24371065
## V1 0.06897051
## V2 0.66252343
## V3 -0.54275157
## V4 -1.13692797
## V5 -0.19143183
## V6 -0.95852340
## V7 .
## V8 -0.56528880
## V9 0.77454157
## V10 -1.45079395
## V11 -0.04363484
## V12 -0.06893736
## V13 .
## V14 .
## V15 .
## V16 0.36685293
## V17 .
## V18 -0.04013948
## V19 .
## V20 .
## V21 .
## V22 0.20882114
## V23 0.34014021
## V24 .
## V25 0.66310204
## V26 -0.33695819
## V27 -0.10570477
## V28 0.24317817
## V29 -0.22445233
## V30 0.11091002predict(cvfit, newx = x[1:10,], s = "lambda.min", type = "class")## 1
## [1,] "0"
## [2,] "1"
## [3,] "1"
## [4,] "0"
## [5,] "1"
## [6,] "0"
## [7,] "0"
## [8,] "0"
## [9,] "1"
## [10,] "1"值得一提的是，预测结果仅仅是针对响应变量的第二个水平。
4.2 多分类逻辑回归
多分类逻辑回归假定响应变量服从多项式分布, 假定响应变量有K个水平：
那么带弹性网惩罚项的非负对数似然函数如下：
这里β是一个p × K维的系数矩阵；βk是其中的第k列，表示第k个模型对应模型的系数；βj是第j行，表示第j个变量前面的系数。
后面的惩罚项||βj ||q有两种情形：当q=1时，它是一个lasso惩罚项；当q=2时，它是一个grouped-lasso惩罚项。
4.2.1 载入示例数据集
首先载入数据集：
data(MultinomialExample)4.2.2 拟合模型
多分类逻辑回归的模型拟和二分类逻辑回归类似，只是这里新增加了一个特殊的参数type.multinomial，当type.multinomial = "grouped"时，模型拟合时会让每个变量前面的系数全为0或者全不为零，其实就是对于每个类建立逻辑回归时所用到的变量完全相同。
fit <- glmnet(x, y, family = "multinomial", type.multinomial = "grouped")4.2.3 查看拟合效果
用plot查看模型拟合结果：
plot(fit, xvar = "lambda", label = TRUE, type.coef = "2norm")这里xvar和label的使用和之前相同，但是多了一个type.coef选项，这个选项仅适用于对多分类逻辑回归以及多响应变量线性回归，若type.coef = "coef"会用多张图分别展示每个响应变量的系数，若type.coef = "2norm"则会展示系数的L2-范数。
4.2.4 交叉验证和预测
我们同样也可以做交叉验证：
cvfit <- cv.glmnet(x, y, family="multinomial", type.multinomial = "grouped", parallel = TRUE)
plot(cvfit)同样地，我们也可以设置parallel = TRUE用并行计算来加速运算。
用拟合的模型来预测：
predict(cvfit, newx = x[1:10,], s = "lambda.min", type = "class")## 1
## [1,] "3"
## [2,] "2"
## [3,] "2"
## [4,] "3"
## [5,] "1"
## [6,] "3"
## [7,] "3"
## [8,] "1"
## [9,] "1"
## [10,] "2"5 泊松回归
泊松回归经常会用到计数模型中，假定其误差满足泊松分布。
经常用其均值的对数来建模：
给定
于是问题变成优化如下带惩罚的对数似然：
5.1 加载数据集
data(PoissonExample)5.2 拟合模型
采用glmnet函数，设置family = "poisson"：
fit = glmnet(x, y, family = "poisson")5.3 查看拟合效果
plot(fit)5.4 预测
用predict做预测,在参数选项的设置中，主要是type存在一些差异，做出说明如下：
“link” 给出线性拟合值
“response” 给出拟合的均值
“coefficients” 计算给定
s下的系数，也可以直接用coef函数“nonzero” 返回一个list对象，该list包含每一个
s对应非零系数的索引
coef(fit, s = 1)## 21 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.61123371
## V1 0.45819758
## V2 -0.77060709
## V3 1.34015128
## V4 0.04350500
## V5 -0.20325967
## V6 .
## V7 .
## V8 .
## V9 .
## V10 .
## V11 .
## V12 0.01816309
## V13 .
## V14 .
## V15 .
## V16 .
## V17 .
## V18 .
## V19 .
## V20 .predict(fit, newx = x[1:5,], type = "response", s = c(0.1,1))## 1 2
## [1,] 2.4944232 4.4263365
## [2,] 10.3513120 11.0586174
## [3,] 0.1179704 0.1781626
## [4,] 0.9713412 1.6828778
## [5,] 1.1133472 1.99345375.5 交叉验证
cvfit <- cv.glmnet(x, y, family = "poisson",type.measure = "mse")cv.glmnet中type.measure的设置：
“deviance” 偏差
“mse” 均方误差
“mae” 平均绝对误差
plot(cvfit)提取最优λ对应的模型系数：
opt.lam <- c(cvfit$lambda.min, cvfit$lambda.1se)
coef(cvfit, s = opt.lam)## 21 x 2 sparse Matrix of class "dgCMatrix"
## 1 2
## (Intercept) 0.017866985 0.129489929
## V1 0.622183525 0.590193059
## V2 -0.988415977 -0.952086326
## V3 1.530986777 1.487203667
## V4 0.233499975 0.214291103
## V5 -0.338975959 -0.315728229
## V6 0.001677868 .
## V7 -0.015428046 .
## V8 -0.001497926 .
## V9 -0.001313504 .
## V10 0.019472194 .
## V11 . .
## V12 0.031754913 0.026671248
## V13 -0.032423900 .
## V14 0.036634724 0.003544390
## V15 -0.008377899 .
## V16 0.020634049 .
## V17 . .
## V18 0.006320470 .
## V19 -0.036964656 -0.004107952
## V20 0.012224147 0.0110702446 Cox模型
这里采用生存分析的分析框架，假定数据的形式为
Cox模型是生存分析中常用的一种半参数回归模型，形式如下：
其中
跟其他模型一样，我们也可以给它加上一个弹性网惩罚项。
6.1 载入数据集
同样地，我们加载预先生成好的样本数据和响应变量，但是这里需要注意的是，这里采用了生存分析的分析框架，输入数据略有差别。首先，我们加载数据集：
data(CoxExample)head(y)## time status
## [1,] 1.76877757 1
## [2,] 0.54528404 1
## [3,] 0.04485918 0
## [4,] 0.85032298 0
## [5,] 0.61488426 1
## [6,] 0.29860939 0可以看到加载的数据还是包含自变量x以及响应变量两部分：x是维的矩阵；不同的是y，y为维的矩阵，其中名为“time”的列是观察时间，名为“status”的列为该观察时间下对应的状态，0表示右删失，1表示死亡。
6.2 拟合模型
同样用glmnet建模:
fit <- glmnet(x, y, family = "cox")6.3 查看效果
plot(fit)提取给定λ下对应的系数：
coef(fit, s = 0.05)## 30 x 1 sparse Matrix of class "dgCMatrix"
## 1
## V1 0.37693638
## V2 -0.09547797
## V3 -0.13595972
## V4 0.09814146
## V5 -0.11437545
## V6 -0.38898545
## V7 0.24291400
## V8 0.03647596
## V9 0.34739813
## V10 0.03865115
## V11 .
## V12 .
## V13 .
## V14 .
## V15 .
## V16 .
## V17 .
## V18 .
## V19 .
## V20 .
## V21 .
## V22 .
## V23 .
## V24 .
## V25 .
## V26 .
## V27 .
## V28 .
## V29 .
## V30 .Cox模型通常不是用来做预测，所以没有给出预测的样例，如果确实需要做预测，参考帮助文档help(predict.glmnet).
6.4 交叉验证
用cv.glmnet做K折交叉验证时，type.measure的选项仅支持“deviance”：
cvfit <- cv.glmnet(x, y, family = "cox")
plot(cvfit)提取最优的λ 值：
cvfit$lambda.min## [1] 0.01594373cvfit$lambda.1se## [1] 0.05343706提取模型系数：
coef(cvfit, s = "lambda.min")## 30 x 1 sparse Matrix of class "dgCMatrix"
## 1
## V1 0.491297030
## V2 -0.174600537
## V3 -0.218648771
## V4 0.175112406
## V5 -0.186673099
## V6 -0.490250398
## V7 0.335197365
## V8 0.091586519
## V9 0.450169329
## V10 0.115921793
## V11 .
## V12 .
## V13 0.017594799
## V14 .
## V15 .
## V16 .
## V17 -0.018364649
## V18 .
## V19 .
## V20 .
## V21 -0.002805776
## V22 -0.001422613
## V23 .
## V24 .
## V25 -0.023429230
## V26 .
## V27 0.001687662
## V28 .
## V29 .
## V30 -0.008236046稀疏矩阵
除了cox模型外，glmnet均支持稀疏矩阵作为输入，它的用法同常规矩阵的用法相同。
我们加载一个稀疏矩阵示例：
data(SparseExample)加载的数据x为100*20的一个稀疏矩阵，y为响应变量。
class(x)## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"创建稀疏矩阵有两种方式，一种方式是采用sparseMatrix生成；还有一种方式是直接采用Matrix来构建。
当输入是稀疏矩阵时，调用glmnet的方式跟普通矩阵没有差别：
fit <- glmnet(x, y)交叉验证也一样：
cvfit <- cv.glmnet(x, y)
plot(cvfit)稀疏矩阵除了可以用作glmnet的输入x，还可以用作predict函数的输入newx，我们来看看如下的例子：
i = sample(1:5, size = 25, replace = TRUE)
j = sample(1:20, size = 25, replace = TRUE)
x = rnorm(25)
nx = sparseMatrix(i = i, j = j, x = x, dims = c(5, 20))
predict(cvfit, newx = nx, s = "lambda.min")## 1
## [1,] 0.10001523
## [2,] 0.07808391
## [3,] -0.01152155
## [4,] -1.92250453
## [5,] 0.06792417自此，glmnet包的所有功能基本介绍完了，你可以用它来做各种回归模型了，考虑到lasso有变量选择的特性，简直不要太好用哦，赶紧试试吧^_^
参考资料
Jerome Friedman, Trevor Hastie and Rob Tibshirani. (2008).
Regularization Paths for Generalized Linear Models via Coordinate Descent
Journal of Statistical Software, Vol. 33(1), 1-22 Feb 2010.
Noah Simon, Jerome Friedman, Trevor Hastie and Rob Tibshirani. (2011).
Regularization Paths for Cox’s Proportional Hazards Model via Coordinate Descent
Journal of Statistical Software, Vol. 39(5) 1-13.
Robert Tibshirani, Jacob Bien, Jerome Friedman, Trevor Hastie, Noah Simon, Jonathan Taylor, Ryan J. Tibshirani. (2010).
Strong Rules for Discarding Predictors in Lasso-type Problems
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(2), 245-266.
Noah Simon, Jerome Friedman and Trevor Hastie (2013).
A Blockwise Descent Algorithm for Group-penalized Multiresponse and Multinomial Regression
(in arXiv, submitted)
《计量经济圈的圈子社群》
写在后面：各位圈友，一个等待数日的好消息，是计量经济圈应圈友提议，09月04日创建了“计量经济圈的圈子”知识分享社群，如果你对计量感兴趣，并且考虑加入咱们这个计量圈子来受益彼此，那看看这篇介绍文章和操作步骤哦（戳这里）。进去之后一定要看“群公告”，不然接收不了群信息。若需要获得计量经济学视频资料，那可以（戳这里）。