仅使用 tidyverse 手动引导置信区间答案

【问题标题】：Manual bootstrapping for confidence intervals using tidyverse only仅使用 tidyverse 手动引导置信区间
【发布时间】：2022-01-11 10:18:53
【问题描述】：

我有一个分组数据集，我有兴趣汇总一列计数（___ 的数量）。要计算摘要的标准误差，我想在组内引导并计算中位数的标准差。我正在努力弄清楚如何在不使用for 循环的情况下手动对此进行编码（用替换重新采样，而不是像boot() 这样的函数）（即，我希望纯粹是tidyverse解决方案）。如果除了使用*apply() 之外还有其他方法，那将是首选。将整个过程封装到一个函数中会很棒——可以在管道中使用，比如summarise()，或者作为一个可以应用于分组数据的独立函数。

临时数据集可以是 mtcars，我已按 gear 分组。我现在有兴趣使用中位数总结hp 列，并获得相同的置信区间。我已经尝试了一些由稍微相关的线程建议的解决方案，例如replicate()+across()、map()/pmap() 等，但无法让它们适用于我的具体情况。

library(tidyverse)

data <- mtcars %>% 
  select(gear, hp) %>% 
  group_by(gear)

> data
# A tibble: 32 x 2
# Groups:   gear [3]
    gear    hp
   <dbl> <dbl>
 1     4   110
 2     4   110
 3     4    93
 4     3   110
 5     3   175
 6     3   105
 7     3   245
 8     4    62
 9     4    95
10     4   123
# ... with 22 more rows

我希望有一种方法可以将引导结果与简单总结整合为另一列（每组 SE）：

data2 <- data %>% 
  summarise(hp = median(hp))

虽然用齿轮数来概括马力可能没有多大意义，而且hp 的分布可能不是典型的泊松分布，但我认为这个示例的编码解决方案仍然适用于我的具体情况。

编辑 1

解决方案不必是干净且健壮的函数。对于这种特定情况，它可能只是获取每组中的自举 SE 值所需的代码行。所需的输出只是 data2 对象，其中 hp 是中位数列，hpse 是 SE 列。


    data2 <- data %>% 
      summarise(hp = median(hp),
            ### hpse = workingcode()
                )

如果不能在 summarise() 调用中直接以这种方式执行此操作，则必须至少可以稍后将值连接到 data2。

相关话题

使用`boot()`

【问题讨论】：

我很困惑说你不想要 *apply() 功能但不介意咕噜声 map() 家庭。它们大多是等价的——使用 base 的解决方案可以很容易地更新为 tidyverse 样式。此外，解决方案可能会涉及类似于 map(1:B, ...) 的内容，这与 for 循环相同。
@kybazzi 是的，我尝试了map() 家庭，但不知道如何让它适用于我的情况，所以我正在寻找其他方法。我不想要*apply() 函数的主要原因是因为这种方法涉及不同的数据对象类型，但我目前正在考虑将其作为一种选择。我在帖子中不够清楚，但如果它们可以适应 tidyverse 风格，我对这些方法持开放态度，但我当然希望避免 for 循环（如果可能的话）。（将编辑帖子以明确这一点。）
另外一点需要澄清的是解决方案应该是什么样子。就目前而言，这个问题似乎有点过于模糊或开放 - 它可以是一个完整的项目来构建可应用于管道的干净引导功能。你能展示一个所需输出的例子吗？
@kybazzi 我已经编辑了带有说明的帖子。我希望它有所帮助。所需的输出只是使用替换组内的抽样为每个组计算的 SE 列。如果有进一步的疑问，请告诉我。谢谢！
谢谢-您的编辑确实澄清了这一点。我已经给你留下了答案，很高兴回答任何问题。

标签： r statistics tidyverse resampling

【解决方案1】：

首先我们可以创建一个引导函数：

boot_fn = function(x, fn = median, B = 1000) {
  1:B %>%
    # For each iteration, generate a sample of x with replacement
    map(~ x[sample(1:length(x), replace = TRUE)]) %>%
    # Obtain the fn estimate for each bootstrap sample
    map_dbl(fn) %>%
    # Obtain the standard error
    sd()
}

注意我是如何给参数fn 赋予默认值median，这样你就有机会将任何你想要的数字函数传递给boot_fn()。

现在我们可以按照您最初的要求使用该功能：

mtcars %>% 
  group_by(gear) %>%
  summarise(
    hp_median = median(hp), 
    se = boot_fn(hp, fn = median)
  )

# A tibble: 3 x 3
   gear hp_median    se
  <dbl>     <dbl> <dbl>
1     3       180  13.2
2     4        94  15.2
3     5       175  70.3

之所以有效，是因为我们的数据是分组的。对于每个组，x 的新值将发送到 boot_fn()。在这种情况下，传递了三个不同的x 值，每个都是与gear 的每个不同值对应的hp 值。

