最近有个小伙伴问到了这样的一个问题，他想对所有自变量的组合进行回归，然后从中筛选系数估计值均显著的。
为此我们使用 auto 数据集进行演示：
library(tidyverse)
library(leaps)
library(modelr)
haven::read_dta("auto.dta") -> auto
|
假如我们想考虑这些变量：headroom trunk length displacement 对 mpg 的影响：
auto %>%
select(mpg, headroom, trunk, length, displacement) -> df
df
|
leaps 包可以实现选择最优模型的操作，不过结果不是很符合我们的需求，但是 leaps 包的 leaps() 函数可以返回所有可能的自变量组合，这样使用即可：
library(leaps)
df %>%
select(headroom, trunk, length, displacement) %>%
as.matrix() -> xm
df %>%
select(mpg) %>%
pull() -> y
leaps(x = xm, y = y) -> leapsdf
leapsdf$which -> regm
head(regm)
|
然后依次循环所有的 15 个模型：
regm %>%
as_tibble() %>%
set_names("headroom", "trunk", "length", "displacement") %>%
mutate(m = row_number()) %>%
mutate(lmfit = map(m, ~lm(y ~ xm[,regm[.x,]]))) -> regdf
|
broom 包的 tidy() 函数可以自动从模型中提取整洁的估计结果：
regdf %>%
mutate(lmfit = map(lmfit, broom::tidy)) %>%
unnest(lmfit) %>%
filter(term != "(Intercept)") %>%
group_by(m) %>%
mutate(pvalmax = max(p.value)) %>%
ungroup() %>%
filter(pvalmax <= 0.05) %>%
select(m, everything())
|
这样我们就得到了我们想要的结果。
变量多一些也不是问题。首先我们生成一些新的随机变量：
df %>%
mutate_at(c("headroom", "trunk", "length", "displacement"),
funs(`2` = . * runif(n = 74))) %>%
mutate_at(c("headroom", "trunk", "length", "displacement"),
funs(`3` = . * runif(n = 74))) -> df2
df2
#> # A tibble: 74 × 13
#> mpg headroom trunk length displacement headroom_2 trunk_2 length_2
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 22 2.5 11 186 121 1.54 2.07 178.
#> 2 17 3 11 173 258 1.35 4.11 152.
#> 3 22 3 12 168 121 2.25 10.2 1.19
#> 4 20 4.5 16 196 196 2.72 11.8 110.
#> 5 15 4 20 222 350 0.650 1.31 101.
#> 6 18 4 21 218 231 2.28 19.8 44.5
#> 7 26 3 10 170 304 2.51 5.08 20.9
#> 8 20 2 16 200 196 1.42 7.27 36.6
#> 9 16 3.5 17 207 231 3.44 0.790 128.
#> 10 19 3.5 13 200 231 2.00 7.18 150.
#> # ℹ 64 more rows
#> # ℹ 5 more variables: displacement_2 <dbl>, headroom_3 <dbl>, trunk_3 <dbl>,
#> # length_3 <dbl>, displacement_3 <dbl>
|
现在总共有这么些变量：
df2 %>%
colnames() %>%
dput()
|
为了不至于组合过多，选择这些变量：”mpg”, “headroom”, “trunk”, “length”, “displacement”, “headroom_2”, “trunk_2”, “length_2”, “displacement_2”, “headroom_3”：
df2 %>%
select("mpg", "headroom", "trunk", "length",
"displacement", "headroom_2",
"trunk_2", "length_2",
"displacement_2", "headroom_3") -> df2
|
获取所有的组合：
df2 %>%
select("headroom", "trunk", "length",
"displacement", "headroom_2",
"trunk_2", "length_2",
"displacement_2", "headroom_3") %>%
as.matrix() -> xm
df2 %>%
select(mpg) %>%
pull() -> y
leaps(x = xm, y = y) -> leapsdf
leapsdf$which -> regm
head(regm)
|
循环回归所有的组合并保留符合要求的：
regm %>%
as_tibble() %>%
set_names("headroom", "trunk", "length",
"displacement", "headroom_2",
"trunk_2", "length_2",
"displacement_2", "headroom_3") %>%
mutate(m = row_number()) %>%
mutate(lmfit = map(m, ~lm(y ~ xm[,regm[.x,]]))) -> regdf
library(broom)
regdf %>%
mutate(lmfit = map(lmfit, broom::tidy)) %>%
unnest(lmfit) %>%
filter(term != "(Intercept)") %>%
group_by(m) %>%
mutate(pvalmax = max(p.value)) %>%
ungroup() %>%
filter(pvalmax <= 0.05) %>%
select(m, everything())
|
这样我们就解决了这个问题。
点击这里跳转到 RStata 短书平台获取附件：如何使用 R 语言生成所有自变量的组合进行回归并筛选自变量系数估计值均显著的模型？
评论