R语言—岭回归实现函数

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函数代码，引用自《线性回归及其优化》。