未嵌套的 glm 模型

Question

我对嵌套的 glm 模型很感兴趣。我嵌套在一个变量 (region) 和运行一个适合模型的函数 region_model 上。

# toy data
test_data = data.frame(region = sample(letters[1:3], 1000, replace = TRUE),
              x = sample(0:1, 1000, replace = TRUE), 
                                   y = sample(1:100, 1000, replace = TRUE), 
                                   z = sample(0:1, 1000, replace = TRUE)) %>% arrange(region)

# nest
by_region = test_data %>%
              group_by(region) %>%
              nest()


# glm function 
region_model  <- function(df) {
 glm(x ~ y + z, data = df, family = "binomial")
}              

# run the model  
    by_region = by_region %>% mutate(mod_rat = data %>% map(region_model))

生成的小标题如下所示：

> by_region
# A tibble: 3 x 3
  region               data   mod_rat
  <fctr>             <list>    <list>
1      a <tibble [352 x 3]> <S3: glm>
2      b <tibble [329 x 3]> <S3: glm>
3      c <tibble [319 x 3]> <S3: glm>

我的目的是解除模型的嵌套以计算边际效应。我试过了，但出现了这个错误：

> unnest(by_region, mod_rat)
Error: Each column must either be a list of vectors or a list of data frames [mod_rat]

我想知道是否可以在这种类型的 objects (<S3: glm>) 上使用 unnest，如果不能，是否有其他方法来获得这些估计值。

Answer 1

碰巧，margins 包最近有一些更新，可以帮助您以整洁的方式完成此操作。特别是添加了一个 margins_summary() 函数，可以映射到嵌套模型对象。

This issue GitHub 上有详细信息。

这是一些适用于您的示例的代码

使用上面的数据

library(tidyverse)
library(magrittr)
library(margins)

# toy data
test_data <- data.frame(region = sample(letters[1:3], 1000, replace = TRUE),
                             x = sample(0:1, 1000, replace = TRUE), 
                             y = sample(1:100, 1000, replace = TRUE), 
                             z = sample(0:1, 1000, replace = TRUE)) %>% 
arrange(region)

# nest
by_region <- 
    test_data %>%
    group_by(region) %>%
    nest()

# glm function 
region_model  <- function(df) {
   glm(x ~ y + z, data = df, family = "binomial")
}              

# run the model  
by_region %<>% 
  mutate(mod_rat = map(data, region_model))

通过 purrr:map2() 使用 margins_summary() 函数计算边际效应（我已经包含了两种计算逻辑回归边际效应的方法，如包小插图中所述）

by_region %<>% 
    mutate(marginals      = map2(mod_rat, data, ~margins_summary(.x, data = .y)),
           marginals_link = map2(mod_rat, data, ~margins_summary(.x, data = .y, type = "link")))

我们现在可以使用边际效应数据取消嵌套任何一个已创建的列表列

by_region %>% 
   unnest(marginals) -> region_marginals
   region_marginals
# A tibble: 6 x 8
  region factor      AME      SE      z     p
    <fct>  <chr>     <dbl>   <dbl>  <dbl> <dbl>
  1 a      y      -9.38e-4 9.71e-4 -0.966 0.334
  2 a      z       3.59e-2 5.55e-2  0.647 0.517
  3 b      y       1.14e-3 9.19e-4  1.24  0.215
  4 b      z      -2.93e-2 5.38e-2 -0.545 0.586 
  5 c      y       4.67e-4 9.77e-4  0.478 0.633
  6 c      z      -3.32e-2 5.49e-2 -0.604 0.546
# ... with 2 more variables: lower <dbl>,
#   upper <dbl>

而且情节很好

region_marginals %>% 
  ggplot(aes(reorder(factor, AME), AME, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, colour = "#AAAAAA") +
  geom_pointrange() +
  facet_wrap(~region) +
  coord_flip()

未嵌套的 glm 模型

Unnest fitted glm models

nested

r

models