В качестве продолжения этого вопрос я хочу запустить сразу несколько уравнений логистической регрессии, а затем отметить, значительно ли группа отличается от контрольной группы. Это решение работает, но оно работает только тогда, когда я не пропускаю значения. Поскольку мои данные состоят из 100 уравнений, они обязательно будут иметь пропущенные значения, поэтому вместо того, чтобы это решение терпело неудачу при обнаружении ошибки, как я могу запрограммировать его на пропуск экземпляров, которые вызывают ошибку?
Вот измененный набор данных, в котором отсутствуют случаи:
library(dplyr)
library(lubridate)
library(broom)
test <- tibble(major = as.factor(c(rep(c("undeclared", "computer science", "english"), 2), "undeclared")),
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)) %>%
mutate(total = rowSums(test[ , c("admit", "reject")], na.rm=TRUE)) %>%
mutate(accept_rate = admit/total)
А вот решение, которое работает, когда в нем есть все случаи (см. набор данных здесь), но когда оно попадает в группу 2019, в которой отсутствуют случаи, оно терпит неудачу:
library(dplyr)
library(lubridate)
library(broom)
library(tidyr)
library(purrr)
test %>%
# create year column
mutate(year = year(time),
major = relevel(major, "undeclared")) %>%
# nest by year
nest(data = -year) %>%
# compute regression
mutate(reg = map(data, ~glm(accept_rate ~ major, data = .,
family = binomial, weights = total, na.action = na.exclude)),
# use broom::tidy to make a tibble out of model object
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))
Я также попытался обернуть решение с помощью пользовательских функций здесь, но я также не смог заставить его работать. Наконец, я также открыт для других идей решения, если оно дает аналогичный результат и устойчиво к этим ошибкам.
purrr::safely
позволяет позаботиться об ошибках. Чтобы обернуть вызов glm
внутри purrr::safely
, я использую вспомогательную функцию glm_safe
. glm_safe
возвращает список с двумя элементами, result
и error
.
Когда ошибки нет, result
содержит объект модели, а element
— это NULL
. В случае ошибки сообщение об ошибке сохраняется в error
, а result
равно NULL
.
Чтобы использовать результаты в вашем конвейере, мы должны извлечь элементы result
, которые можно получить с помощью transpose(reg)$result
.
library(dplyr)
library(lubridate)
library(broom)
library(tidyr)
library(purrr)
test <- tibble(
major = as.factor(c(rep(c("undeclared", "computer science", "english"), 2), "undeclared")),
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)
)
test <- test %>%
mutate(total = rowSums(test[, c("admit", "reject")], na.rm = TRUE)) %>%
mutate(accept_rate = admit / total)
glm_safe <- purrr::safely(
function(x) {
glm(accept_rate ~ major,
data = x,
family = binomial, weights = total, na.action = na.exclude
)
}
)
test %>%
# create year column
mutate(
year = year(time),
major = relevel(major, "undeclared")
) %>%
# nest by year
nest(data = -year) %>%
# compute regression
mutate(reg = map(data, glm_safe),
reg = transpose(reg)$result) |>
mutate(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))
#> # A tibble: 4 × 8
#> year major time admit reject total accept_rate significant
#> <dbl> <fct> <date> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 2021 computer science 2021-01-01 1000 300 1300 0.769 Yes
#> 2 2021 english 2021-01-01 450 1000 1450 0.310 No
#> 3 2020 computer science 2020-01-01 300 100 400 0.75 No
#> 4 2020 english 2020-01-01 100 900 1000 0.1 Yes
Хм. Я только что внес изменения, чтобы добавить полный репрекс, включая данные. Который, как вы могли видеть, работает без ошибок. Может быть, вы могли бы попробовать еще раз.
да, он все еще дает мне ошибку даже после обновления. Возможно, это что-то с моим компьютером. Тем не менее, спасибо, что нашли время, чтобы представить это!
Чтобы игнорировать ошибки, используйте эту функцию:
get_model <- function(df) {
tryCatch(
glm(accept_rate ~ major, data = df, family = binomial, weights = total, na.action = na.exclude),
error = function(e) NULL, warning=function(w) NULL)
}
Используйте его там, где вы вызываете mutate(reg=map()...)
:
# compute regression
mutate(reg = map(data, get_model), reg_tidy = map(reg, tidy))
Выход:
# A tibble: 4 x 8
year major time admit reject total accept_rate significant
<dbl> <fct> <date> <dbl> <dbl> <dbl> <dbl> <chr>
1 2021 computer science 2021-01-01 1000 300 1300 0.769 Yes
2 2021 english 2021-01-01 450 1000 1450 0.310 No
3 2020 computer science 2020-01-01 300 100 400 0.75 No
4 2020 english 2020-01-01 100 900 1000 0.1 Yes
Спасибо за вашу помощь! Когда я пытаюсь запустить это, я получаю сообщение об ошибке: Ошибка: Проблема со столбцом
mutate()
reg
. ℹreg = transpose(reg)$result
. x Элемент 1 ввода списка не является атомарным вектором