在 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 中,我看到了如何使用 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 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
但我想一次完成这一切,并使用 purr
和 broom
将结果 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
在 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 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
但我想一次完成这一切,并使用 purr
和 broom
将结果 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