带有两个子组的嵌套 glm答案

【问题标题】：r nested glm with two subgroups带有两个子组的嵌套 glm
【发布时间】：2022-11-02 02:57:58
【问题描述】：

我正在尝试像这样运行一个简单的 glm 模型。

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

         data("mtcars")
         head(mtcars)
         
         mtcars$Name <- row.names(mtcars)
         row.names(mtcars) <- NULL

         glm(mpg ~ wt, data=mtcars)

到目前为止没有任何问题。

接下来我试图在gear 的每个子组上运行这个模型，即gear=3, gear=4, gear=5，所以我在这样的 dlply 函数中运行我的 glm 模型。


Model1 <- plyr::dlply(mtcars, "gear",
            function(x)
              tryCatch(
                glm(mpg ~ wt,
                    data =x ),
                error = function(e) NA), .drop = TRUE)

SummaryCars <- map2_df(Model1,
                        names(Model1),
                          ~broom::tidy(.x, confint = TRUE)[2,] %>%
                              mutate(gear = .y))

现在我有了第三个子组carb。这个变量有 6 个级别

table(mtcars$carb)

1  2  3  4  6  8 
7 10  3 10  1  1

排除碳水化合物水平 6 和 8。我喜欢在 carb 水平 1、2、3 和 4 上运行我的模型。对于gear 的每个级别。但我喜欢在每次迭代中排除一级碳水化合物。

型号1- 碳水化合物水平 1,2,3 (从 carb= 4 中排除数据)

 ```
 Model1 <- plyr::dlply(mtcars, "gear",
        function(x)
          tryCatch(
            glm(mpg ~ wt,
                data =x ),
            error = function(e) NA), .drop = TRUE)

  SummaryCars <- map2_df(Model1,
                    names(Model1),
                      ~broom::tidy(.x, confint = TRUE)[2,] %>%
                          mutate(gear = .y))
   ```

**Model2 ** - 碳水化合物等级 1,2,4 (从 carb= 3 中排除数据)

 ```
 Model1 <- plyr::dlply(mtcars, "gear",
        function(x)
          tryCatch(
            glm(mpg ~ wt,
                data =x ),
            error = function(e) NA), .drop = TRUE)

  SummaryCars <- map2_df(Model1,
                    names(Model1),
                      ~broom::tidy(.x, confint = TRUE)[2,] %>%
                          mutate(gear = .y))
   ```

**Model3 ** - 碳水化合物等级 1,3,4 (从 carb= 2 中排除数据)

 ```
 Model1 <- plyr::dlply(mtcars, "gear",
        function(x)
          tryCatch(
            glm(mpg ~ wt,
                data =x ),
            error = function(e) NA), .drop = TRUE)

  SummaryCars <- map2_df(Model1,
                    names(Model1),
                      ~broom::tidy(.x, confint = TRUE)[2,] %>%
                          mutate(gear = .y))
   ```

**Model4l ** - 碳水化合物等级 2,3,4 (从 carb= 1 中排除数据)

 ```
 Model1 <- plyr::dlply(mtcars, "gear",
        function(x)
          tryCatch(
            glm(mpg ~ wt,
                data =x ),
            error = function(e) NA), .drop = TRUE)

  SummaryCars <- map2_df(Model1,
                    names(Model1),
                      ~broom::tidy(.x, confint = TRUE)[2,] %>%
                          mutate(gear = .y))
   ```

如您所见，我可以在 Gear (3,4,5) 的每个级别中运行模型，但我不确定如何在此之上添加另一个循环，其中一个级别的数据被排除并考虑休息。

预期的最终结果

         Model            SubModel     Estimate    Lower(CI)      Upper(CI)    stdError   p
         Exclude carb= 4  Gear = 3     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 4  Gear = 4     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 4  Gear = 5     xxx         xxxx           xxxx         xxx        x

         Exclude carb= 3  Gear = 3     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 3  Gear = 4     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 3  Gear = 5     xxx         xxxx           xxxx         xxx        x

         Exclude carb= 2  Gear = 3     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 2  Gear = 4     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 2  Gear = 5     xxx         xxxx           xxxx         xxx        x

         Exclude carb= 1  Gear = 3     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 1  Gear = 4     xxx         xxxx           xxxx         xxx        x
         Exclude carb= 1  Gear = 5     xxx         xxxx           xxxx         xxx        x

任何帮助深表感谢。提前致谢。

我在我的问题中包含了代码，显示了我尝试过的内容。

【问题讨论】：

标签： r dplyr apply glm

【解决方案1】：

你可以这样做：

mtcars %>%
  filter(!carb %in% c(6,8))%>%
  mutate(Exclude = list(sort(unique(carb))))%>%
  group_by(gear)%>%
  summarise(.groups = 'drop',
            Exclude = Exclude[[1]],
            result = map(Exclude, 
                  ~tidy(glm(mpg~wt, data = cur_data(), subset = carb!=.x))))%>%
  unnest(result)
