最近在看特征的选择,看到lasso对特征选择不错,下面直接上干货
数据为广州统计年检2015年数据
目标:
1) 梳理影响地方财政收入的关键特征,分析、识别影响地方财政收入的关键特征的选择模型;
2) 结合目标1的因素分析,对广州市2015年的财政总收入及各个类别收入进行预测。
下面为R语言代码部分
head(data)
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 3831732 181.54 448.19 7571.00 6212.70 6370241 525.71 985.31 60.62
## 2 3913824 214.63 549.97 9038.16 7601.73 6467115 618.25 1259.20 73.46
## 3 3928907 239.56 686.44 9905.31 8092.82 6560508 638.94 1468.06 81.16
## 4 4282130 261.58 802.59 10444.60 8767.98 6664862 656.58 1678.12 85.72
## 5 4453911 283.14 904.57 11255.70 9422.33 6741400 758.83 1893.52 88.88
## 6 4548852 308.58 1000.69 12018.52 9751.44 6850024 878.26 2139.18 92.85
## x10 x11 x12 x13 y
## 1 65.66 120.0 1.029 5321 64.87
## 2 95.46 113.5 1.051 6529 99.75
## 3 81.16 108.2 1.064 7008 88.11
## 4 91.70 102.2 1.092 7694 106.07
## 5 114.61 97.7 1.200 8027 137.32
## 6 152.78 98.5 1.198 8549 188.14
Data <- data
###数据概括性度量
Min=sapply(Data,min) #最小值
Max=sapply(Data,max) #最大值
Mean=sapply(Data,mean) #均值
SD=sapply(Data,sd) #方差
cbind(Min,Max,Mean,SD)
## Min Max Mean SD
## x1 3831732.000 7599295.000 5579519.9500 1.262195e+06
## x2 181.540 2110.780 765.0350 5.956983e+02
## x3 448.190 6882.850 2370.8250 1.919167e+03
## x4 7571.000 42049.140 19644.6850 1.020302e+04
## x5 6212.700 33156.830 15870.9480 8.199771e+03
## x6 6370241.000 8323096.000 7350513.6000 6.213419e+05
## x7 525.710 4454.550 1712.2390 1.184714e+03
## x8 985.310 15420.140 5705.7990 4.478400e+03
## x9 60.620 228.460 129.4935 5.050983e+01
## x10 65.660 852.560 340.2165 2.515779e+02
## x11 97.500 120.000 103.3050 5.513283e+00
## x12 1.029 1.906 1.4222 2.532348e-01
## x13 5321.000 41972.000 17273.8000 1.110919e+04
## y 64.870 2088.140 618.0840 6.092545e+02
#pearson相关系数,保留两位小数
round(cor(Data,method = c("pearson")),2)
## x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11
## x1 1.00 0.95 0.95 0.97 0.97 0.99 0.95 0.97 0.98 0.98 -0.29
## x2 0.95 1.00 1.00 0.99 0.99 0.92 0.99 0.99 0.98 0.98 -0.13
## x3 0.95 1.00 1.00 0.99 0.99 0.92 1.00 0.99 0.98 0.99 -0.15
## x4 0.97 0.99 0.99 1.00 1.00 0.95 0.99 1.00 0.99 1.00 -0.19
## x5 0.97 0.99 0.99 1.00 1.00 0.95 0.99 1.00 0.99 1.00 -0.18
## x6 0.99 0.92 0.92 0.95 0.95 1.00 0.93 0.95 0.97 0.96 -0.34
## x7 0.95 0.99 1.00 0.99 0.99 0.93 1.00 0.99 0.98 0.99 -0.15
## x8 0.97 0.99 0.99 1.00 1.00 0.95 0.99 1.00 0.99 1.00 -0.15
## x9 0.98 0.98 0.98 0.99 0.99 0.97 0.98 0.99 1.00 0.99 -0.23
## x10 0.98 0.98 0.99 1.00 1.00 0.96 0.99 1.00 0.99 1.00 -0.17
## x11 -0.29 -0.13 -0.15 -0.19 -0.18 -0.34 -0.15 -0.15 -0.23 -0.17 1.00
## x12 0.94 0.89 0.89 0.91 0.90 0.95 0.89 0.90 0.91 0.90 -0.43
## x13 0.96 1.00 1.00 1.00 0.99 0.94 1.00 1.00 0.99 0.99 -0.16
## y 0.94 0.98 0.99 0.99 0.99 0.91 0.99 0.99 0.98 0.99 -0.12
## x12 x13 y
## x1 0.94 0.96 0.94
## x2 0.89 1.00 0.98
## x3 0.89 1.00 0.99
## x4 0.91 1.00 0.99
## x5 0.90 0.99 0.99
## x6 0.95 0.94 0.91
## x7 0.89 1.00 0.99
## x8 0.90 1.00 0.99
## x9 0.91 0.99 0.98
## x10 0.90 0.99 0.99
## x11 -0.43 -0.16 -0.12
## x12 1.00 0.90 0.87
## x13 0.90 1.00 0.99
## y 0.87 0.99 1.00
#加载adapt-lasso源代码
lasso.adapt.bic2<-function(x,y){
require(lars)
ok<-complete.cases(x,y)
x<-x[ok,]
y<-y[ok]
m<-ncol(x)
n<-nrow(x)
x<-as.matrix(x)
one <- rep(1, n)
meanx <- drop(one %*% x)/n
xc <- scale(x, meanx, FALSE)
normx <- sqrt(drop(one %*% (xc^2)))
names(normx) <- NULL
xs <- scale(xc, FALSE, normx)
out.ls=lm(y~xs)
beta.ols=out.ls$coeff[2:(m+1)]
w=abs(beta.ols)
xs=scale(xs,center=FALSE,scale=1/w)
object=lars(xs,y,type="lasso",normalize=