在 purrr、broom 和 dplyr1.0.0 中使用几个不同的预测因子将 glm 应用于相同的结果变量

applying a glm to the same outcome variable using several different predictors in purrr, broom, and dplyr1.0.0

post 中,我看到了如何使用 purrrbroom 中的功能,使用相同的预测器对多个结果进行 运行 多个 glms包裹。现在我想反其道而行之,对具有不同预测变量但使用相同结果(即一系列单变量测试)的多个模型分别应用一个 glm

# data
set.seed(1234)
df <- data.frame(out = c(rbinom(50, 1, prob = 0.2),
                          rbinom(50, 1, prob = 0.8)),
                 pred1 = factor(rep(letters[1:2], each = 50)),
                 pred2 = factor(rep(letters[1:2], times = 50)))

这就是每个 glm return 如果我们 运行 他们孤立

summary(mod1 <- glm(out ~ pred1, data = df, family = binomial))  

# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  -1.6582     0.3858  -4.299 1.72e-05 ***
# pred1b        3.4735     0.5612   6.190 6.03e-10 ***


summary(mod2 <- glm(out ~ pred2, data = df, family = binomial))

# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)  0.08004    0.28307   0.283    0.777
# pred2b      -0.08004    0.40016  -0.200    0.841

但我想一次完成这一切,并使用 purrbroom 将结果 return 编辑到一个对象中。我尝试了以下语法,

df %>% 
  select(starts_with("pred")) %>%
    map_df(~broom::tidy(glm(out ~ .,
                            data = df,
                            family = binomial)))

# output  
#   term        estimate std.error statistic  p.value
#   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
# 1 (Intercept)   -1.58      0.469    -3.38  7.36e- 4
# 2 pred1b         3.48      0.562     6.19  6.19e-10
# 3 pred2b        -0.157     0.561    -0.280 7.79e- 1
# 4 (Intercept)   -1.58      0.469    -3.38  7.36e- 4
# 5 pred1b         3.48      0.562     6.19  6.19e-10
# 6 pred2b        -0.157     0.561    -0.280 7.79e- 1

现在这个输出是我想要的格式(即数据帧),但结果本身很奇怪,一些系数(例如pred1)报告正确但报告了两次,一些报告两次但不正确(例如模型 1 的截距),还有一些被忽略(例如模型 2 的截距和系数。

非常感谢任何帮助

有人可以

如果你想坚持你的方法,我建议你使用 map2 而不是 map:

library(dplyr)
library(purrr)
library(broom)

df %>%
  select(starts_with("pred")) %>%
  map2_dfr(names(df)[-1], ~ tidy(glm(out ~ .x, data = df, family = "binomial")) %>%
            mutate(term = c("(Intercept)", .y)))

# A tibble: 4 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  -1.66       0.386    -4.30  1.72e- 5
2 pred1         3.47       0.561     6.19  6.03e-10
3 (Intercept)   0.0800     0.283     0.283 7.77e- 1
4 pred2        -0.0800     0.400    -0.200 8.41e- 1

这可能是另一种有点冗长但有效的方法:

library(purrr)
library(broom)

fn <- function(n) {
  tidy(glm(df[["out"]] ~ df[[n]], data = df, family = "binomial")) %>%
    mutate(term = c("(Intercept)", names(df)[n]))
}

seq_len(ncol(df))[-c(1, 2)] %>%
  reduce(~ .x %>% 
           bind_rows(fn(.y)), .init = fn(seq_len(ncol(df))[2]))

# A tibble: 4 x 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)  -1.66       0.386    -4.30  1.72e- 5
2 pred1         3.47       0.561     6.19  6.03e-10
3 (Intercept)   0.0800     0.283     0.283 7.77e- 1
4 pred2        -0.0800     0.400    -0.200 8.41e- 1