重复测量自助统计，按多个因素分组

Question

我有一个数据框，看起来像这样，但显然有更多的行等：

df <- data.frame(id=c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2),
                 cond=c('A', 'A', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'B', 'B'),
                 comm=c('X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y','X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y'),
                 measure=c(0.8, 1.1, 0.7, 1.2, 0.9, 2.3, 0.6, 1.1, 0.7, 1.3, 0.6, 1.5, 1.0, 2.1, 0.7, 1.2))

所以我们有两个因素（每个因素都有两个水平，因此有4种组合）和一个连续的测量指标。我们还有重复测量设计，因为我们在每个单元格内有多个与相同id对应的测量。我已经尝试先解决groupby问题，然后是bootstrap问题，最后将两者结合起来，但基本上卡住了...

按照两个因素分组的统计数据

我可以通过以下方式获取每个4个单元格的多个摘要统计数据：

summary_stats <- aggregate(df$measure, 
                           by = list(df$cond, df$comm),
                           function(x) c(mean = mean(x), median = median(x), sd = sd(x)))
print(summary_stats)

导致结果为：

  Group.1 Group.2     x.mean   x.median       x.sd
1       A       X 0.85000000 0.85000000 0.12909944
2       B       X 0.65000000 0.65000000 0.05773503
3       A       Y 1.70000000 1.70000000 0.58878406
4       B       Y 1.25000000 1.20000000 0.17320508

这很棒，因为我们可以获得每个4个单元格的多个统计数据。

但是我真正想要的是95%的bootstrap置信区间，对于每个统计数据和每个4个单元格。如果必须运行最终解决方案以获取统计数据（例如均值、中位数等），那也没关系，但一次性完成会更好。

重复测量的Bootstrap 我无法完全理解这个问题，但我想要的是95%的bootstrap置信区间，以适应此重复测量设计。除非我弄错了，否则我希望根据id选择bootstrap样本（而不是基于数据框的行），然后计算每个4个单元格的摘要统计量（例如平均值）。

library(boot)
myfunc <- function(data, indices) {
   # select bootstrap sample to index into `id`
   d <- data[data$id==indicies,]
   return(c(mean=mean(d), median=median(d), sd = sd(d)))
}

bresults <- boot(data = CO2$uptake, statistic = myfunc, R = 1000)

问题1：我在选择Bootstrap样本时遇到了错误，即行d <- data[ data$id==indicies, ]

结合Bootstrap和groupby 2个因素

问题2：我不知道如何将这两种方法结合起来实现最终的期望结果。我的唯一想法是将aggregate调用放在myfunc中，以便在每个Bootstrap复制下重复计算单元格统计信息，但是对于R来说，我已经超出了自己的舒适区。

- Ben Vincent

Answer 1

通过你的两个问题，你有两个问题：

1. 如何引导（重新采样）你的数据，以便你可以基于 id 而不是行进行重新采样
2. 如何对 2x2 设计中的四个组执行单独的引导

一个简单的方法是使用以下软件包（所有这些都是 tidyverse 的一部分）：

• dplyr 用于操作数据（特别是汇总每个 id 的数据），还可用于整洁的 %>% 前向管道运算符，它将表达式的结果作为下一个表达式的第一个参数，因此您可以链接命令
• broom 用于在数据框中为每个组执行操作
• boot（您已经使用）用于引导

---

加载软件包：

library(dplyr)
library(broom)
library(boot)

首先，为了确保在重新采样时是否包括一个主题，我会将每个主题的各种值保存为列表：

df <- df %>%
    group_by(id, cond, comm) %>%
    summarise(measure=list(measure)) %>%
    ungroup()

现在数据框的行数减少了（每个ID有4行），而变量measure不再是数字（而是一个列表）。这意味着我们可以只使用boot提供的索引（解决问题1），但实际进行计算时，我们也必须对其进行"unlist"操作，所以您的函数现在如下：

myfunc <- function(data, indices) {
    data <- data[indices,]
    return(c(mean=mean(unlist(data$measure)),
             median=median(unlist(data$measure)),
             sd = sd(unlist(data$measure))))
}

