1.利用GCV(广义交叉验证)实现最优岭回归参数选择
#根据GCV方法,获取最佳岭参数k
#x:自变量的数据矩阵
#y:响应变量向量或矩阵
#kmax:岭参数的最大值
#qnum:根据0~kmax分成qnum等分
#intT:是否计算矩阵
getBestK <-function(X,Y,kMax=1,qnum=10,intT=TRUE){
if(intT)
X <-cbind(t0=1,X)
kvals=c(0,(1:qnum)*(kMax/qnum))
n=nrow(X)
glist <-NULL
for (k in kvals) {
mk=X%*%solve(t(X)%*%X+k*diag(ncol(X)))%*%t(X)
yk=mk%*%Y
Xs<-svd(X)
d <-Xs$d
dx <-length(d)
div <-d^2 +rep(k,rep(dx,1))
GCV <-sum((Y-yk)^2)/(n-sum(matrix(d^2/div,dx)))^2
glist <-c(glist,GCV)}
return(list(k=kvals[which.min(glist)],gval=min(glist),glist=glist))
}
引用getBestK函数,实现最优岭回归参数k的选择
library(MASS)
x=as.matrix(data[,1:3])
y=data[,4]
m0=getBestK(x,t(t(y)),kMax = 1,qnum = 2000)
m0$k
结果如下:
2.实现岭回归分析的函数
在进行岭回归时,发现可以使用两个函数,lm.Ridge及linearRidge,都可以实现。
首先读入数据:
data <-read.csv("D:\\sample\\sample_2.csv",header=T)
(1)lm.Ridge()函数
实现的代码分别如下:
Library(MASS)
lr<-lm.ridge(N ~ SHDI + MIDU + LSI + CONTAG,data=data5,lambda = 0.0085,model = TRUE) #把lambda设置为0.0085,是通过下面的代码求得的最佳岭回归参数,但获取岭回归参数的方法较多,同时也会造成结果不同,如果lambda不进行设置,会默认为0
lr #显示模型的参数系数,不能显示R^2等误差的情况,也不能使用summary()函数
结果:
对于k值得选择除了,最上面的getBestK函数外,还可以通过岭迹曲线得到,但是在k值选择时,受人为主观意识影响大,获取lm.ridge函数得到的模型对应的岭迹曲线,代码如下:
library(MASS)
ridge.sol <- lm.ridge(N ~ SHDI + MIDU + LSI + CONTAG, lambda = seq(0,1, length = 2000), data = data5, model = TRUE) #landad为使用seq函数把0-1范围均等分割为2000分,得到不同的lambda
画出图形,得到四个自变量的系数在lambad变化的情况下的变化情况,当所有系数趋于平滑时,可以选取这时的lambad值:
matplot(x=ridge.sol$lambda, y=t(ridge.sol$coef), xlab = expression(lamdba), ylab= "Cofficients", type = "l", lty = 1:20) # lty = 1:20可加可不加,设置线的形状.
#作出lambdaGCV取最小值时的那条竖直线
abline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])
结果:
发现,上图在lambad在0.8左右,自变量的系数值趋于稳定。
下面的语句绘出lambda同GCV之间关系的图形:
plot(ridge.sol$lambda, ridge.sol$GCV, type = "l", xlab = expression(lambda), ylab = expression(beta))
abline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)]) #语句ridge.sol$coef[which.min(ridge.sol$GCV)] 为找到GCV最小时对应的系数
结果:
在上面的代码中,还可以调整lambda = seq(0, 1, length = 2000)的范围及大小,如果lambad只是一个值,则得到的图像为空,只有在lambad变化的时候,才能得到岭迹曲线。
(2)linearRidge()函数
linearRidge()函数也可以用于求岭回归,如果lambad属性默认,则该函数可以自动选取岭回归参数,同时也可以自己通过其他的方式选择好,再进行设置;lm.Ridge函数不能自动选取岭回归参数,如果不对lambad进行设置,则默认为0. linearRidge函数实现岭回归的代码如下:
library(ridge)
mod <- linearRidge(N ~ SHDI + MIDU+LSI+CONTAG,lambda = 0.3,data = data5)
summary(mod) #可以查看判断结果,lm.Ridge使用该函数只能得带自变量系数;
结果如下:
从模型运行结果看,测岭回归参数值为0.3(即k值),是自己设置的lambad,且选择恰当的k值,进行岭回归会提高结果的显著性,因此在进行岭回归时,k值得选择还是挺重要的。
针对linearRidge函数,得到岭迹曲线,因在lambad取默认值的时候,lambad是会变化的,因此,岭迹曲线比较容易得到,代码如下:
plot(mod <- linearRidge(N ~ SHDI + MIDU+LSI+CONTAG, data = data5)) #没有设置lambad,默认
结果:
同时,你也可以在得到岭迹曲线时,自己设置lambad的变化范围,如下:
plot(mod <- linearRidge(N ~ SHDI + MIDU+LSI+CONTAG,lambda = seq(0,1,length=2000), data = data5))
结果:
如果lambad等于一个固定的值,则图像中只有几个点,点数对应于自变量的数目:
plot(mod <- linearRidge(N ~ SHDI + MIDU+LSI+CONTAG,lambda = 0.028, data = data5))
结果:
总结:
lm.Ridge和linearRidge是有区别的,再具体的差异就不是很清楚了,只能从上面几点上进行总结。且两者在lambad设置相同的情况下,得到的自变量系数结果也是不一样的。
注:
data为自己的数据集,getBestK函数代码,引用自《线性回归及其优化》。