Извлечение значений NA в сводный вывод glm - PullRequest
Новая
       21

Извлечение значений NA в сводный вывод glm

0 голосов
/ 04 марта 2020

Я использую регрессию logisti c для множества различных подмножеств большого фрейма данных. Для этого я использую следующий код (используя dplyr и purrr): 

# define model to be run
mod_fun <- function(df) {
  glm(presence ~ transect, data = df, family = "binomial")
}

# nest data and run model
mod.glm <- dat %>%
  nest(-c(region, fYear, species, road)) %>%
  mutate(model = map(data, mod_fun))

# define functions to extract model coefficients
b_fun <- function(mod) {
  coef(summary(mod))[2]
}
p_fun <- function(mod) {
  coef(summary(mod))[8]
}

# extract coefficients
slope<-mod.glm %>% group_by(species, region, fYear, road) %>%
  transmute(beta = map_dbl(model, b_fun),
            p_val = map_dbl(model, p_fun))

Как видите, я хочу только извлечь оценку и p-значение уклона (называемое transect). Для этого я использую индексирование coef(summary(mod))[2] et c. Проблема в том, что в моем фрейме данных также есть некоторые подмножества, что приводит к переопределенным системам, в которых некоторые коэффициенты будут установлены в NA. Использование coef(summary(mod))[2] извлекает 2-е значение из вывода coef(), и поскольку NA игнорируются в coef(), это больше не будет оценкой transect, которую я хочу извлечь. До сих пор я пытался unsing coef(summary(mod_2), complete = TRUE) (-> ничего не меняется, NA все еще не отображается) и напрямую связываться со значениями coef(summary(mod_2), complete = TRUE)["transect","Estimate"] (-> выдает ошибку). Кто-нибудь знает, как я мог решить эту проблему?

Что я пробовал до сих пор: 

# two example models; mod_2 will result in NAs
mod_1 <- glm(presence ~ transect, data = dat[dat$fYear == 1&  dat$species=="Plantago lanceolata",], family = "binomial")
mod_2 <- glm(presence ~ transect, data = dat[dat$fYear == 2&  dat$species=="Plantago lanceolata",], family = "binomial")

coef(summary(mod_1))[2] # works fine
coef(summary(mod_2))[2] # not the value I want

coef(summary(mod_1), complete = TRUE)["transect","Estimate"] # works fine
coef(summary(mod_2), complete = TRUE)["transect","Estimate"] # error

coef(summary(mod_2), complete = TRUE) # NAs for transect are still not displayed

summary(mod_2)$coefficients["transect","Estimate"] # is not working either

Данные: 

dput(dat)
structure(list(region = c("HWI", "HWI", "HWI", "HWI", "HWI", 
"HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", 
"HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", "HWI", 
"HWI", "HWI", "HWI", "HWI", "HWI", "HWI"), road = c("MK", "MK", 
"MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", 
"MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", "MK", 
"MK", "MK", "MK", "MK", "MK"), transect = c(1L, 1L, 2L, 2L, 3L, 
3L, 4L, 4L, 4L, 4L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 
11L, 11L, 12L, 12L, 13L, 13L, 15L, 15L, 1L), fYear = c(1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L), species = c("Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata", 
"Poa pratensis", "Plantago lanceolata", "Poa pratensis", "Plantago lanceolata"
), presence = c(1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -29L), groups = structure(list(
    fYear = c(1L, 1L, 2L), road = c("MK", "MK", "MK"), species = c("Plantago lanceolata", 
    "Poa pratensis", "Plantago lanceolata"), .rows = list(c(1L, 
    3L, 5L, 7L, 9L, 11L, 13L, 15L, 17L, 19L, 21L, 23L, 25L, 27L
    ), c(2L, 4L, 6L, 8L, 10L, 12L, 14L, 16L, 18L, 20L, 22L, 24L, 
    26L, 28L), 29L)), row.names = c(NA, -3L), class = c("tbl_df", 
"tbl", "data.frame"), .drop = TRUE))

Спасибо за вашу помощь!

Ответы [ 2 ]

2 голосов
/ 06 марта 2020

Переходя на совершенно другой маршрут, вы можете изменить свои функции извлечения, чтобы они возвращали NA, когда они дают ошибки. Это работа для таких функций, как tryCatch(), но я нашел possibly() из purrr , что очень удобно для такого рода задач.

possibly() оборачивает функцию , Аргументы otherwise указывают, какое значение возвращать при возникновении ошибки при использовании функции.

Вот две ваши функции, заключенные в possibly(). Я изменил их, чтобы они специально работали со строкой "transect" в сводке коэффициентов, так что это приведет к ошибке, если ее не существует. 

b_fun <- possibly(
     function(mod) {
          coef(summary(mod))["transect", 1]
          }, otherwise = NA)

p_fun <- possibly(
     function(mod) {
          coef(summary(mod))["transect", 4]
          }, otherwise = NA)

# extract coefficients
mod.glm %>%
     group_by(species, region, fYear, road) %>%
     transmute(beta = map_dbl(model, b_fun),
               p_val = map_dbl(model, p_fun) ) %>%
     ungroup() 

# A tibble: 3 x 6
  species             region fYear road     beta  p_val
  <chr>               <chr>  <int> <chr>   <dbl>  <dbl>
1 Plantago lanceolata HWI        1 MK    -42.7    0.999
2 Poa pratensis       HWI        1 MK     -0.206  0.272
3 Plantago lanceolata HWI        2 MK     NA     NA
2 голосов
/ 06 марта 2020

Я не вижу, как извлечь строки NA из таблицы коэффициентов. Вместо этого ответ может заключаться в добавлении в строки NA после извлечения нужных элементов. Это можно сделать с помощью complete().

В этом примере я использую broom::tidy(), поскольку он позволяет легко получить коэффициенты, меру неопределенности и результаты тестов. Это означает, что я должен отфильтровывать перехваты после факта, но вы наверняка могли бы сделать что-то похожее с вашими функциями, если бы вы добавили столбец к complete() on. 

library(purrr)
library(tidyr)
library(dplyr)

mod.glm %>%
     group_by(species, region, fYear, road) %>%
     transmute(results = map(model, broom::tidy) ) %>%
     unnest(results) %>%
     complete(term = "transect") %>%
     filter(term != "(Intercept)") %>%
     ungroup()

# A tibble: 3 x 9
  species          region fYear road  term    estimate std.error statistic p.value
  <chr>            <chr>  <int> <chr> <chr>      <dbl>     <dbl>     <dbl>   <dbl>
1 Plantago lanceo~ HWI        1 MK    transe~  -42.7   31127.     -0.00137   0.999
2 Plantago lanceo~ HWI        2 MK    transe~   NA        NA      NA        NA    
3 Poa pratensis    HWI        1 MK    transe~   -0.206     0.188  -1.10      0.272
...