R语言笔记之广义线性模型压缩方法1

glmnet包可以对一系列调优参数值同时计算参数估计。
该包可以用于线性回归，也可以拟合广义线性模型，如逻辑回归，多项式回归，泊松回归，cox回归。

初始glmnet

> install.packages("glmnet")
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.4/glmnet_2.0-13.zip'
Content type 'application/zip' length 1769153 bytes (1.7 MB)
downloaded 1.7 MB

package ‘glmnet’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\LLJiang\AppData\Local\Temp\RtmpCswyie\downloaded_packages
> library(glmnet)
载入需要的程辑包：Matrix
载入需要的程辑包：foreach
Loaded glmnet 2.0-13

Warning messages:
1: 程辑包‘glmnet’是用R版本3.4.3 来建造的 
2: 程辑包‘foreach’是用R版本3.4.3 来建造的 
> dat=read.csv("https://raw.githubusercontent.com/happyrabbit/DataScientistR/master/Data/SegData.csv")
#对数据进行一些清理，删除错误的样本观测，消费金额不能为负数
> dat=subset(dat,store_exp>0&online_exp>0)
> trainx=dat[,grep("Q",names(dat))]
#将实体店消费量和在线消费之和当作因变量
#得到总消费量=实体店消费+在线消费
> trainy=dat$store_exp+dat$online_exp
> glmfit=glmnet(as.matrix(trainx),trainy)
>

我们可以通过plot(),coef(),predict()等函数来得到相应的信息。

> plot(glmfit,label=T)
>

这里写图片描述
图中每种颜色的线代表对应一个自变量，展示的是随着Lasso罚函数对应调优参数的变化，各个变量对应的参数估计路径。
图中有上下两个X轴标度，下X轴是调优参数变化对应最优解的一阶范数值，上X轴是调优参数对应的非0参数估计个数，也就是Lasso模型的自由度。

> print(glmfit)

Call:  glmnet(x = as.matrix(trainx), y = trainy) 

      Df   %Dev   Lambda
 [1,]  0 0.0000 3042.000
 [2,]  2 0.1038 2771.000
 [3,]  2 0.1919 2525.000
 [4,]  2 0.2650 2301.000
 [5,]  3 0.3264 2096.000
 [6,]  3 0.3894 1910.000
 [7,]  3 0.4417 1741.000
 [8,]  3 0.4852 1586.000
 [9,]  3 0.5212 1445.000
[10,]  3 0.5512 1317.000
[11,]  3 0.5760 1200.000

如上，第一例Df表示非0估计的参数个数
第二列%Dev解释方差百分比
第三列表示Lambda调优参数的取值。
当%Dev百分比只发生微小变化时，算法会提前停止。

> coef(glmfit,s=1200)
11 x 1 sparse Matrix of class "dgCMatrix"
                    1
(Intercept) 2255.2221
Q1          -390.9214
Q2           653.6437
Q3           624.4068
Q4             .     
Q5             .     
Q6             .     
Q7             .     
Q8             .     
Q9             .     
Q10            .     
>

如上，当最优参数为1200时，只有3个变量的参数估计非0.

> newdat=matrix(sample(1:9,30,replace=T),nrow=3)
> predict(glmfit,newdat,s=c(1741,2000))
            1        2
[1,] 6331.299 6142.796
[2,] 3662.438 3921.191
[3,] 6242.257 6445.598
>

如上，结果中每列分别对应一个最优参数取值的预测。
需要通过交互校验进行参数调优。

> cvfit=cv.glmnet(as.matrix(trainx),trainy)
> plot(cvfit)
>

这里写图片描述
如上图，红色的点是不同调优参数取值对应的交互校验均方误差，灰色的线是相应置信区间。
两条虚线表示选中的两个调优参数。左边的那个调优参数值对应的是最小交互校验均方误差，右边的那个调优参数值是离最小均方误差一个标准差的调优参数值。

> cvfit$lambda.min
[1] 13.79349
> cvfit$lambda.1se
[1] 1316.656
>

也可以采用如下方式查看不同调优参数值对应的回归参数估计

> coef(cvfit,s="lambda.1se")
11 x 1 sparse Matrix of class "dgCMatrix"
                    1
(Intercept) 2221.0392
Q1          -340.6045
Q2           629.6860
Q3           585.1639
Q4             .     
Q5             .     
Q6             .     
Q7             .     
Q8             .     
Q9             .     
Q10            .     
>

收缩线性回归

线性回归又两种，一种是高斯家族模型，其中因变量是一个向量。
另一种是多元高斯，也就是多元响应变量情况，此时因变量是一个矩阵，参数也是矩阵。

> fit=glmnet(as.matrix(trainx),trainy,alpha=0.2,nlambda=20)
> print(fit,digits=4)

