多元回归中常见的变量选择方法及其在R中实现过程

多元回归中常用的变量选择方法及其在R中的实现

在多元回归中，有时预测变量太多，需要想办法减少预测变量的数量。一般来说，减少预测变量的数量可能有两个不冲突的原因：（1）寻求简约模型，有利于解释变量之间的关系；（2）预测变量过多）可能会导致模型混乱，例如某些预测变量可能存在强线性相关，即共线性问题，可能导致回归系数不稳定。

变量选择是从所有预测变量的集合中选择一个预测变量子集进行回归的过程，其目标是在不损失尽可能多的解释力的情况下减少候选预测变量的数量。

对于这个过程，可以根据回归中每个预测变量的回归系数的显着性来手动判断选择。还有一些自动机器选择方法，如逐步法、全子集回归等。

下面举例介绍这些变量选择方法在R语言中的实现过程。

注意：下面提到的变量选择方法可以用于多种类型的回归，例如广义线性模型。

示例数据，R代码的百度盘链接（提取码，8x3y）：

如果百度盘出现故障，也可以从GitHub备份中获取：

多元线性回归示例

先看一个多元线性回归的例子，来自前面的“”。

摘自马里兰州生物资源调查 ()，该调查记录了马里兰州河流中每 75 米水域的鱼种 Rhinichthys cataractae 的丰度，并测量了相应的环境特征。

第一列代表被测河段的位置信息，其余列分别为：

fish，鱼种Rhinichthys cataractae在水中的个体数量，代表物种丰富度；

英亩，流域面积（英亩，英亩）；

do2，水中溶解氧含量（mg/L，mg/L）；

depth，最大水深（cm，cm）；

no3，水中硝酸盐浓度（mg/L，mg/L）；

so4，水中硫酸盐浓度（mg/L，mg/L）；

温度，水温（摄氏度，°C）。

希望通过多元线性回归，探索可能影响白内鲟鱼种丰度的环境因素，并解释鱼种丰度变化的原因。

#读取鱼类物种丰度和水体环境数据
dat <- read.delim('fish_data.txt', sep = 't', row.names = 1)
 
#首先使用全部环境变量拟合与鱼类物种丰度的多元线性回归
fit1 <- lm(fish~acre+do2+depth+no3+so4+temp, data = dat)
summary(fit1)  #展示拟合方程的简单统计

每个结果项的描述见上面的“”。简而言之，英亩（流域面积）和no3（硝酸盐浓度）的增加有助于物种丰度的增加；深度（水深）的增加也对物种丰度的增加略有贡献；没有观察到物种丰度之间的差异。 do2（水中溶解氧含量）、so4（水中硫酸盐浓度）和temp（水中温度）呈线性响应。

忽略全模型后修正后的R2只有0.247，不算很高。这是由于物种丰度与大多数环境因素之间的关系不是线性和简单的。这里选择线性模型。原因只是更容易理解和简化问题。

接下来，一些变量选择方法有望去除一些“无效”的预测因子，即对物种丰度变化贡献低的环境因素，从而优化这个线性模型。

手动选择变量

这点在上一篇文章“”中也有提及，这里就展示一下。

可以根据回归中每个预测变量的p值来判断，不显着的预测变量可以手动去除。例如，do2（水中的溶解氧）、so4（水中的硫酸盐浓度）和 temp（水中的温度）非常不显着，表明这些环境因素的变化不会对物种丰度的增加或减少产生线性响应。可能还有其他非线性响应模式，但由于这里预计只关注线性关系，复杂情况暂不考虑），所以考虑去掉。深度（水深）的临界意义，简单保留，重新建立与显着英亩（流域面积）和no3（硝酸盐浓度）的线性回归。

#只考虑使用 3 个和鱼类物种丰度线性关系较为明显的环境变量拟合多元线性回归
#这里将处于线性关系临界值的 depth 也算在内
fit2 <- lm(fish~acre+no3, data = dat)
summary(fit2)  #展示拟合方程的简单统计
 
#anova() 检验前后两个嵌套模型的拟合优度
#比较去除 3 个不显著的预测变量后，回归整体预测性能是否一样好
anova(fit2, fit1)
 
#AIC 评估前后两个回归复杂性与拟合优度的关系
#AIC 值较小的回归优先选择，表明较少的预测变量已经获得了足够的拟合度
AIC(fit2, fit1)

这是手动选择变量后计算的多元线性回归，比较变量选择前后两次回归的优度。

总体而言，在从回归中移除 3 个无关紧要的预测变量后，它变得精简且易于解释；同时，准确率没有下降。以修正后的R2为例，之前是0.247，这里0.246，是一致的。

