逻辑回归：如何尝试 R 中的每种预测变量组合？答案

【问题标题】：Logistic regression: how to try every combination of predictors in R?逻辑回归：如何尝试 R 中的每种预测变量组合？
【发布时间】：2017-07-23 22:51:47
【问题描述】：

这是https://stats.stackexchange.com/questions/293988/logistic-regression-how-to-try-every-combination-of-predictors 的副本。

我想执行逻辑回归：我有 1 个因变量和约 10 个预测变量。

我想对每个组合进行详尽的搜索，例如更改顺序和添加/删除预测变量等。例如：

y ~ x1 + x2 + x3 + x4 + x5
y ~ x2 + x1 + x3 + x4 + x5
y ~ x1 + x2 + x3
y ~ x5 + x1 + x2 + x3 + x4
y ~ x4 + x2
...等等。

在这种情况下，计算时间对我来说并不是一个停止的问题：这主要是一种教育练习。

你知道我该怎么做吗？我用 R。

编辑：要明确：这主要是一个教育练习：我想测试每个模型，以便我可以根据一些索引（例如 AUC 或伪 R²）按顺序对它们进行排序向我的“学生”展示哪些预测变量看起来很有趣，但没有科学意义。我计划执行引导重采样以进一步测试“fishiest”模型。

【问题讨论】：

订单在这里不重要吗？
正如其他人在您在 stackexchange 中的帖子中评论的那样，y ~ x1 + x2 + x3 + x4 + x5 与 y ~ x2 + x1 + x3 + x4 + x5 相同
谢谢，但这主要是一个教育练习：即使顺序无关紧要，我需要尝试每个组合
为什么您需要尝试每个订单......您的编辑也没有解释它，因为无论订单如何，AUC 和其他值都是相同的
如果顺序无关紧要，那么从 10 个变量开始，您可以拟合 2^10 = 1024 个可能的模型。如果顺序确实很重要，那么您还有...更多。

标签： r combinatorics logistic-regression

【解决方案1】：

我不确定这个“教育练习”的价值，但为了编程，这是我的方法：

首先，让我们创建一些示例预测器名称。我在您的示例中使用 5 个预测变量，但对于 10 个，您显然需要将 5 替换为 10。

X = paste0("x",1:5)
X
[1] "x1" "x2" "x3" "x4" "x5"

现在，我们可以得到combn的组合。

例如，一次只有一个变量：

 t(combn(X,1))
     [,1]
[1,] "x1"
[2,] "x2"
[3,] "x3"
[4,] "x4"
[5,] "x5"

一次两个变量：

> t(combn(X,2))
      [,1] [,2]
 [1,] "x1" "x2"
 [2,] "x1" "x3"
 [3,] "x1" "x4"
 [4,] "x1" "x5"
 [5,] "x2" "x3"
 [6,] "x2" "x4"
 [7,] "x2" "x5"
 [8,] "x3" "x4"
 [9,] "x3" "x5"
[10,] "x4" "x5"

等等

我们可以使用lapply 连续调用这些函数，同时考虑越来越多的变量，并将结果捕获到一个列表中。例如，看看lapply(1:5, function(n) t(combn(X,n))) 的输出。要将这些组合转化为公式，我们可以使用以下方法：

out <- unlist(lapply(1:5, function(n) {
  # get combinations
  combinations <- t(combn(X,n))
  # collapse them into usable formulas:
  formulas <- apply(combinations, 1, 
                    function(row) paste0("y ~ ", paste0(row, collapse = "+")))}))

或等效地使用combn 的FUN 参数（正如user20650 所指出的那样）：

out <- unlist(lapply(1:5, function(n) combn(X, n, FUN=function(row) paste0("y ~ ", paste0(row, collapse = "+")))))

这给出了：

out
 [1] "y ~ x1"             "y ~ x2"             "y ~ x3"             "y ~ x4"             "y ~ x5"            
 [6] "y ~ x1+x2"          "y ~ x1+x3"          "y ~ x1+x4"          "y ~ x1+x5"          "y ~ x2+x3"         
