为什么 R 在使用 group_nest() 时忽略 relevel？

Question

作为这个的延续，我正在尝试有效地执行许多逻辑回归，以便生成一个列来说明某个组是否与我的参考组有显着差异。

当我尝试将数据嵌套在一列中时，效果很好。但是，现在我需要按两列分组，代码运行了，但我无法更改引用组。我尝试了以下方法：

添加 relevel 参数（如下所示）
在自定义函数本身中添加 relevel 参数（如下所示）
将所需的引用组重命名为以 'AAA' 开头以欺骗 R 使其成为第一个选项

这是一个示例数据集：

library(dplyr)
library(lubridate)
library(tidyr)
library(purrr)
library(broom)

test <- tibble(
  major = as.factor(c(rep(c("undeclared", "computer science", "english"), 2), "undeclared")),
  app_deadline = ymd(c(rep("'2021-04-04", 3), rep("'2020-03-23", 3), rep("'2019-05-23", 1))),
  time = ymd(c(rep("'2021-01-01", 3), rep("'2020-01-01", 3), rep("'2019-01-01", 1))),
  admit = c(500, 1000, 450, 800, 300, 100, 1000),
  reject = c(1000, 300, 1000, 210, 100, 900, 1500)
)

test2 <- test %>%
  mutate(total = rowSums(test[ , c("admit", "reject")], na.rm=TRUE)) %>%
  mutate(accept_rate = admit/total)

这是不允许我更改参考级别的代码：

#Custom function --note that english has been set as reference level
library(tidyr)
library(dplyr)
library(purrr)
library(broom)


get_model_t <- function(df) {
  tryCatch(
    expr = glm(accept_rate ~ relevel(major, ref = "english"), data = df, family = binomial, weights = total, na.action = na.exclude),
    error = function(e) NULL, warning=function(w) NULL)
}

#putting it altogether--note again that english has been marked as reference level
test2 %>% 
  # create year column
  mutate(year = year(time), 
         major = relevel(major, "english")) %>%
  
  # nest by year
  group_nest(year, app_deadline) %>% 
  
  # compute regression
  mutate(reg = map(data, get_model_t), reg_tidy = map(reg, tidy)) %>% 
  
  
  # get data and regression results back to tibble form
  unnest(c(data, reg_tidy)) %>% 
  filter(term != "(Intercept)") %>%
  
  # create the significant yes/no column
  mutate(significant = ifelse(p.value < 0.05, "Yes", "No")) %>% 
  
  # remove the unnecessary columns
  select(-c(term, estimate, std.error, statistic, p.value, reg)) %>%
  full_join(test2)

#Note that, based on the significance column, it's clear that 'undeclared' is being used as the reference group

为什么会这样？对于解决方案，我更喜欢它是否可以灵活——即，不仅适用于 'english'，而且也可以切换为适用于 'computer science'。

Answer 1

它确实尊重 relevel() 函数，问题是返回的结果与 major 列不匹配。看看如果你停在 unnest() 函数会发生什么：

test2 <- test %>%
  mutate(total = rowSums(test[ , c("admit", "reject")], na.rm=TRUE)) %>%
  mutate(accept_rate = admit/total)


get_model_t <- function(df) {
  tryCatch(
    expr = glm(accept_rate ~ relevel(major, ref = "english"), data = df, family = binomial, weights = total, na.action = na.exclude),
    error = function(e) NULL, warning=function(w) NULL)
}

#putting it altogether--note again that english has been marked as reference level
tmp <- test2 %>% 
  # create year column
  mutate(year = year(time), 
         major = relevel(major, "english")) %>%
  
  # nest by year
  group_nest(year, app_deadline) %>% 
  
  # compute regression
  mutate(reg = map(data, get_model_t), reg_tidy = map(reg, tidy)) %>% 
  
  
  # get data and regression results back to tibble form
  unnest(c(data, reg_tidy))

现在，看看 major 和 term

tmp %>% select(major, term)
# # A tibble: 6 × 2
#   major            term                                               
#   <fct>            <chr>                                              
# 1 undeclared       "(Intercept)"                                      
# 2 computer science "relevel(major, ref = \"english\")computer science"
# 3 english          "relevel(major, ref = \"english\")undeclared"      
# 4 undeclared       "(Intercept)"                                      
# 5 computer science "relevel(major, ref = \"english\")computer science"
# 6 english          "relevel(major, ref = \"english\")undeclared"

可以看到major为"english"的那几行，其实是针对"undeclared"参数估计的。根据以上结果，我认为您可以通过以下方式捕获您想要的内容：

tmp %>% 
  filter(term != "(Intercept)") %>% 
  mutate(major = gsub(".*\)(.*)", "\1", term)) %>% 
  # create the significant yes/no column
  mutate(significant = ifelse(p.value < 0.05, "Yes", "No")) %>% 
  
  # remove the unnecessary columns
  select(year, app_deadline, major, time, significant) %>%
  full_join(test2)

# Joining, by = c("app_deadline", "major", "time")
# # A tibble: 7 × 9
#   year app_deadline major            time       significant admit reject total accept_rate
#   <dbl> <date>       <chr>            <date>     <chr>       <dbl>  <dbl> <dbl>       <dbl>
# 1  2020 2020-03-23   computer science 2020-01-01 Yes           300    100   400       0.75 
# 2  2020 2020-03-23   undeclared       2020-01-01 Yes           800    210  1010       0.792
# 3  2021 2021-04-04   computer science 2021-01-01 Yes          1000    300  1300       0.769
# 4  2021 2021-04-04   undeclared       2021-01-01 No            500   1000  1500       0.333
# 5    NA 2021-04-04   english          2021-01-01 NA            450   1000  1450       0.310
# 6    NA 2020-03-23   english          2020-01-01 NA            100    900  1000       0.1  
# 7    NA 2019-05-23   undeclared       2019-01-01 NA           1000   1500  2500       0.4

为什么 R 在使用 group_nest() 时忽略 relevel？

Why does R ignore relevel when using group_nest()?

r

nested

tidyr

purrr

broom