什么是正则化
线性模型的建模为了提高模型的泛化能力,一般会进行正则化处理,也就是在损失函数的构造上加上正则化项,如L1正则化项或者L2正则化项,L1正则化也就是常说的Lasso回归,将损失函数加上了L1范数,L2正则化就是Ridge回归,损失函数加上了L2范数。正则化项的大小是通过一个超参数(一般命名为lambda)控制,lambda越大则正则化项作用越强,拟合的模型系数会变小或变成0,这个超参数一般使用Cross-validation交叉验证来获取。
Lasso回归的特点是可以将模型中的一些参数系数缩小到0,起到筛选特征参数的作用,而Ridge回归则不会将任何模型项的系数降为0,但是Lasso回归有一个缺点,若变量中存在高度相关的变量组,则Lasso回归仅选择一个而忽视其他变量,就这一点而言,Ridge回归要优于Lasso回归。
为了同时保留Lasso的筛选模型参数的优点和Ridge回归会保留模型参数的优点,可以使用弹性网络(Elastic Net)回归进行兼顾,它使用一个超参数(一般都是命名为alpha)。alpha为0时,模型退化为Ridge回归,alpha为1时,模型退化为Lasso回归。同样的如果需要进行弹性网络拟合,则这个参数一般使用Cross-validation交叉验证来确定。
上式就是glmnet进行正则化拟合时使用的损失函数,关注一下式子中第二部分的正则化项,可以发现它是通过lambda来控制正则化项的大小,而具体的正则化项是一个alpha控制的L1和L2混合的正则化项,如果alpha等于1,则正则化项就是模型系数的L2范数,即为Ridge回归,如果alpha等于0,则正则化项就是模型系数的L1范数,即为Lasso回归。
什么是广义
最开始接触的线性回归的思想是从最小二乘法解决一个连续响应变量y和一个连续预测变量x发端,也就是一元线性回归,这种情况还是非常常见的,比如测定物质浓度时常用的标准曲线就是拟合一个浓度和吸光度的模型。
而这个思路可以很容易的推广到多元回归的,就是预测变量x是有多个特征,特征就是指的自变量,比如预测一个学生的数据成绩,可以使用的预测特征有学生做题时间、习题完成度、课堂注意时间等等。
再往后拓展就是如何拟合多次模型,比如平方项、立方项、交互作用项等等,其实有了多元回归的概念,平方项等高次项是很好解决的,先将相应的自变量运算得到相应的高次项,再将它也作为一个特征即可,比如需要拟合x1平方项,则可以先将x1的值进行平方,然后将其命名为一个新的特征如x1^2,令其参与到多元线性回归即可。
到目前为止响应变量和预测变量都是连续变量,如果预测变量是分类变量应该如何做,比如临床的风险因素:吸烟与否和饮酒与否都是分类变量?这个时候可以将分类变量编码为0 1等之类的数值变量,又叫做哑变量。
如果响应变量也不是连续变量,又要如何解决?这种情况下,同样会编码变量成哑变量,然后使用特定的连接函数来处理它,将其处理为连续变量。以logistics回归为例,它的连接函数是
,如果绘制这个函数的图像,则可以发现它的自变量在(0,1),函数值是(-Inf,Inf)。除了二分类,还有多分类、cox回归等各种情况,都可以通过使用连接函数变换后去使用线性回归。
使用glmnet进行正则化广义线性回归
代码语言:javascript复制library(glmnet)
library(tidyverse)
library(patchwork)
library(ggthemes)
data(BinomialExample)
x <- BinomialExample$x
y <- BinomialExample$y
导入必要的R包,使用glmnet自带的二分类测试数据集:BinomialExample进行logistics回归。
R代码很简单,使用glmnet函数,将family参数调整为binomial即可。
代码语言:javascript复制fit <- glmnet(x, y, family = "binomial")
plot(fit)
默认alpha值为1,也就是Loass回归,默认最大尝试100个lambda值,可以使用nlambda 参数控制最大尝试次数。
如果要挑选最佳lambda值,可以使用cv.glmnet函数进行交叉验证。交叉验证可以返回两种lambda值:lambda.min和lambda.1se,lambda.1se是指的在错误度量值最低的1个标准差内的最大lambda值。lambda越大,Lasso回归会倾向于简化更多的模型参数,可以获得更加精简的模型。
代码语言:javascript复制cvfit <- cv.glmnet(x, y, family = "binomial")
plot(cvfit)
cvfit$lambda.min
# [1] 0.02140756
cvfit$lambda.1se
# [1] 0.05427596
# 使用lambda.min重新拟合获得模型
opti_fit <- glmnet(x, y, family = "binomial", lambda = cvfit$lambda.min)
# 模型系数,以下两种方式均可
opti_fit$beta
coef(opti_fit)
如果需要挑选合适的alpha值,也就是进行弹性网络拟合,可以手动进行交叉验证(glmnet不支持对alpha自动交叉验证):
代码语言:javascript复制# 参数搜索
alpha_seq <- seq(0, 1, by = 0.1)
# 使用10-fold交叉验证,因此将样本分配为10个fold编号
foldid <- sample(1:10, nrow(x), replace = TRUE)
# 1. 每个alpha值进行一次交叉验证
# 返回结果:
# cvm:就是这10次交叉验证的错误度量平均值,常规线性模型默认使用Deviance,也就是MSE(平均标准误差),logistics回归是使用Bionomical Deviance
# cvsd:10次交叉验证的错误度量标准差
# lambda: 尝试的lambda值
# index_min:最低错误率模型对应的lambda值
# index_1se:错误率可接受的最精简模型对应的lambda值
cv_model_ls <-
alpha_seq %>%
set_names(., .) %>% # 对向量添加names为自身,保证map返回的列表也是有names的
map(function(alpha){
cv.model <- cv.glmnet(x, y, family = "binomial", foldid = foldid, alpha = alpha)
list(
cvm = cv.model$cvm,
cvsd = cv.model$cvsd,
lambda = cv.model$lambda,
index_min = cv.model$index["min",1],
index_1se = cv.model$index["1se",1]
)
}
)
# 2. 手动绘制11个alpha值下的lambda分布曲线,见下图
cv_model_ls %>%
map2(names(.), function(model, alpha){
labmda_min <- model$lambda[model$index_min] %>% log
labmda_1se <- model$lambda[model$index_1se] %>% log
data.frame(x = log(model$lambda),
y = model$cvm,
std = model$cvsd) %>%
ggplot(aes(x = x, y = y))
geom_point(color = "Red")
geom_errorbar(aes(ymin = y - std, ymax = y std), color = "gray")
geom_vline(xintercept = c(labmda_min, labmda_1se) , linetype = 3)
labs(x = expression(paste("Log(", lambda, ")")), y = "Binomial Deviance")
theme_bw()
ggtitle(paste0("alpha: ", alpha ))
}) %>%
wrap_plots(
ncol = 4
)
将11个模型绘制到同一张图中展示,见下图,从错误度量来说,alpha=1时的模型最佳,具有最低的错误程度。由于alpha=1恰好就是上面的Lasso交叉验证回归模型opti_fit,所以就不需要再进行一次glmnet拟合了,一般情况下需要根据最佳alpha和lambda值重新进行一次glmnet获取模型。
代码语言:javascript复制# 3. 将11个alpha的交叉验证结果进行展示
p <- cv_model_ls %>%
map(function(x) list(cvm = x$cvm, lambda = x$lambda) ) %>%
map(as_tibble) %>%
enframe(name = "alpha") %>%
unnest(value) %>%
mutate(lambda = log(lambda)) %>%
ggplot(aes(x = lambda, y = cvm))
geom_line(aes(color = alpha), size = 1)
labs(x = expression(paste("Log(", lambda, ")")), y = "Binomial Deviance")
ggthemes::theme_base()
# 展示原图及局部放大图
p (p xlim(c(-5, 0)) ylim(c(0.7, 1)))
guide_area()
plot_layout(
guides = 'collect',
ncol = 3,
widths = c(5, 5, 1)
)
# 最佳模型为alpha=1时的lambda为0.02349477的模型
cv_model_ls$`1` %>% {.$lambda[.$index_min]}
# [1] 0.02349477
本例中的自变量x的各个特征的相关性并不强,见下图,因此也并非一定要使用弹性网络或者Ridge回归进行拟合,Lasso回归的模型是比较不错的,它可以简化模型,筛去一部分不重要的特征。
代码语言:javascript复制cor(x) %>%
corrplot::corrplot(
method = "color",
order = "hclust",
tl.pos = "n",
cl.pos = "n"
)