[11] "y ~ x2+x4"          "y ~ x2+x5"          "y ~ x3+x4"          "y ~ x3+x5"          "y ~ x4+x5"         
[16] "y ~ x1+x2+x3"       "y ~ x1+x2+x4"       "y ~ x1+x2+x5"       "y ~ x1+x3+x4"       "y ~ x1+x3+x5"      
[21] "y ~ x1+x4+x5"       "y ~ x2+x3+x4"       "y ~ x2+x3+x5"       "y ~ x2+x4+x5"       "y ~ x3+x4+x5"      
[26] "y ~ x1+x2+x3+x4"    "y ~ x1+x2+x3+x5"    "y ~ x1+x2+x4+x5"    "y ~ x1+x3+x4+x5"    "y ~ x2+x3+x4+x5"   
[31] "y ~ x1+x2+x3+x4+x5"

现在可以将其传递给您的逻辑回归函数。

例子：

让我们使用mtcars 数据集，以mpg 作为因变量。

X = names(mtcars[,-1])
X
 [1] "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear" "carb"

现在，让我们使用前面提到的函数：

out <- unlist(lapply(1:length(X), function(n) combn(X, n, FUN=function(row) paste0("mpg ~ ", paste0(row, collapse = "+")))))

它为我们提供了所有组合的向量作为公式。

要运行相应的模型，我们可以这样做

mods = lapply(out, function(frml) lm(frml, data=mtcars))

由于您想捕获特定的统计数据并相应地对模型进行排序，我会使用broom::glance。 broom::tidy 将 lm 输出转换为数据帧（如果您想比较系数等很有用），broom::glance 将例如r-squared、sigma、F-statistic、logLikelihood、AIC、BIC 等到数据框中。例如：

library(broom)
library(dplyr)
tmp = bind_rows(lapply(out, function(frml) {
  a = glance(lm(frml, data=mtcars))
  a$frml = frml
  return(a)
}))

head(tmp)
  r.squared adj.r.squared    sigma statistic      p.value df    logLik      AIC      BIC deviance df.residual       frml
1 0.7261800     0.7170527 3.205902 79.561028 6.112687e-10  2 -81.65321 169.3064 173.7036 308.3342          30  mpg ~ cyl
2 0.7183433     0.7089548 3.251454 76.512660 9.380327e-10  2 -82.10469 170.2094 174.6066 317.1587          30 mpg ~ disp
3 0.6024373     0.5891853 3.862962 45.459803 1.787835e-07  2 -87.61931 181.2386 185.6358 447.6743          30   mpg ~ hp
4 0.4639952     0.4461283 4.485409 25.969645 1.776240e-05  2 -92.39996 190.7999 195.1971 603.5667          30 mpg ~ drat
5 0.7528328     0.7445939 3.045882 91.375325 1.293959e-10  2 -80.01471 166.0294 170.4266 278.3219          30   mpg ~ wt
6 0.1752963     0.1478062 5.563738  6.376702 1.708199e-02  2 -99.29406 204.5881 208.9853 928.6553          30 mpg ~ qsec

您可以根据需要对其进行排序。

【讨论】：

你为什么要加载combinat？ combn 包含在 utils 包中。
只是为了小伙伴们，看看如果你也包括变量的所有不同顺序会发生什么。
好吧，我相信 OP 最终会意识到这对于给定的任务并不是很有洞察力。
查看：rstudio-pubs-static.s3.amazonaws.com/…，包括 glmulti 的链接：jstatsoft.org/v34/i12/paper
谢谢。即使没有重复，这就是我的意思。

【解决方案2】：

有一个包可以做到这一点，MuMIn（multimodel inference），作为更有原则的 multi -model 方法（即它不只是选择最好的模型并忽略选择已经完成的事实）：

设置数据和完整模型：

set.seed(101)
d <- data.frame(replicate(5,rnorm(100)))
d$y <- rbinom(100,size=1,prob=0.5)
full <- glm(y~.,data=d,na.action=na.fail)

“疏通”结果：

library(MuMIn)
allfits <- dredge(full)

结果（也包含所有拟合参数）：

head(allfits[,7:11])
##    df    logLik     AICc    delta     weight
## 3   3 -69.66403 145.5781 0.000000 0.15916685
## 11  4 -69.22909 146.8792 1.301191 0.08304293
## 19  4 -69.30856 147.0382 1.460123 0.07669921
## 7   4 -69.31233 147.0457 1.467655 0.07641093
## 4   4 -69.40589 147.2328 1.654775 0.06958615
## 1   2 -72.07662 148.2769 2.698896 0.04128523

【讨论】：

谢谢，非常有趣的包