机器算法验证 - 使用 H0 下的 bootstrap 对两种均值的差异进行检验：组内替换或合并样本内的替换 - 吾爱随笔录

使用 H0 下的 bootstrap 对两种均值的差异进行检验：组内替换或合并样本内的替换

机器算法验证 r 假设检验引导程序小样本置换检验

2022-01-19 14:45:03

假设我有一个包含两个独立组的数据：

g1.lengths <- c (112.64, 97.10, 84.18, 106.96, 98.42, 101.66)

g2.lengths <- c (84.44, 82.10, 83.26, 81.02, 81.86, 86.80, 
                     85.84, 97.08, 79.64, 83.32, 91.04, 85.92,
                     73.52, 85.58, 97.70, 89.72, 88.92, 103.72,
                     105.02, 99.48, 89.50, 81.74)

group = rep (c ("g1", "g2"), c (length (g1.lengths), length (g2.lengths)))

lengths = data.frame( lengths = c(g1.lengths, g2.lengths), group)

很明显，每组的样本量是有偏差的，其中g1 有 6 个观测值，而g2 有 22个观测值。传统的 ANOVA 表明，当临界值设置为 0.05（p 值为0.0044）时，组具有不同的均值。

summary (aov (lengths~group, data = lengths))

鉴于我的目标是比较均值差异，这种不平衡和小样本数据可能会与传统方法给出不合适的结果。因此，我想执行置换测试和引导。

排列测试

零假设 (H0) 表明组的均值相同。通过将组合并到一个样本中来证明置换检验中的这一假设是合理的。这确保了两组的样本来自相同的分布。通过从合并数据中重复采样（或更准确地说 - 重新洗牌），观察以新的方式重新分配（洗牌）到样本，并计算测试统计量。执行 n 次，将在 H0 为 TRUE 的假设下给出测试统计的抽样分布。最后，在 H0 下，p 值是检验统计量等于或超过观察值的概率。

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.p <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.p <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool))

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.p[i] = mean (g1.perm) - mean (g2.perm) 
}
p.permute <- (sum (abs (sampl.dist.p) >= abs(obs.diff.p)) + 1)/ (iterations+1)

置换检验的报告 p 值为0.0053。好的，如果我做得正确，排列和参数方差分析会给出几乎相同的结果。

引导程序

首先，我知道当样本量太小时，bootstrap 无济于事。这篇文章表明它可能会更糟糕和误导。此外，第二个强调当假设检验是主要目标时，置换检验通常比引导法更好。尽管如此，这篇伟大的文章解决了计算机密集型方法之间的重要差异。然而，在这里我想提出（我相信）一个不同的问题。

让我首先介绍最常见的引导方法（Bootstrap1：在池样本中重新采样）：

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b1 <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.b1 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool), replace = TRUE) 
        # "replace = TRUE" is the only difference between bootstrap and permutations

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.b1[i] = mean (g1.perm) - mean (g2.perm) 
}
p.boot1 <- (sum (abs (sampl.dist.b1) >= obs.diff.b1) + 1)/ (iterations+1)

以这种方式执行的引导的 P 值为0.005。即使这听起来很合理并且与参数方差分析和置换检验几乎相同，基于我们刚刚合并从中抽取后续样本的样本来证明这个引导程序中的 H0 是否合适？

我在几篇科学论文中发现了不同的方法。具体来说，我看到研究人员修改数据以在引导之前满足 H0。环顾四周，我在 CV中发现了非常有趣的帖子， @jan.s在帖子问题中解释了引导程序的异常结果，其目的是比较两种方法。但是，在那篇文章中没有介绍在引导之前修改数据时如何执行引导。在引导之前修改数据的方法如下所示：

H0 表示两组的均值相同
如果我们从合并样本的平均值中减去单个观察值，H0 成立

在这种情况下，数据的修改应影响组均值，从而影响其差异，但不会影响组内（和组间）的变化。

修改后的数据将成为进一步引导的基础，但需要注意的是，抽样是在每个组内分别进行的。
计算 g1 和 g2 的自举平均值之间的差异，并与组间观察到的（未修改的）差异进行比较。
与观察到的值相等或更多的极值除以迭代次数的比例将给出 p 值。

这是代码（Bootstrap2：在修改 H0 为 TRUE 后在组内重新采样）：

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b2 <- mean (g1.lengths) - mean (g2.lengths)

