為什麼對重採樣數據集的假設檢驗經常拒絕空值？
tl;dr：從在 null 下生成的數據集開始，我對帶有替換的案例進行了重新採樣，並對每個重新採樣的數據集進行了假設檢驗。這些假設檢驗在超過 5% 的時間內拒絕零。
在下面，非常簡單的模擬中，我生成數據集，並且我為每個模型擬合了一個簡單的 OLS 模型。然後，對於每個數據集，我通過使用替換對原始數據集的行進行重新採樣來生成 1000 個新數據集（Davison & Hinkley 的經典文本中專門描述的適用於線性回歸的算法）。對於其中的每一個，我都適合相同的 OLS 模型。最終，bootstrap 樣本中大約 16% 的假設檢驗拒絕了 null，而我們應該得到 5%（就像我們在原始數據集中所做的那樣）。
我懷疑這與重複觀察導致關聯膨脹有關，因此為了比較，我在下面的代碼中嘗試了其他兩種方法（已註釋掉）。在方法 2 中，我修復, 然後替換在原始數據集上使用來自 OLS 模型的重採樣殘差。在方法 3 中，我繪製了一個沒有放回的隨機子樣本。這兩種選擇都有效，即他們的假設檢驗在 5% 的時間裡拒絕了零。
**我的問題：重複觀察是罪魁禍首，我說得對嗎？**如果是這樣，鑑於這是一種標準的自舉方法，我們究竟在哪裡違反了標準的自舉理論？
更新 #1：更多模擬
我嘗試了一個更簡單的場景，一個只截取的回歸模型. 出現同樣的問題。
# note: simulation takes 5-10 min on my laptop; can reduce boot.reps # and n.sims.run if wanted # set the number of cores: can change this to match your machine library(doParallel) registerDoParallel(cores=8) boot.reps = 1000 n.sims.run = 1000 for ( j in 1:n.sims.run ) { # make initial dataset from which to bootstrap # generate under null d = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) ) # fit OLS to original data mod.orig = lm( Y1 ~ X1, data = d ) bhat = coef( mod.orig )[["X1"]] se = coef(summary(mod.orig))["X1",2] rej = coef(summary(mod.orig))["X1",4] < 0.05 # run all bootstrap iterates parallel.time = system.time( { r = foreach( icount( boot.reps ), .combine=rbind ) %dopar% { # Algorithm 6.2: Resample entire cases - FAILS # residuals of this model are repeated, so not normal? ids = sample( 1:nrow(d), replace=TRUE ) b = d[ ids, ] # # Method 2: Resample just the residuals themselves - WORKS # b = data.frame( X1 = d$X1, Y1 = sample(mod.orig$residuals, replace = TRUE) ) # # Method 3: Subsampling without replacement - WORKS # ids = sample( 1:nrow(d), size = 500, replace=FALSE ) # b = d[ ids, ] # save stats from bootstrap sample mod = lm( Y1 ~ X1, data = b ) data.frame( bhat = coef( mod )[["X1"]], se = coef(summary(mod))["X1",2], rej = coef(summary(mod))["X1",4] < 0.05 ) } } )[3] ###### Results for This Simulation Rep ##### r = data.frame(r) names(r) = c( "bhat.bt", "se.bt", "rej.bt" ) # return results of each bootstrap iterate new.rows = data.frame( bt.iterate = 1:boot.reps, bhat.bt = r$bhat.bt, se.bt = r$se.bt, rej.bt = r$rej.bt ) # along with results from original sample new.rows$bhat = bhat new.rows$se = se new.rows$rej = rej # add row to output file if ( j == 1 ) res = new.rows else res = rbind( res, new.rows ) # res should have boot.reps rows per "j" in the for-loop # simulation rep counter d$sim.rep = j } # end loop over j simulation reps ##### Analyze results ##### # dataset with only one row per simulation s = res[ res$bt.iterate == 1, ] # prob of rejecting within each resample # should be 0.05 mean(res$rej.bt); mean(s$rej)更新#2：答案
評論和答案中提出了幾種可能性，我做了更多的模擬來經驗性地測試它們。事實證明，JWalker 是正確的，問題在於我們需要通過原始數據的估計來集中引導統計數據，以便在下得到正確的抽樣分佈. 但是，我也認為 whuber 關於違反參數測試假設的評論也是正確的，儘管在這種情況下，當我們解決 JWalker 的問題時，我們實際上確實得到了名義上的誤報。
當您對空值重新採樣時，回歸係數的預期值為零。當您對某些觀察到的數據集重新採樣時，預期值是該數據的觀察係數。如果在對觀察到的數據重新採樣時 P <= 0.05，則不是 I 型錯誤。實際上，如果 P > 0.05，則屬於 II 型錯誤。
您可以通過計算 abs(b) 和 mean(P) 之間的相關性來獲得一些直覺。這是複制您所做的更簡單的代碼，並在一組模擬中計算 b 和“I 型”錯誤之間的相關性
boot.reps = 1000 n.sims.run = 10 n <- 1000 b <- matrix(NA, nrow=boot.reps, ncol=n.sims.run) p <- matrix(NA, nrow=boot.reps, ncol=n.sims.run) for(sim_j in 1:n.sims.run){ x <- rnorm(n) y <- rnorm(n) inc <- 1:n for(boot_i in 1:boot.reps){ fit <- lm(y[inc] ~ x[inc]) b[boot_i, sim_j] <- abs(coefficients(summary(fit))['x[inc]', 'Estimate']) p[boot_i, sim_j] <- coefficients(summary(fit))['x[inc]', 'Pr(>|t|)'] inc <- sample(1:n, replace=TRUE) } } # note this is not really a type I error but whatever type1 <- apply(p, 2, function(x) sum(x <= 0.05))/boot.reps # correlation between b and "type I" cor(b[1, ], type1)通過 grand_chat更新答案不是 P <= 0.05 的頻率 > 0.05 的原因。答案很簡單，正如我上面所說的——每個重採樣的平均值的期望值是原始的、觀察到的平均值。這是 bootstrap 的全部基礎，它的開發是為了在觀察到的平均值上生成標準誤差/置信限，而不是作為假設檢驗。由於期望不為零，當然“I 型錯誤”會大於 alpha。這就是為什麼係數的大小（離零有多遠）與“I 型錯誤”與 alpha 的偏差大小之間存在相關性的原因。