现在我们可以使用 boot 来对每一行进行重新抽样，我们可以考虑如何整洁地按组进行操作。这就是broom包的用处：您可以要求它对数据框中的每个组执行一个操作，并将其存储在一个整洁的数据帧中，每个分组有一行，函数生成值的列。因此，我们再次对数据框进行分组，然后调用do（tidy(...)），用 . 代替变量名。这样，您就可以解决问题2了！

bootresults <- df %>%
    group_by(cond, comm) %>%
    do(tidy(boot(data = ., statistic = myfunc, R = 1000)))

这会产生以下结果：

# Groups:   cond, comm [4]
     cond   comm   term  statistic         bias    std.error
   <fctr> <fctr>  <chr>      <dbl>        <dbl>        <dbl>
 1      A      X   mean 0.85000000  0.000000000 5.280581e-17
 2      A      X median 0.85000000  0.000000000 5.652979e-17
 3      A      X     sd 0.12909944 -0.004704999 4.042676e-02
 4      A      Y   mean 1.70000000  0.000000000 1.067735e-16
 5      A      Y median 1.70000000  0.000000000 1.072347e-16
 6      A      Y     sd 0.58878406 -0.005074338 7.888294e-02
 7      B      X   mean 0.65000000  0.000000000 0.000000e+00
 8      B      X median 0.65000000  0.000000000 0.000000e+00
 9      B      X     sd 0.05773503  0.000000000 0.000000e+00
10      B      Y   mean 1.25000000  0.001000000 7.283065e-02
11      B      Y median 1.20000000  0.027500000 7.729634e-02
12      B      Y     sd 0.17320508 -0.030022214 5.067446e-02

希望这是您想要看到的内容！

---

如果您想进一步使用数据框中的值，可以使用其他dplyr函数来选择查看此表中的哪些行。例如，要查看条件A / X的测量标准差的自助法标准误差，可以执行以下操作：

bootresults %>% filter(cond=='A', comm=='X', term=='sd') %>% pull(std.error)

希望这有所帮助！

- janfreyberg

---

非常好，但我认为Bootstrap对于我特定的重复测量情况不太合适。目前它是按每个id``cond``comm组合进行选择。这比在原始df中按行选择要好，但与仅按id选择略有不同。您的答案是否可以修改，以便在列表中包括cond和comm以及measure，使得修改后的df中的每一行现在都对应一个单独的id？ - Ben Vincent

我认为它已经是这样了 - 目前它的作用是对每个id重复的测量进行总结，但保留cond和comm的单独行。然后，它按cond和comm对数据进行分组，有效地为两个因素的每种组合创建一个新的数据框架。然后，它针对这四个数据框架中的每一个进行自助法，并根据行进行重新采样（一旦它们被分组，每个子组只有2行，一个用于每个id）。 - janfreyberg

除非您希望确保在每次引导迭代中，每个条件下都重新采样相同的主题组？ - janfreyberg

是的，这就是它的作用。最初调试dplyr管道可能有点棘手，但检查分组是否有效的一个好方法是检查每个组有多少行：df %>% group_by(cond, comm) %>% count()。这使我得到每个组2行，这些行将是boot在do(tidy(...))调用内部重新采样的行。实际上，当您通过group_by(...) %>% do(tidy(...))进行操作时，do(tidy())调用内部的函数只会看到数据框的每个组，而不是整个数据框。 - janfreyberg

@BenVincent 只是提醒一下，我最终在 tidy 中发现了一个 bug，并且已经 开了一个 issue 来尝试修复它。使用这种方法，现在有点难以得到正确的置信区间，但肯定是可能的。稍后我会再看一下！ - janfreyberg

Answer 2

针对带有集群变量的引导，以下是一种不需要额外软件包的解决方案。我没有使用boot软件包。

第一部分：引导

该函数从一组聚类观察中随机抽取样本。

.clusterSample <- function(x, id){

  boot.id <- sample(unique(id), replace=T)
  out <- lapply(boot.id, function(i) x[id%in%i,])

  return( do.call("rbind",out) )

}

