0 r combinatorics logistic-regression
我想进行逻辑回归:我有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²)对它们进行排序,以便向我的" 学生 "展示哪些预测器似乎有趣但没有科学意义.我打算执行bootstrap重采样以进一步测试"最狂野 "的模型.
我不确定这种"教育练习"的价值,但为了编程,这将是我的方法:
首先,让我们创建一些示例预测变量名称.我在你的例子中使用了5个预测变量,但是对于10,你显然需要用10替换5.
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 = "+")))}))
或等价使用FUN的参数combn(如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平方,sigma,F统计量,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
你可以按照自己的意愿排序.