逻辑回归：如何尝试 R 中预测变量的每个组合？
Logistic regression: how to try every combination of predictors in R?
我想执行逻辑回归：我有 1 个因变量和 ~10 个预测变量。
我想尝试每个组合执行详尽搜索，例如更改顺序和 adding/deleting 预测变量等。例如：
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²) 以便向我的“students”展示哪些预测变量看起来很有趣但没有科学意义。我计划执行 bootstrap 重采样以进一步测试“fishiest”模型。
我不确定这个 "educational exercise" 的价值，但为了编程，我的做法是：
首先，让我们创建一些示例预测变量名称。我在您的示例中使用了 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, the F-statistic, the logLikelihood, AIC, BIC etc into a dataframe.例如：
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
您可以随意排序。
有一个包可以做到这一点，MuMIn (multimodel in参考），作为更有原则的多模型方法的一部分（即它不只是选择最好的模型而忽略已经完成选择的事实）：
设置数据和完整模型：
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)
"dredge" 结果：
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
我想执行逻辑回归：我有 1 个因变量和 ~10 个预测变量。
我想尝试每个组合执行详尽搜索，例如更改顺序和 adding/deleting 预测变量等。例如：
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²) 以便向我的“students”展示哪些预测变量看起来很有趣但没有科学意义。我计划执行 bootstrap 重采样以进一步测试“fishiest”模型。
我不确定这个 "educational exercise" 的价值，但为了编程，我的做法是：
首先，让我们创建一些示例预测变量名称。我在您的示例中使用了 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, the F-statistic, the logLikelihood, AIC, BIC etc into a dataframe.例如：
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
您可以随意排序。
有一个包可以做到这一点，MuMIn (multimodel in参考），作为更有原则的多模型方法的一部分（即它不只是选择最好的模型而忽略已经完成选择的事实）：
设置数据和完整模型：
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)
"dredge" 结果：
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