时间戳

首页 > 解决方案 > 使用 purrr、broom 和 dplyr1.0.0 中的几个不同预测器将 glm 应用于相同的结果变量

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

问题描述

在这篇purrr文章中，我展示了如何使用和broom包中的功能使用相同的预测器对多个结果运行多个 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 将返回的内容

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

但我想一次完成所有这些，并在单个对象中返回结果，使用purrand broom。我尝试了以下语法，

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 的截距和系数。

非常感谢任何帮助

有人可以

标签： rdplyrtidyversepurrrbroom

解决方案

如果你想坚持你的方法，我建议你使用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

推荐阅读