第二部分：Bootstrap估计和置信区间 下一步的函数将绘制多个样本，并对每个样本应用相同的aggregate语句。通过使用mean和quantile方法，可以得到Bootstrap估计和置信区间。

clusterBoot <- function(data, formula, cluster, R=1000, alpha=.05, FUN){

  # cluster variable
  cls <- model.matrix(cluster,data)[,2]

  template <- aggregate(formula, .clusterSample(data,cls), FUN)
  var <- which( names(template)==all.vars(formula)[1] )
  grp <- template[,-var,drop=F]
  val <- template[,var]

  x <- vapply( 1:R, FUN=function(r) aggregate(formula, .clusterSample(data,cls), FUN)[,var],
               FUN.VALUE=val )

  if(is.vector(x)) dim(x) <- c(1,1,length(x))
  if(is.matrix(x)) dim(x) <- c(nrow(x),1,ncol(x))

  # bootstrap estimates
  est <- apply( x, 1:2, mean )
  lo <- apply( x, 1:2, function(i) quantile(i,alpha/2) )
  up <- apply( x, 1:2, function(i) quantile(i,1-alpha/2) )
  colnames(lo) <- paste0(colnames(lo), ".lo")
  colnames(up) <- paste0(colnames(up), ".up")

  return( cbind(grp,est,lo,up) )

}

请注意使用vapply。我使用它是因为我更喜欢使用数组而不是列表。还要注意，我使用了formula接口进行聚合，这也更好用。

第三部分：示例

基本上可以与任何类型的统计数据一起使用，甚至没有分组变量。一些例子包括：

myStats <- function(x) c(mean = mean(x), median = median(x), sd = sd(x))

clusterBoot(data=df, formula=measure~cond+comm, cluster=~id, R=10, FUN=myStats)
#   cond comm mean median         sd mean.lo median.lo      sd.lo mean.up median.up      sd.up
# 1    A    X 0.85  0.850 0.11651125    0.85      0.85 0.05773503    0.85      0.85 0.17320508
# 2    B    X 0.65  0.650 0.05773503    0.65      0.65 0.05773503    0.65      0.65 0.05773503
# 3    A    Y 1.70  1.700 0.59461417    1.70      1.70 0.46188022    1.70      1.70 0.69282032
# 4    B    Y 1.24  1.215 0.13856406    1.15      1.15 0.05773503    1.35      1.35 0.17320508

clusterBoot(data=df, formula=measure~cond+comm, cluster=~id, R=10, FUN=mean)
#   cond comm  est  .lo  .up
# 1    A    X 0.85 0.85 0.85
# 2    B    X 0.65 0.65 0.65
# 3    A    Y 1.70 1.70 1.70
# 4    B    Y 1.25 1.15 1.35

clusterBoot(data=df, formula=measure~1, cluster=~id, R=10, FUN=mean)
#      est    .lo    .up
# 1 1.1125 1.0875 1.1375

- SimonG

---

干得好。为了完整性，需要添加 library(boot) 来访问 boot 函数。 - Ben Vincent

谢谢你发现了这个问题，已经添加了。 - SimonG

这很棒。但我正在尝试弄清楚引导是否基于原始数据框中的id进行抽样。我的初学者眼睛怀疑不是，因为aggregate仅将df$measure作为数据传递给bootGroup函数，这意味着在.bootStats中选择是基于行而不是实际的id？ - Ben Vincent

你的猜测是正确的：它确实在组内进行引导。具体来说，aggregate按组拆分数据，然后调用bootGroup，该函数返回每个分组变量组合的引导估计和置信区间。这是否是期望的行为（或者不是）？ - SimonG

接近了，但我们真正需要根据个体进行引导选择，标记为“id”。这是否需要将aggregate函数更改为bootGroup函数？ - Ben Vincent

那很有道理。看看我的更新答案，检查它是否符合你的要求。这个新答案只使用基础 R，因为我不知道如何使用boot包进行聚类抽样。此外，顶层函数现在使用公式来控制aggregate的操作。 - SimonG