使用 dplyr::do 按组拟合许多模型时处理错误的正确方法答案

【问题标题】：Proper way of dealing with errors when fitting many models by group with dplyr::do使用 dplyr::do 按组拟合许多模型时处理错误的正确方法
【发布时间】：2018-09-27 21:30:54
【问题描述】：

使用dplyr::do，可以非常简单地按组拟合多个模型，如下所示：

library(tidyverse)
set.seed(100)
tbl <- tibble(
  group_id = rep(1:3, each = 10),
  y1 = rnorm(30),
  y2 = runif(30),
  x1 = rnorm(30),
  x2 = runif(30)
)

tbl %>%
  group_by(group_id) %>%
  do(
    model1 = lm(y1 ~ x1 + x2, data = .),
    model2 = lm(y2 ~ x1 + x2, data = .)
  )
#> Source: local data frame [3 x 3]
#> Groups: <by row>
#> 
#> # A tibble: 3 x 3
#>   group_id model1   model2  
#> *    <int> <list>   <list>  
#> 1        1 <S3: lm> <S3: lm>
#> 2        2 <S3: lm> <S3: lm>
#> 3        3 <S3: lm> <S3: lm>

这是broom::tidy 和broom::glance 用于按组提取r.squared 和系数的理想格式。但是，当一组（此处为 group_id == 3）具有所有缺失值时，就会出现问题：

tbl2 <- mutate(tbl, y2 = c(runif(20), rep(NA, 10)))

tbl2 %>%
  group_by(group_id) %>%
  do(
    model1 = lm(y1 ~ x1 + x2, data = .),
    model2 = lm(y2 ~ x1 + x2, data = .)
  )
#> Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...): 0 (non-NA) cases

正如预期的那样，因为group_id == 3 没有y2 的非缺失值，所以model2 不适合任何东西。我发现的其他问题建议在拟合之前简单地删除具有NA 值的行，但是我不想这样做，因为那样我会失去model1 的成功拟合。我想到的另一种方法是使用try 捕获错误，但我无法仅用缺失值替换错误。我尝试了以下使用purrr::modify_if 的代码的许多变体，但不知道为什么没有替换该值（例如，

modify_if(list(1, "a", TRUE), ~ inherits(., "numeric"), `is.na<-`)

工作正常。）您可以看到，使用map 和inherits 正确地发现了哪些单元格是try-error 类，但是将其包裹在modify_if 中使其不再被发现。

tbl2 %>%
  group_by(group_id) %>%
  do(
    model1 = lm(y1 ~ x1 + x2, data = .),
    model2 = try(
      lm(y2 ~ x1 + x2, data = .),
      silent = TRUE
    )
  ) %>%
  ungroup() %>%
  mutate_all(
    function(col) map_lgl(col, function(cell) inherits(cell, "try-error"))
  )
#> # A tibble: 3 x 3
#>   group_id model1 model2
#>   <lgl>    <lgl>  <lgl> 
#> 1 FALSE    FALSE  FALSE 
#> 2 FALSE    FALSE  FALSE 
#> 3 FALSE    FALSE  TRUE

tbl2 %>%
  group_by(group_id) %>%
  do(
    model1 = lm(y1 ~ x1 + x2, data = .),
    model2 = try(
      lm(y2 ~ x1 + x2, data = .),
      silent = TRUE
    )
  ) %>%
  ungroup() %>%
  mutate_at(
    .vars = vars(starts_with("model_")),
    .funs = function(col) {
      modify_if(
        .x = col,
        .p = function(cell) inherits(cell, "try-error"),
        .f = function(cell) unclass(`is.na<-`(cell)))
    }
  )
#> # A tibble: 3 x 3
#>   group_id model1   model2         
#> *    <int> <list>   <list>         
#> 1        1 <S3: lm> <S3: lm>       
#> 2        2 <S3: lm> <S3: lm>       
#> 3        3 <S3: lm> <S3: try-error>

由reprex package (v0.2.0) 于 2018 年 4 月 17 日创建。

我的实际数据有约 80k 组和约 10 个模型供参考。任何有关改进此代码或捕获错误的更好方法的建议将不胜感激。

【问题讨论】：

标签： r dplyr lm

【解决方案1】：

我认为这是我发现的解决此问题的最佳方法。与其使用modify 来尝试替换错误模型，不如将它们过滤掉并替换glance 之后的缺失行。这是因为glance 无论如何都不能很好地处理格式错误的lm 输出。

tbl2 %>%
  group_by(group_id) %>%
  do(
    model1 = lm(y1 ~ x1 + x2, data = .),
    model2 = try(
      lm(y2 ~ x1 + x2, data = .),
      silent = TRUE
    )
  ) %>%
  ungroup() %>%
  gather(model, lm, starts_with("model")) %>%
  mutate(error = map_lgl(lm, ~inherits(., "try-error"))) %>%
  filter(error == FALSE) %>%
  rowwise() %>%
  glance(lm) %>%
  ungroup() %>%
  complete(group_id = 1:3, model = c("model1", "model2"))
#> # A tibble: 6 x 14
#>   group_id model  error r.squared adj.r.squared  sigma statistic p.value
#>      <int> <chr>  <lgl>     <dbl>         <dbl>  <dbl>     <dbl>   <dbl>
#> 1        1 model1 FALSE    0.0215       -0.258   0.629    0.0769   0.927
#> 2        1 model2 FALSE    0.107        -0.149   0.329    0.418    0.674
#> 3        2 model1 FALSE    0.208        -0.0184  0.868    0.919    0.442
#> 4        2 model2 FALSE    0.0808       -0.182   0.362    0.308    0.745
#> 5        3 model1 FALSE    0.0707       -0.195   0.738    0.266    0.774
#> 6        3 model2 NA      NA            NA      NA       NA       NA    
#> # ... with 6 more variables: df <int>, logLik <dbl>, AIC <dbl>, BIC <dbl>,
#> #   deviance <dbl>, df.residual <int>

【讨论】：