重复测量引导统计，按多个因素分组

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))

所以我们有 2 个因子（每个因子有 2 个水平，因此有 4 个组合）和一个连续测量。我们还有一个重复测量设计，每个单元格中有多个measure，对应于相同的id。

我尝试先解决 groupby 问题，然后解决引导问题，然后将两者结合起来，但几乎陷入困境......

统计数据，按 2 个因素分组

我可以通过以下方式获得 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% 的引导 CI，对于每个统计数据，对于 4 个单元格中的每一个。 我不介意我是否必须运行一次最终解决方案统计数据（例如平均值、中位数等），但一次性完成所有操作的奖励积分。

引导重复测量

不能完全完成这项工作，但我想要的是 95% 的引导 CI，以适合这种重复测量设计的方式完成。除非我弄错了，否则我想根据id（不是基于数据框的行）选择引导样本，然后计算一个汇总度量（例如mean） 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)

Q1：我在通过id 选择引导样本时遇到错误，即d <- data[ data$id==indicies, ] 行

结合 bootstrap 和 groupby 2 个因素

Q2：我不知道如何将这两种方法结合在一起以达到最终的预期结果。我唯一的想法是将aggregate 调用放在myfunc 中，以重复计算每个引导复制下的单元格统计信息，但我在这里使用R 超出了我的舒适区。

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 数据框中，每个组都有一行，以及您的函数产生的值的列。所以我们简单地再次对数据框进行分组，然后调用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)

希望对你有帮助！

Answer 2

对于带有集群变量的引导程序，这是一个无需额外软件包的解决方案。不过我没有使用boot 包。

第 1 部分：引导

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

.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) )

}

第 2 部分：Boostrap 估计和 CI

下一个函数抽取多个样本并将相同的aggregate 语句应用于每个样本。然后通过mean 和quantile 获得引导估计和CI。

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 接口进行聚合，我也更喜欢它。

第 3 部分：示例

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

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