如果我们只是在函数中添加cat() 语句，这很容易确认：

boot_fn = function(x, fn = median, B = 1000, verbose = FALSE) {
  if (verbose) cat("Hello, x is ", x, "\n")
  1:B %>%
    # For each iteration, generate a sample of x with replacement
    map(~ x[sample(1:length(x), replace = TRUE)]) %>%
    # Obtain the fn estimate for each bootstrap sample
    map_dbl(fn) %>%
    # Obtain the standard error
    sd()
}

data %>%
  summarise(
    hp_median = median(hp), 
    se = boot_fn(hp, fn = median, verbose = TRUE)
  )

输出：

Hello, x is 110 175 105 245 180 180 180 205 215 230 97 150 150 245 175 
Hello, x is 110 110 93 62 95 123 123 66 52 65 66 109 
Hello, x is 91 113 264 175 335 
# A tibble: 3 x 3
   gear hp_median    se
  <dbl>     <dbl> <dbl>
1     3       180  13.5
2     4        94  14.9
3     5       175  69.6

此函数在用于实际数据时可能会出现故障（由于 NAs 之类的原因），但这是一个好的开始。

【讨论】：

感谢您，它适用于当前目的。您在第一步中使用map() 而不是map_dbl() 有什么特别的原因吗？（后者也可以接受向量输入，对吧？）
另外，我知道这超出了我最初的问题的范围，但是否可以将您的解决方案的性能与类似的内容进行简要比较：boot_se <- function(x, fn = median, B = 1000){replicate(B, do.call("fn", list(sample(x, n(), replace = T))), simplify = F) %>% unlist() %>% sd()}？
关于 map() 与 map_dbl()：我建议您将管道更新为 map() 并运行到该行的末尾以查看差异。 map_dbl() 返回一个数字向量 - 每个元素代表从 1 到 B 的引导样本的中位数。使用map() 会改为给出一个列表，但每个元素都是引导样本的中值。在这种情况下，我们想使用map_dbl() 来输入sd()。
关于性能 - 您可以自己尝试这样做，如果出现问题，我很乐意回答任何问题。在您发送的代码中，您可能需要使用sample(x, length(x), replace = TRUE)，因为n() 在 dplyr 管道之外不可用。以下是如何对两个函数进行基准测试的示例：stackoverflow.com/questions/22515525/…
关于map()与map_dbl()：你必须在那里使用map()，因为从1到B的每个复制都需要有自己的x样本。这只能存储在一个列表中，即map() 的输出。另一方面，map_dbl() 应该在每个复制的输出是单个数字时使用。这就是为什么在下一行中使用它的原因 - 从 1 到 B 的每个复制都会计算其原始样本的中值，我们将结果保存在由 map_dbl() 给出的单个向量中。如果还不清楚，请告诉我。

【解决方案2】：

@kybazzi 的解决方案适合流水线工作流程的替代方案是：

boot_se <- function(x, fn = median, B = 100){
  replicate(B,
            do.call("fn", list(sample(x, n(), replace = T))),
            simplify = F) %>% 
    unlist() %>% 
    sd()
}

有时似乎比较慢：


boot_fn = function(x, fn = median, B = 100) {
  1:B %>%
    # For each iteration, generate a sample of x with replacement
    map(~ x[sample(1:length(x), replace = TRUE)]) %>%
    # Obtain the fn estimate for each bootstrap sample
    map_dbl(fn) %>%
    # Obtain the standard error
    sd()
}


data1 <- mtcars %>% 
  select(gear, hp) %>% 
  group_by(gear)

data2 <- data %>% 
  summarise(hpmed = median(hp),
            hpse = boot_se(hp))

data3 <- data %>% 
  summarise(hpmed = median(hp),
            hpse = boot_fn(hp))

#######################################

library(microbenchmark)

microbenchmark((data %>% 
                 summarise(hpmed = median(hp),
                           hpse = boot_fn(hp))),
               (data %>% 
                  summarise(hpmed = median(hp),
                            hpse = boot_se(hp))))

# Output:

Unit: milliseconds
                                                          expr     min       lq
  (data %>% summarise(hpmed = median(hp), hpse = boot_fn(hp))) 14.5737 15.63690
  (data %>% summarise(hpmed = median(hp), hpse = boot_se(hp))) 20.6675 21.64715
     mean   median       uq     max neval
 22.23120 16.78140 25.85675 91.4154   100
 29.15338 22.68525 32.01430 87.6299   100

#######################################

microbenchmark(data2, data3, times = 1000)

# Output:

Unit: nanoseconds
  expr min    lq   mean median  uq  max neval
 data2   0 100.0 95.986    101 101 3501  1000
 data3   0   1.5 92.318    101 101 2700  1000

【讨论】：

编辑 1

相关话题

使用boot()

使用*apply()

使用for 循环

其他

使用`boot()`

使用`*apply()`

使用`for` 循环