R语言实现岭回归

##########岭回归##########
岭回归的关键是寻找合适的岭参数k，下面介绍了三种选择k值的方法。并且介绍了如何进行岭回归。需要注意的是，这里的参数lamdba就代表k值，当你进行回归时，需要对这个参数进行变动。

##########选择合适的k值
#方法一
#getbestk函数参数解释如下：
#根据GCV方法，获取最佳岭参数k
#x:自变量的数据矩阵
#y：响应变量向量或矩阵
#kmax:岭参数的最大值
#qnum：根据0~kmax分成qnum等分
#intT:是否计算矩阵getBestK <-function(X,Y,kMax=1,qnum=10,intT=TRUE){
 

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))
}
x=as.matrix(data[,1:16])#根据你的数据选择自变量与因变量
y=data[,17]#根据你的数据选择自变量与因变量
m0=getBestK(x,t(t(y)),kMax = 1,qnum = 2000)
m0$k

#方法二
library(MASS)
library(ggplot2)

library(MASS)
library(ggplot2)

ridge.sol <- lm.ridge(Churn~., lambda = seq(0,1, length = 2000), data = data, model = TRUE) 
#landad为使用seq函数把0-1范围均等分割为2000分，得到不同的lambda
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)])
#下面的语句绘出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最小时对应的系数

#方法三

#这个方法可以自动寻找最优的k值

library(ridge)
mod <- linearRidge(Churn~.,data = data)
summary(mod) #可以查看判断结果，lm.Ridge使用该函数只能得带自变量系数；

##########进行回归##########

# 岭回归
X <- as.matrix(data[,1:16])#生成自变量与因变量
y <- data[, 17]
fit <- linearRidge(Churn~., data =data, lambda =0.984)

# 计算系数的标准误
fit.se <- sqrt(diag(solve(t(X) %*% X + fit$lambda * diag(ncol(X)))))
fit$se <- fit.se
# 将标准误为0的值替换为1e-10
fit.se[fit.se==0] <- 1e-10

# 计算系数的P值
t_value <- fit$coef / fit$se
df <- nrow(X) - ncol(X)
p_value <- 2 * pt(abs(t_value), df = df, lower.tail = FALSE)

# 将系数、标准误和P值合并到一个数据框中
result <- data.frame(lambda = fit$lambda, coef = fit$coef, se = fit$se, p_value = p_value)

result

# 绘制岭回归系数的棒图
ggplot(result, aes(x = 1:length(fit$coef), y = coef, fill = p_value < 0.05)) +
  geom_bar(stat = "identity", width = 0.5, color = "black") +
  geom_errorbar(aes(ymin = coef - se, ymax = coef + se), width = 0.2) +
  theme_bw() +
  labs(title = "Ridge Regression Coefficients and Standard Errors",
       x = "Variable Index", y = "Coefficient")

#利用岭回归进行预测
prediction<- predict(fit, newx = newx, type = "response", s=0.984)