最终选定的影响鱼类丰度的环境因素的多元回归公式为：

公式中保留了深度（水的深度），也被认为对物种丰度的贡献很小（虽然不是很显着用r语言做多元线性回归，但也很接近）。你也可以选择去掉它，重新拟合回归，但是这个变量选择存在一些人为的主观性。

公式中的截距项不显着，所以不保留。当其他预测变量为零时，截距项代表物种丰度也是可以解释的。毕竟acre（流域面积）或depth（水深）的值为0时，说明河流已经枯竭，此时鱼类不可能生存，所以截距项的存在为毫无意义。

逐步回归选择变量

当初始回归模型中预测变量较多时，人工判断和选择变量比较麻烦。这时候可以尝试一些机器自动变量选择的方法。

首先介绍逐步法。在逐步回归中，回归模型一次添加或删除一个预测变量，直到达到某个停止标准。逐步回归方法的实现根据添加和删除变量的标准而有所不同。常见的有前向逐步回归、后向逐步回归和双向逐步回归三种。

在R语言中，MASS包中的stepAIC()可以实现逐步回归（前向、后向和双向），基于AIC（Akaike Information Criterion，Akaike Information Criterion）。通过选择最小AIC统计量对应的回归模型来达到删除或添加变量的目的。接下来，利用这种自动选择方法，对上述物种丰度与环境因素关系的多元线性回归进行变量选择。

前向逐步回归

前向逐步回归，或前向选择，一次向回归添加一个预测变量，直到添加变量不会提高回归准确度（总 R2） 直到得到改善。

因此，首先需要建立一个只包含1-2个少数预测变量的起始回归模型（建议找最显着的预测变量作为起始模型），或者直接使用一个没有任何预测变量的空模型，然后使用前向逐步回归，将其他预测变量依次添加到回归中，直到模型达到最优。

library(MASS)
 
#首先构建了一个只包含 1 个非常显著的预测变量的起始回归模型
fit0 <- lm(fish~acre, data = dat)
 
#如果使用不含任何预测变量的空模型
#fit0 <- lm(fish~1, data = dat)
 
#随后应用前向选择逐步往里添加新的预测变量，直到模型最优或达到指定的最大变量数量
stepAIC(fit0, direction = 'forward',
    scope = list(lower = fit0, upper = ~acre+do2+depth+no3+so4+temp))

初始回归模型仅包含一个预测变量 acre（流域面积）。然后在每一步，AIC 列提供模型添加新预测变量后的 AIC 值，其中 AIC 值表示未添加变量时模型的 AIC 值。当添加一个新的预测变量后有一个 AIC 值小于 AIC 值时，选择可以使 AIC 值最小的预测变量添加到原始回归中。以此类推，直到没有添加新的预测器，AIC值小于AIC值终止选择。

在前向逐步回归中选择变量后，最终保留了 3 个有效预测因子，acre（流域面积）、no3（硝酸盐浓度）和 depth（水深）。因此，最终影响鱼类丰度的环境因素多元回归方程为：

与上述结果一致。

向后逐步回归

Backward stepwise regression，或者backward selection，从完整的模型（包括所有的预测变量）中逐一删除预测变量，如果删除后回归准确率下降最少的变量，则执行删除，直到当变量被删除时，回归准确率显着下降已删除。

因此，首先考虑建立一个包含所有预测变量的初始回归模型，然后通过反向逐步回归，依次对回归中的预测变量进行逐一剪枝，直到模型最优。

初始回归模型包括所有 6 个预测变量，acre（流域面积）、no3（硝酸盐浓度）、depth（水深）、do2（水中溶解氧含量）、so4（水中硫酸盐浓度）和 temp（水温） 然后在每一步中，AIC 列提供模型在修剪现有预测变量后的 AIC 值，其中 AIC 值表示未修剪变量时模型的 AIC 值。当删除某个预测变量后存在小于AIC值的AIC值时，在原回归中选择能够使AIC值最小的预测变量并删除。以此类推，直到没有一个AIC值小于删除任意预测变量后的AIC值，才终止选择。

在反向逐步回归中选择变量后，最终保留了 3 个有效预测因子，acre（流域面积）、no3（硝酸盐浓度）和 depth（水深）。因此，最终影响鱼类丰度的环境因素多元回归方程为：

与上述结果一致。

双向逐步回归

逐步逐步回归，或称前向-后向选择，结合了前向选择和后向选择的方法。一次将一个预测变量添加到回归中，并对其进行分析以查看是否可以删除包含非贡献预测变量以改进模型。可以多次添加和删除预测变量，直到获得最佳模型。

