使用 H0 下的 bootstrap 對兩種均值的差異進行檢驗：組內替換或合併樣本內的替換
假設我有一個包含兩個獨立組的數據：
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’的解釋，但我無法弄清楚關鍵在哪裡……
這是我對它的看法，基於 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一種方法在精神上類似於更知名的置換測試：樣本是在整個觀察集合中獲取的，同時忽略分組標籤。然後第一個被分配到第一組，其餘的到第二組。
# 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 的分佈是否相同。如果我們只是對它們的總體均值是否相等感興趣，而不對它們的方差做任何假設，我們應該在下面生成數據以稍微不同的方式。你的方法是正確的，但你的翻譯是與教科書中提出的有點不同。生成我們需要從第一組的觀察中減去第一組的平均值，然後加上共同或合併的平均值. 對於第二組，我們做同樣的事情。
當您計算新變量的均值時，這變得更加直觀. 通過首先減去它們各自的組均值，變量以零為中心。通過添加整體平均值我們最終得到一個以整體平均值為中心的觀察樣本。換句話說，我們對觀測值進行了變換，使它們具有相同的均值，這也是兩組加在一起的總體均值，也就是.
# 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 值。希望這可以幫助！