# A tibble: 24 × 7
    gear Exclude term        estimate std.e…¹ stati…² p.value
   <dbl>   <dbl> <chr>          <dbl>   <dbl>   <dbl>   <dbl>
 1     3       1 (Intercept)    25.1    3.51     7.14 3.15e-5
 2     3       1 wt             -2.44   0.842   -2.90 1.59e-2
 3     3       2 (Intercept)    28.9    3.01     9.58 5.11e-6
 4     3       2 wt             -3.28   0.733   -4.47 1.55e-3
 5     3       3 (Intercept)    28.4    3.15     9.04 3.96e-6
 6     3       3 wt             -3.18   0.786   -4.04 2.36e-3
 7     3       4 (Intercept)    29.8    5.18     5.76 4.25e-4
 8     3       4 wt             -3.43   1.47    -2.33 4.82e-2
 9     4       1 (Intercept)    37.8    4.39     8.61 1.35e-4
10     4       1 wt             -5.37   1.49    -3.60 1.14e-2
# … with 14 more rows, and abbreviated variable names
#   ¹std.error, ²statistic
# ℹ Use `print(n = ...)` to see more rows

【讨论】：

就像我对 zephryl 提到的一样，您的代码在这个玩具数据集上完美运行，但是如果有任何子组没有数据 (NA)，那么这个代码就会失败。我之前在我的代码中使用了trycatch 函数来处理这些情况，但是如果某些子组没有数据，我不确定如何包含trycatch 类型的错误处理。
@EagleHawk 在运行整洁部分之前使用 if else 。即if(nrow(na.omit(filter(cur_data(), carb!=.x)))>2) glm(mpg~wt, data = cur_data(), subset = carb!=.x) else NA
哇，你真快。谢谢，我会试试的。

【解决方案2】：

该解决方案使用嵌套的purrr::map_dfr() 调用围绕broom::tidy(lm())。

编辑：OP 要求提供一种处理没有有效数据的子集的解决方案。我现在通过将tidy(lm()) 调用包装在purrr::possibly() 中来做到这一点。输出将包括无效模型的单行 NA 模型参数。

library(dplyr)
library(purrr)
library(broom)

# rmv carb == 6 or 8, add test case with all missing data, then split by gear
mtc_gears <- filter(mtcars, carb <= 4) %>% 
  mutate(wt = if_else(gear == 4 & carb != 4, NA_real_, wt))
mtc_gears <- split(mtc_gears, mtc_gears$gear)

safe_tidy_lm <- possibly(
  function(formula, data) tidy(lm(formula, data = data), conf.int = TRUE),
  tibble(.rows = 1)
)

map_dfr(
  1:4, 
  (carb_drop) map_dfr(
      mtc_gears, 
      (mtc_gear) safe_tidy_lm(
        mpg ~ wt, 
        data = filter(mtc_gear, carb != carb_drop)
      ), 
      .id = "Gear"
    ), 
  .id = "Exclude Carb"
)

#> Warning in qt(a, object$df.residual): NaNs produced
#> Warning in qt(a, object$df.residual): NaNs produced
#> # A tibble: 23 × 9
#>    `Exclude Carb` Gear  term     estim…¹ std.e…² stati…³ p.value conf.…⁴ conf.…⁵
#>    <chr>          <chr> <chr>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
#>  1 1              3     (Interc…   25.1    3.51     7.14 3.15e-5   17.2   32.9  
#>  2 1              3     wt         -2.44   0.842   -2.90 1.59e-2   -4.32  -0.563
#>  3 1              4     (Interc…   30.2    3.61     8.37 1.40e-2   14.7   45.7  
#>  4 1              4     wt         -3.38   1.16    -2.92 1.00e-1   -8.37   1.61 
#>  5 1              5     (Interc…   44.4    1.82    24.3  2.62e-2   21.2   67.5  
#>  6 1              5     wt         -8.92   0.769  -11.6  5.47e-2  -18.7    0.846
#>  7 2              3     (Interc…   28.9    3.01     9.58 5.11e-6   22.1   35.7  
#>  8 2              3     wt         -3.28   0.733   -4.47 1.55e-3   -4.93  -1.62 
#>  9 2              4     (Interc…   30.2    3.61     8.37 1.40e-2   14.7   45.7  
#> 10 2              4     wt         -3.38   1.16    -2.92 1.00e-1   -8.37   1.61 
#> # … with 13 more rows, and abbreviated variable names ¹estimate, ²std.error,
#> #   ³statistic, ⁴conf.low, ⁵conf.high

^{创建于 2022-11-01，reprex v2.0.2}

【讨论】：

这非常适用于玩具数据集。但是，如果某些级别具有 NA 或某些子组没有数据，则代码不起作用。我在 trycatch 函数中调用我的 glm 以排除没有任何数据的子组 tryCatch( glm(mpg ~ wt, data =x ), error = function(e) NA), .drop = TRUE) 但我不确定如何在您的解决方案中实现类似的错误处理机制。
请参阅我编辑的答案。
我会测试一下，谢谢。