问题描述
我怀疑这个问题可能是重复的，但是，我没有发现任何令人满意的东西。想象一个具有如下结构的简单数据集：
set.seed(123)
df <- data.frame(cov_a = rbinom(100,1,prob = 0.5),cov_b = rbinom(100,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() 并通过使用某种更直接或更优雅的方法来解决它。
解决方法
我不确定这是否会被认为更直接/更优雅，但这是我不使用双 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_",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_",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 而是：可以使用包 fixest 一次执行多个估计。
语法是明确的，问题的解决方案可以写在一行代码中：
library(fixest)
set.seed(123)
n = 100
df = data.frame(cov_a = rbinom(n,1,prob = 0.5),cov_b = rbinom(n,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
Here's the documentation 在多重估计上。您也可以立即对因变量/固定效应（如果适用）/拆分样本应用多个估计——所有这些都是紧凑的。
最后一件事：多重估计得到了高度优化，因此相对于基于循环的解决方案（如地图），您将获得高性能提升。
,我提出了两种方法：
第一个比较无聊。我们可以使用 {dplyr} 的 rowwise 符号代替 purrr::map。这种方法有两种风格。在 rowwise 之后，我们可以使用 (A) mutate %>% unnest 或 (B)，我们可以使用 group_map。在这两种方法中，我都避免复制数据，但如果需要，我们可以轻松地将数据包含在每一行中（在设置 tibble 时，我们可以执行 tibble(myvar = .,data = list(df))）。虽然 (A) 为我们提供了一个包含所有数据的 tibble，但 (B) 中的 group_map 返回一个类似于原始示例中的“双映射”方法的列表。
我认为第二种方法相当“新鲜”（虽然不那么直接），因为它既不使用 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,cov_b = rbinom(100,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 package (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
这是另一个版本：
library(broom)
purrr::map_df(grep("cov_",value = TRUE),~tidy(lm(reformulate(c('cont_a','cont_b',.x),'dep'),data = df)))
# 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
#5 (Intercept) 0.415 0.0787 5.27 0.000000846
#6 cont_a -0.0874 0.0984 -0.888 0.377
#7 cont_b 0.181 0.0980 1.84 0.0682
#8 cov_b 0.0482 0.0576 0.837 0.405
如果您想要返回列表输出，请使用 map 而不是 map_df。
再补充一个。如果您的所有协变量都相同并且您只想更改 cov_ 变量，您可以通过先选择变量然后在 ~. 中使用 lm() 表示法来简化它。正如您可能猜到的，这意味着“包括数据集中所有不是 dv 的其他变量”。
purrr::map(grep("cov_",function(x){
df2 <- select(df,x,cont_a,cont_b,dep)
tidy(lm(dep ~ .,df2))
})
有人用 data.table 尝试过吗？
library("data.table")
covariables <- c("cov_a","cov_b")
# or,according to what you want to do :
covariables <- names(df)[grepl(names(df),pattern = "cov_")]
formulas <- paste0("dep ~ cont_a + cont_b + ",covariables )
res <- df[,lapply(formulas,function(x) list(lm(x,data=.SD)))]
summary.lm(res$V1[[1]])
summary.lm(res$V2[[1]])
这是一个尝试（在之前的答案中这里和那里收集了其他东西以进行一些比较）（请注意，library broom 和 tidyverse 仅用于比较）：
library(tidyverse)
library(broom)
library("data.table")
library(microbenchmark)
test <- microbenchmark(dt = {
covariables <- c("cov_a","cov_b")
covariables <- names(df)[grepl(names(df),data=.SD)))]
summary.lm(res$V1[[1]])
summary.lm(res$V2[[1]])
},broom = {
df_tbl <- as_tibble(df)
df_tbl %>%
summarise(across(starts_with("cov_"),values_to = "res") %>%
unnest(res)
},gep_cov = {
gen_model_expr <- function(var) {
form = paste("dep ~ cont_a + cont_b +",var)
tidy(lm(as.formula(form),data = df))
}
grep("cov_",value = TRUE) %>%
map(gen_model_expr)
},times = 100)
结果：
> test
Unit: milliseconds
expr min lq mean median uq max neval cld
dt 2.3403 2.63330 3.116303 2.8403 3.21620 10.5219 100 a
broom 18.5069 20.40585 22.711605 21.3388 24.04675 46.8398 100 c
gep_cov 7.7221 8.39180 10.086074 9.4315 10.54345 21.0377 100 b