# make H0 to be true (no difference between means of two groups)
H0 <- pool - mean (pool)

# g1 from H0 
g1.H0 <- H0[1:s.size.g1] 

# g2 from H0
g2.H0 <- H0[(s.size.g1+1):length(pool)]

iterations <- 10000
sampl.dist.b2 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        # Sample with replacement in g1
        g1.boot = sample (g1.H0, replace = T)

        # Sample with replacement in g2
        g2.boot = sample (g2.H0, replace = T)

        # bootstrapped difference
        sampl.dist.b2[i] <- mean (g1.boot) - mean (g2.boot)  
}
p.boot2 <- (sum (abs (sampl.dist.b2) >= obs.diff.b2) + 1)/ (iterations+1)

这种执行的引导程序将给出0.514的 p 值，这与之前的测试相比有很大的不同。我相信这必须处理@jan.s'的解释，但我无法弄清楚关键在哪里......

1个回答

这是我对它的看法，基于 Efron 和 Tibshirani 的An Introduction to the bootstrap的第 16 章（第 220-224 页）。简而言之，您的第二个引导算法实施错误，但总体思路是正确的。

在进行引导测试时，必须确保重采样方法生成与原假设相对应的数据。我将使用 R 中的睡眠数据来说明这篇文章。请注意，我使用的是学生化测试统计数据，而不仅仅是教科书推荐的均值差异。

经典 t 检验使用分析结果来获取有关 t 统计量的抽样分布的信息，产生以下结果：

x <- sleep$extra[sleep$group==1]
y <- sleep$extra[sleep$group==2]
t.test(x,y)
t = -1.8608, df = 17.776, p-value = 0.07939

一种方法在精神上类似于更知名的置换测试：样本是在整个观察集合中获取的，同时忽略分组标签。然后第一个 $n1$ 被分配到第一组，其余的 $n2$ 到第二组。

# pooled sample, assumes equal variance
pooled <- c(x,y)
for (i in 1:10000){
  sample.index <- sample(c(1:length(pooled)),replace=TRUE)
  sample.x <- pooled[sample.index][1:length(x)]
  sample.y <- pooled[sample.index][-c(1:length(y))]
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.pooled <-  (1 + sum(abs(boot.t) >= abs(t.test(x,y)$statistic))) / (10000+1) 
p.pooled
[1] 0.07929207

但是，这个算法实际上是在测试 x 和 y 的分布是否相同。如果我们只是对它们的总体均值是否相等感兴趣，而不对它们的方差做任何假设，我们应该在下面生成数据 $H_0$ 以稍微不同的方式。你的方法是正确的，但你的翻译是 $H_0$ 与教科书中提出的有点不同。生成 $H_0$ 我们需要从第一组的观察中减去第一组的平均值，然后加上共同或合并的平均值 $\bar{z}$ . 对于第二组，我们做同样的事情。

{\tilde{x}}_{i} = x_{i} - \bar{x} + \bar{z}

$\tilde{x}_i = x_i - \bar{x} + \bar{z}$

{\tilde{y}}_{i} = y_{i} - \bar{y} + \bar{z}

$\tilde{y}_i = y_i - \bar{y} + \bar{z}$

当您计算新变量的均值时，这变得更加直观 $\tilde{x}/\tilde{y}$ . 通过首先减去它们各自的组均值，变量以零为中心。通过添加整体平均值 $\bar{z}$ 我们最终得到一个以整体平均值为中心的观察样本。换句话说，我们对观测值进行了变换，使它们具有相同的均值，这也是两组加在一起的总体均值，也就是 $H_0$ .

# sample from H0 separately, no assumption about equal variance
xt <- x - mean(x) + mean(sleep$extra)
yt <- y - mean(y) + mean(sleep$extra)

boot.t <- c(1:10000)
for (i in 1:10000){
  sample.x <- sample(xt,replace=TRUE)
  sample.y <- sample(yt,replace=TRUE)
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.h0 <-  (1 + sum(abs(boot.t) >= abs(t.test(x,y)$statistic))) / (10000+1) 
p.h0
[1] 0.08049195

这一次，我们最终得到了三种方法的相似 p 值。

其它你可能感兴趣的问题

上一篇新闻中的方程式：将多层次模型翻译给普通观众下一篇后验与先验和可能性非常不同