tl; dr: Zaczynając od zestawu danych wygenerowanego pod wartością zerową, ponownie próbkowałem przypadki z zamianą i przeprowadzałem test hipotez na każdym zestawie danych ponownie próbkowanym. Te testy hipotez odrzucają zero w ponad 5% przypadków.
W poniższej, bardzo prostej symulacji, generuję zestawy danych za pomocą i dopasowuję do nich prosty model OLS. Następnie dla każdego zestawu danych generuję 1000 nowych zestawów danych poprzez ponowne próbkowanie wierszy oryginalnego zestawu danych z zastąpieniem (algorytm opisany szczegółowo w klasycznym tekście Davisona i Hinkleya jako odpowiedni dla regresji liniowej). Do każdego z nich pasuję do tego samego modelu OLS. Ostatecznie około 16% testów hipotez w próbkach bootstrap odrzuca wartość zerową , podczas gdy powinniśmy otrzymać 5% (tak jak robimy to w oryginalnych zestawach danych).
Podejrzewałem, że ma to coś wspólnego z powtarzającymi się obserwacjami powodującymi zawyżone skojarzenia, więc dla porównania wypróbowałem dwa inne podejścia w poniższym kodzie (skomentowałem). W metodzie 2 naprawiam , a następnie zamieniam resztki próbek z modelu OLS w oryginalnym zestawie danych. W metodzie 3 rysuję losową podpróbkę bez zamiany. Obie te alternatywy działają, tj. Ich testy hipotez odrzucają zerowy 5% czasu.
Moje pytanie: czy mam rację, że winowajcami są powtarzające się obserwacje? Jeśli tak, biorąc pod uwagę, że jest to standardowe podejście do ładowania, to gdzie dokładnie naruszamy standardową teorię ładowania?
Aktualizacja nr 1: Więcej symulacji
Próbowałem jeszcze prostszy scenariusz, model regresji przechwycenia tylko dla . Ten sam problem występuje.
# 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)Aktualizacja nr 2: odpowiedź
W komentarzach i odpowiedziach zaproponowano kilka możliwości i wykonałem więcej symulacji, aby je empirycznie przetestować. Okazuje się, że JWalker ma rację, że problem polega na tym, że musieliśmy wyśrodkować statystyki bootstrapu według szacunków oryginalnych danych, aby uzyskać prawidłowy rozkład próbkowania w . Myślę jednak również, że komentarz Whubera dotyczący naruszenia założeń testu parametrycznego jest również poprawny, chociaż w tym przypadku faktycznie otrzymujemy nominalne fałszywe alarmy, gdy naprawimy problem JWalkera.
idsids <- unique(ids)