Call:  glmnet(x = as.matrix(trainx), y = trainy, alpha = 0.2, nlambda = 20) 

      Df   %Dev    Lambda
 [1,]  0 0.0000 15210.000
 [2,]  4 0.2502  9366.000
 [3,]  6 0.4590  5768.000
 [4,]  7 0.5848  3552.000
 [5,]  9 0.6502  2188.000
 [6,]  9 0.6823  1347.000
 [7,]  9 0.6967   829.700
 [8,]  9 0.7033   511.000
 [9,]  9 0.7064   314.700
[10,]  9 0.7080   193.800
[11,]  9 0.7088   119.300
[12,]  9 0.7093    73.500
[13,]  9 0.7095    45.270
[14,]  9 0.7096    27.880
[15,]  9 0.7096    17.170
[16,]  9 0.7096    10.570
[17,] 10 0.7096     6.511
[18,] 10 0.7096     4.010
> plot(fit,xvar="lambda",label=T)
>

这里写图片描述
由上图，可以看到，随着调优参数的增大，参宿逐步向0收缩。
Q1，Q2，Q3，Q8参数估计差不多在最后才收缩为0.
说明这几个变量对解释因变量最重要。

> plot(fit,xvar="dev",label=T)
>

这里写图片描述
上图为解释方差的百分比的路径图
其中横坐标为解释方差百分比时的路径图和之前不太一样。
非0参数个数越多，即使方差百分比就越大。
到最右端的时候，解释方差百分比变化很小，但是参数估计却急剧增大。
由上图我们可以看出，Q7对模型拟合几乎不起作用。
通过以下代码可以查看你想要的取值是不是已经调优拟合过

> any(fit$lambda==1000)
[1] FALSE
>

上面结果表明调优参数=1000这个取值不在调优取值内。
通常情况下用线性插值就可以，不需要重新拟合模型。
我们可以在新的数据集中预测结果

> newdat=matrix(sample(1:9,30,replace=T),nrow=3)
> predict(fit,newdat,type="response",s=1000)
             1
[1,]  5674.221
[2,] 10990.417
[3,]  3684.514

> predict(fit,newdat,
#得到的是线性插值参数估计
type="coefficients",s=1000)
11 x 1 sparse Matrix of class "dgCMatrix"
                     1
(Intercept) 1177.37262
Q1          -483.29175
Q2           850.50376
Q3           645.89528
Q4            60.82823
Q5           163.98503
Q6           192.23258
Q7             .      
Q8          -193.76827
Q9           198.42621
Q10         -132.20608
>

type=”nonzero”将会返回一个向量告诉你，那些自变量的参数估计非0

> predict(fit,newdat,type="nonzero",s=1000)
  X1
1  1
2  2
3  3
4  4
5  5
6  6
7  8
8  9
9 10
>

type.measure=”class”仅仅适用于二项回归和多项逻辑回归，其使用的是误判率。
type.measure=”auc”只针对二分类逻辑回归，使用的是ROC线下的面积。
type.measure=”mae”绝对误差均值可以用于除了cox模型之外的其他所有模型。

> cvfit=cv.glmnet(as.matrix(trainx),trainy,type.measure="mae",nfolds=20)
> cvfit$lambda.min
[1] 7.191939
> predict(cvfit,newx=newdat,s="lambda.min")
             1
[1,]  8992.229
[2,] 11106.268
[3,]  7238.749
> coef(cvfit,s="lambda.min")
11 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept) -2410.85923
Q1           -301.81607
Q2            675.95354
Q3            689.31828
Q4            218.49555
Q5            600.75548
Q6            673.01104
Q7              .      
Q8            -90.21178
Q9            359.73608
Q10          -299.33332
>

此外，也可以自己指定交互校验的样本分层情况

> foldid=sample(1:10,size=length(trainy),replace=T)
> cv1=cv.glmnet(as.matrix(trainx),trainy,foldid=foldid,alpha=1)
> cv.2=cv.glmnet(as.matrix(trainx),trainy,foldid=foldid,alpha=.2)
> cv.5=cv.glmnet(as.matrix(trainx),trainy,foldid=foldid,alpha=.5)
> cv0=cv.glmnet(as.matrix(trainx),trainy,foldid=foldid,alpha=0)
> plot(log(cv1$lambda),cv1$cvm,pch=19,col=2,xlab="log(Lambda)",ylab=cv1$name)
>

这里写图片描述

> points(log(cv.5$lambda),cv.5$cvm,col=1)
> points(log(cv.2$lambda),cv.2$cvm,col=3)
> points(log(cv0$lambda),cv0$cvm,col=4)
> legend("topleft",legend=c("alpha=1","alpha=0.5","alpha=0.2","alpha=0"),pch=19,col=c(2,1,3,4))
>

这里写图片描述
由上图所示，不同alpha的取值对应的最优均方误差几乎相同。