因此，首先需要建立一个只包含1-2个少数预测变量的起始回归模型（建议找最显着的预测变量作为起始模型），或者直接使用一个没有任何预测变量的空模型，然后使用双向逐步回归将其他预测变量依次添加到回归中，并考虑是否可以删除现有预测变量之一以提高模型的优度，直到模型达到最优。

#首先构建了一个只包含 1 个非常显著的预测变量的起始回归模型
fit0 <- lm(fish~acre, data = dat)
 
#如果使用不含任何预测变量的空模型
#fit0 <- lm(fish~1, data = dat)
 
#随后执行双向选择，逐步往里添加新的有效预测变量，并考虑删除已存在的无效预测变量，直到模型最优或达到指定的最大变量数量
stepAIC(fit0, direction = 'both',
    scope = list(lower = fit0, upper = ~acre+do2+depth+no3+so4+temp))

初始回归模型仅包含一个预测变量 acre（流域面积）。然后在每一步，AIC 列提供模型在添加或删除预测变量后的 AIC 值，其中 AIC 值表示未添加或删除变量时模型的 AIC 值。当存在AIC值小于添加或删除后的AIC值的预测变量时，在原回归中选择能够使AIC值最小的预测变量进行添加或删除。以此类推，直到没有新的预测变量添加或删除，AIC值小于AIC值终止选择。

在选择变量进行双逐步回归后用r语言做多元线性回归，最终保留了 3 个有效预测因子，acre（流域面积）、no3（硝酸盐浓度）和 depth（水深）。因此，最终影响鱼类丰度的环境因素多元回归方程为：

与上述结果一致。

所有子集回归选择变量

逐步回归方法有一些局限性，例如无法保证每个可能的预测变量组合都参与评估。相比之下，全子集回归能够考虑所有可能的模型组合，综合评价N个不同子集大小（一个、两个或多个预测变量）模型的最佳组合。

在 R 语言中，可以使用leaps 包中的 regsubsets() 来实现全子集回归。接下来，利用这种自动选择方法，对上述物种丰度与环境因素关系的多元线性回归进行变量选择。

library(leaps)
 
#全子集回归评估最佳的 N 变量组合
leap <- regsubsets(fish~acre+do2+depth+no3+so4+temp, data = dat, nbest = 2)
plot(leap, scale = 'adjr2')  #以校正后 R2 代表模型精度

regsubsets() 设置nbest=2，即在每个N个变量组合下的所有情况下，只显示两个最好的N个变量组合。

在这张图中，横轴代表预测变量，纵轴代表回归精度（这里使用的校正后的R2）表示考虑不同预测变量组合时回归精度的变化。例如，只考虑no3作为预测变量时，回归的R2.adj=0.082；只考虑acre作为预测变量时，回归的修正R2=0.11；考虑acre 和 no3 作为预测变量，回归 R2.adj=0.22；当考虑 acre、no3 和 depth 作为预测变量时，回归 R2.adj= 0.25。

从这张图来看，acre（流域面积）、no3（硝酸盐浓度）和 depth（水深）是三个预测变量的最合适的组合。其余情况，要么是修正后的R2过低，要么是预测变量的数量过多，不利于解释。

此外，car包还提供了帮助在全子集回归结果中选择最佳变量子集的方法。

#辅助选择预测变量组合
library(car)
 
subsets(leap, statistic = 'cp')
abline(1,1,lty = 2,col = 'red')

预测变量组合越好，越接近图中的红线（截距项和斜率为 1 的线）。在此基础上，再次判断acre（流域面积）、no3（硝酸盐浓度）和depth（水深）是三个预测因子最合适的组合。

那么，通过这三个环境预测因子，就足以建立与鱼类物种丰度的多元线性回归，不用多说，与上述结果是一致的。

关于变量选择的一些注意事项

对于变量的自动机器选择，全子集回归优于逐步回归，因为它考虑了更多的预测变量组合。但是，当有大量预测变量时，全子集回归可能会很慢。一般来说，自动变量选择应被视为回归模型选择的辅助手段，而不是直接的方法，因为可能会选择数学上有意义但生物学上不重要的模型。

当多元回归涉及的变量较少时，手动选择方法优于机器自动选择。虽然主观性很强，但从经验上建立易于解释的变量关系是有益的。

参考文献

罗伯特·卡巴科夫。 R语言实践（第二版）（王小宁、刘协新、黄俊文等译）。人民邮电出版社，2016.