机器学习3-课后题:使用岭回归与lasso算法选择变量

1.题目

分别使用岭回归和Lasso解决薛毅书第279页(PDF为p331)例6.10的回归问题

2.准备

2.1.准备安装和加载包

使用到R语言的函数和对应包

说明：如何找哪个函数由哪个包提供：http://cran.rstudio.com/ ->TASK Views->Machine Learning-> search “关键字，比如lars”

执行代码如下

install.packages("lars") #http://cran.rstudio.com/ ->TASK Views->Machine Learning-> search lars
library(lars) #
library(ridge) #lm.ridge函数在ridge包
library(MASS) ##linearRidge函数在MASS包

2.2.读入数据

cement<-data.frame(
X1=c(7, 1,11,11,7,11,3,1,2,21,1,11,10),
X2=c(26, 29, 56, 31,52,55,71,31,54,47,40,66,68),
X3=c(6,15, 8, 8,6,9,17,22,18,4,23,9,8),
X4=c(60, 52, 20, 47,33,22,6,44,22,26,34,12,12),
Y =c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5,
93.1,115.9, 83.8, 113.3, 109.4) 
)

3.多重共线性检查

3.1.全部变量参与线性回归

summary(lm(Y~.,data=cement))

根据观察，截距与变量均不能通过系数检验，且残差标准差比较大

3.2.全部变量参与线性回归

求解样本方阵的条件数(kappa值) ，可以判断变量间具有严重的多重共线性。（k<100,说明共线性程度小,如果100< k< 1000,有较强的多重共线性,k>1000,在严重的多重共线）

>kappa(cor(cement),exact=TRUE)
[1] 2009.756

4.岭回归

4.1.全部变量做岭回归

k取值[0,0.1],步长0.001.求解不同k值下的岭回归估计族

library(MASS) #
lm.ridge(Y~.,cement,lambda=seq(0,0.1,0.001))
plot(lm.ridge(Y~.,cement,lambda=seq(0,0.1,0.001))) 

图形如下，有一条线系数绝对值始终很小，而且从正穿越到负，需要舍弃。
由于图形上没有标记这条线是哪个变量，结合下图各变量k[0,0.1]上的系数β值列表，只有X3是从正到负，可以判断这条线就是X3.

去掉X3后的线性回归：
summary(lm(Y~.-X3,data=cement))

根据观察，X2和X4的系数检验仍然不好，且残差标准差比较大，仍不满意，需要进一步做岭回归

4.1.去掉X3再做岭回归

根据岭回归的变量选择

• 中心化标准化后，可以剔除系数很小
• 随着增大，不稳定(正负穿梭)，震动趋于零的变量可以剔除

去掉X3再做岭回归：

结合图和β值列表，发现X4系数较小，可以删除X4变量

去掉X3和X4后后的线性回归：
summary(lm(Y~.-X3-X4,data=cement))

发现截距和变量的系数检验都是三颗星。
最终模型为：Y^ = 52.577 + 1.4683X1 + 0.6622X2.

5.使用linearRidge自动做岭回归

自动取了岭参数Ridge parameter: 0.01473162,舍弃X3

> summary(linearRidge(Y~.,data=cement));#自动取了岭参数Ridge parameter: 0.01473162,舍弃X3，

Call:
linearRidge(formula = Y ~ ., data = cement)


Coefficients:
            Estimate Scaled estimate Std. Error (scaled) t value (scaled) Pr(>|t|)    
(Intercept)  83.7040              NA                  NA               NA       NA    
X1            1.2922         26.3321              3.6721            7.171 7.45e-13 ***
X2            0.2977         16.0463              3.9883            4.023 5.74e-05 ***
X3           -0.1478         -3.2790              3.5979            0.911    0.362    
X4           -0.3506        -20.3290              3.9963            5.087 3.64e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ridge parameter: 0.01473162, chosen automatically, computed using 2 PCs

Degrees of freedom: model 3.01 , variance 2.837 , residual 3.183 

去掉X3后，发现X4还只有一个点

> summary(linearRidge(Y~.-X3,data=cement));#发现X4还只有一个点

Call:
linearRidge(formula = Y ~ . - X3, data = cement)


Coefficients:
            Estimate Scaled estimate Std. Error (scaled) t value (scaled) Pr(>|t|)    
(Intercept)  72.8790              NA                  NA               NA       NA    
X1            1.4436         29.4160              2.2433           13.113  < 2e-16 ***
X2            0.4003         21.5778              7.8141            2.761  0.00576 ** 
X4           -0.2501        -14.5015              7.8464            1.848  0.06458 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ridge parameter: 0.005856655, chosen automatically, computed using 2 PCs

Degrees of freedom: model 2.812 , variance 2.656 , residual 2.968 

去掉X3和X4后

> summary(linearRidge(Y~.-X3-X4,data=cement));#发现X4还只有一个点

Call:
linearRidge(formula = Y ~ . - X3 - X4, data = cement)


Coefficients:
            Estimate Scaled estimate Std. Error (scaled) t value (scaled) Pr(>|t|)    
(Intercept)  52.7494              NA                  NA               NA       NA    
X1            1.4629         29.8090              2.3450            12.71   <2e-16 ***
X2            0.6595         35.5511              2.3450            15.16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ridge parameter: 0.004852612, chosen automatically, computed using 2 PCs

Degrees of freedom: model 1.99 , variance 1.98 , residual 2 

结论：各变量的系数检验都是三颗星，最终的回归模型：
Y^ = 52.5749 + 1.4629X1 + 0.6595X2.

6.lasso

6.1.变量选择顺序

执行代码如下：

> w=as.matrix(cement)
> library(lars) #
> lars(w[,1:4],w[,5])

Call:
lars(x = w[, 1:4], y = w[, 5])
R-squared: 0.982 
Sequence of LASSO moves:
     X4 X1 X2 X3
Var   4  1  2  3
Step  1  2  3  4

6.2.选择哪些变量

plot(lars(w[,1:4],w[,5])),通过图表发现X3一直为零，结合summary，可以去掉x3

> summary(lars(w[,1:4],w[,5]))
LARS/LASSO
Call: lars(x = w[, 1:4], y = w[, 5])
  Df     Rss       Cp
0  1 2715.76 442.9167
1  2 2219.35 361.9455
2  3 1917.55 313.5020
3  4   47.97   3.0184
4  5   47.86   5.0000

在第三步cp指标(Mallows’s Cp)值最小,且残差RSS和比较小。
第三步的结果是x4, x1, x2变量，所以最终的变量选择是x4, x1, x2变量。

去掉X3后的再进行线性回归
summary(lm(Y~.-X3,data=cement))

发现X4的系数检验非常差，再去掉X4后(去掉X4不是通过lasso方法，而是lasso筛选过后，再通过线性回归的系数检验人工判断)，然后得到最终公式。