tmf*_*mnk 11 regression r dplyr tidyverse
我怀疑这个问题可能是重复的，但是，我发现没有什么令人满意的。想象一个具有如下结构的简单数据集：
set.seed(123)
df <- data.frame(cov_a = rbinom(100, 1, prob = 0.5),
cov_b = rbinom(100, 1, prob = 0.5),
cont_a = runif(100),
cont_b = runif(100),
dep = runif(100))
cov_a cov_b cont_a cont_b dep
1 0 1 0.238726027 0.784575267 0.9860542973
2 1 0 0.962358936 0.009429905 0.1370674714
3 0 0 0.601365726 0.779065883 0.9053095817
4 1 1 0.515029727 0.729390652 0.5763018376
5 1 0 0.402573342 0.630131853 0.3954488591
6 0 1 0.880246541 0.480910830 0.4498024841
7 1 1 0.364091865 0.156636851 0.7065019011
8 1 1 0.288239281 0.008215520 0.0825027458
9 1 0 0.170645235 0.452458394 0.3393125802
10 0 0 0.172171746 0.492293329 0.6807875512
我正在寻找的是一个优雅的dplyr/tidyverse选项，可以为每个cov_变量拟合一个单独的回归模型，同时包括相同的一组附加变量和相同的因变量。
我能够使用此代码来解决这个问题（要求purrr，dplyr，tidyr和broom）：
map(.x = names(df)[grepl("cov_", names(df))],
~ df %>%
nest() %>%
mutate(res = map(data, function(y) tidy(lm(dep ~ cont_a + cont_b + !!sym(.x), data = y)))) %>%
unnest(res))
[[1]]
# A tibble: 4 x 6
data term estimate std.error statistic p.value
<list> <chr> <dbl> <dbl> <dbl> <dbl>
1 <tibble [100 × 5]> (Intercept) 0.472 0.0812 5.81 0.0000000799
2 <tibble [100 × 5]> cont_a -0.103 0.0983 -1.05 0.296
3 <tibble [100 × 5]> cont_b 0.172 0.0990 1.74 0.0848
4 <tibble [100 × 5]> cov_a -0.0455 0.0581 -0.783 0.436
[[2]]
# A tibble: 4 x 6
data term estimate std.error statistic p.value
<list> <chr> <dbl> <dbl> <dbl> <dbl>
1 <tibble [100 × 5]> (Intercept) 0.415 0.0787 5.27 0.000000846
2 <tibble [100 × 5]> cont_a -0.0874 0.0984 -0.888 0.377
3 <tibble [100 × 5]> cont_b 0.181 0.0980 1.84 0.0682
4 <tibble [100 × 5]> cov_b 0.0482 0.0576 0.837 0.405
但是，我想避免使用 double-map()并通过使用某种更直接或更优雅的方法来解决它。
我不确定这是否会被认为更直接/优雅，但这是我不使用 double 的解决方案map：
library(tidyverse)
library(broom)
gen_model_expr <- function(var) {
form = paste("dep ~ cont_a + cont_b +", var)
tidy(lm(as.formula(form), data = df))
}
grep("cov_", names(df), value = TRUE) %>%
map(gen_model_expr)
输出（注意不保留数据列）：
[[1]]
# A tibble: 4 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.472 0.0812 5.81 0.0000000799
2 cont_a -0.103 0.0983 -1.05 0.296
3 cont_b 0.172 0.0990 1.74 0.0848
4 cov_a -0.0455 0.0581 -0.783 0.436
[[2]]
# A tibble: 4 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.415 0.0787 5.27 0.000000846
2 cont_a -0.0874 0.0984 -0.888 0.377
3 cont_b 0.181 0.0980 1.84 0.0682
4 cov_b 0.0482 0.0576 0.837 0.405
编辑
为了提高速度性能（受到@TimTeaFan 的启发），下面显示了比较检索协变量名称的不同方法的基准测试。grep("cov_", names(df), value = TRUE)是最快的
# A tibble: 4 x 13
expression min median `itr/sec` mem_alloc
<bch:expr> <bch:> <bch:> <dbl> <bch:byt>
1 names(df)[grepl("cov_", names(df))] 7.59µs 8.4µs 101975. 0B
2 grep("cov_", colnames(df), value = TRUE) 8.21µs 8.96µs 103142. 0B
3 grep("cov_", names(df), value = TRUE) 6.96µs 7.43µs 128694. 0B
4 df %>% select(starts_with("cov_")) %>% colnames 1.17ms 1.33ms 636. 5.39KB
不是直接tidyverse基于而是：可以使用包fixst一次执行多个估计。
语法是明确的，问题的解决方案可以写在一行代码中：
library(fixest)
set.seed(123)
n = 100
df = data.frame(cov_a = rbinom(n, 1, prob = 0.5),
cov_b = rbinom(n, 1, prob = 0.5),
cont_a = runif(n), cont_b = runif(n),
dep = runif(n))
# Estimation: sw means stepwise
res = feols(dep ~ sw(cov_a, cov_b) + cont_a + cont_b, df)
就是这样：两个估计都完成了。（请注意，feols就像lm。）
然后你可以显示结果：
# Display the results
etable(res, order = "Int|cov")
#> model 1 model 2
#> Dependent Var.: dep dep
#>
#> (Intercept) 0.4722*** (0.0812) 0.4148*** (0.0788)
#> cov_a -0.0455 (0.0581)
#> cov_b 0.0482 (0.0576)
#> cont_a -0.1033 (0.0983) -0.0874 (0.0984)
#> cont_b 0.1723. (0.0990) 0.1808. (0.0980)
#> _______________ __________________ __________________
#> S.E. type Standard Standard
#> Observations 100 100
#> R2 0.04823 0.04910
#> Adj. R2 0.01849 0.01938
并轻松地将它们带到 tidyverse：
# Get the results into tidyverse
library(broom)
lapply(as.list(res), function(x) tidy(x))
#> [[1]]
#> # A tibble: 4 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.472 0.0812 5.81 0.0000000799
#> 2 cov_a -0.0455 0.0581 -0.783 0.436
#> 3 cont_a -0.103 0.0983 -1.05 0.296
#> 4 cont_b 0.172 0.0990 1.74 0.0848
#>
#> [[2]]
#> # A tibble: 4 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.415 0.0787 5.27 0.000000846
#> 2 cov_b 0.0482 0.0576 0.837 0.405
#> 3 cont_a -0.0874 0.0984 -0.888 0.377
#> 4 cont_b 0.181 0.0980 1.84 0.0682
这是关于多重估计的文档。您也可以立即对因变量/固定效应（如果适用）/拆分样本应用多个估计——所有这些都非常紧凑。
最后一件事：多重估计是高度优化的，因此与基于循环的解决方案（如地图）相比，您将获得高性能收益。
我提出两种方法：
第一个比较无聊。我们可以使用 {dplyr} 的 rowwise符号代替purrr::map。这种方法有两种风格。在rowwise我们可以使用 (A)mutate %>% unnest或 (B) 之后，我们可以使用group_map. 在这两种方法中，我都避免复制数据，但如果需要，我们可以轻松地将数据包含在每一行中（在设置tibble我们可以做的时候tibble(myvar = ., data = list(df))）。虽然 (A) 给了我们一个包含所有数据的 tibble，但group_map(B) 中的返回一个类似于原始示例中的“双映射”方法的列表。
第二种方法我认为是比较“新鲜”（虽然不太直接的），因为它既不使用rowwise也不是map。这里我们使用{} dplyr的across功能结合cur_column()，为每个输出新列，然后pivot_longer和unnest有一个所有的结果tibble。
最终基准显示：“doulbe_map”最慢（因为重复数据列），其次是“across”和“row_unnest”，而“row_group_map”则相当快。尽管如此，最快的方法是@latlio 使用map和自定义函数（以下称为“map_custom_fun”）的方法，但尽管它正在使用purrr，但它可能不那么“dplyr-ish”。
library(tidyverse)
library(broom)
set.seed(123)
df <- data.frame(cov_a = rbinom(100, 1, prob = 0.5),
cov_b = rbinom(100, 1, prob = 0.5),
cont_a = runif(100),
cont_b = runif(100),
dep = runif(100))
# first approach: using dplyr rowwise
# Variation A:
# rowwise %>% mutate %>% unnest
df %>%
select(starts_with("cov_")) %>%
colnames %>%
tibble(myvar = .) %>%
rowwise() %>%
mutate(res = list(tidy(lm(dep ~ cont_a + cont_b + eval(sym(.data$myvar)), data = df)))) %>%
unnest(res)
#> # A tibble: 8 x 6
#> myvar term estimate std.error statistic p.value
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 cov_a (Intercept) 0.472 0.0812 5.81 0.0000000799
#> 2 cov_a cont_a -0.103 0.0983 -1.05 0.296
#> 3 cov_a cont_b 0.172 0.0990 1.74 0.0848
#> 4 cov_a eval(sym(.data$myvar)) -0.0455 0.0581 -0.783 0.436
#> 5 cov_b (Intercept) 0.415 0.0787 5.27 0.000000846
#> 6 cov_b cont_a -0.0874 0.0984 -0.888 0.377
#> 7 cov_b cont_b 0.181 0.0980 1.84 0.0682
#> 8 cov_b eval(sym(.data$myvar)) 0.0482 0.0576 0.837 0.405
# Variation B:
# rowwise %>% group_map
df %>%
select(starts_with("cov_")) %>%
colnames %>%
tibble(myvar = .) %>%
rowwise() %>%
group_map(.keep = TRUE,
~ tidy(lm(dep ~ cont_a + cont_b + eval(sym(myvar)), data = df)))
#> [[1]]
#> # A tibble: 4 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.472 0.0812 5.81 0.0000000799
#> 2 cont_a -0.103 0.0983 -1.05 0.296
#> 3 cont_b 0.172 0.0990 1.74 0.0848
#> 4 eval(sym(.x$myvar)) -0.0455 0.0581 -0.783 0.436
#>
#> [[2]]
#> # A tibble: 4 x 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 0.415 0.0787 5.27 0.000000846
#> 2 cont_a -0.0874 0.0984 -0.888 0.377
#> 3 cont_b 0.181 0.0980 1.84 0.0682
#> 4 eval(sym(.x$myvar)) 0.0482 0.0576 0.837 0.405
# second approach: using summarise(across)
# we need a `tibble` here, otherwise printing gets messed up
df_tbl <- as_tibble(df)
df_tbl %>%
summarise(across(starts_with("cov_"),
~ list(tidy(lm(
reformulate(c("cont_a","cont_b", cur_column()), "dep"),
data = df_tbl))))) %>%
pivot_longer(cols = everything(),
names_to = "var",
values_to = "res") %>%
unnest(res)
#> # A tibble: 8 x 6
#> var term estimate std.error statistic p.value
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 cov_a (Intercept) 0.472 0.0812 5.81 0.0000000799
#> 2 cov_a cont_a -0.103 0.0983 -1.05 0.296
#> 3 cov_a cont_b 0.172 0.0990 1.74 0.0848
#> 4 cov_a cov_a -0.0455 0.0581 -0.783 0.436
#> 5 cov_b (Intercept) 0.415 0.0787 5.27 0.000000846
#> 6 cov_b cont_a -0.0874 0.0984 -0.888 0.377
#> 7 cov_b cont_b 0.181 0.0980 1.84 0.0682
#> 8 cov_b cov_b 0.0482 0.0576 0.837 0.405
由reprex 包(v0.3.0)于 2021 年 1 月 10 日创建
基准
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> double_map 20.92847 22.238626 23.645280 22.726705 25.002493 34.29351 100
#> row_unnest 15.30179 15.835506 16.714873 16.358134 17.314802 20.60496 100
#> row_groupmap 10.10168 10.490374 11.237398 10.709524 11.452677 20.40186 100
#> across 16.47369 17.117178 18.593908 17.945136 19.431190 29.29384 100
#> map_custom_fun 6.85758 7.311608 7.953066 7.611394 8.305757 19.57006 100