¿Por qué las pruebas de hipótesis de los conjuntos de datos remuestreados rechazan el nulo con demasiada frecuencia?

Question

tl;dr: Comenzando con un conjunto de datos generados bajo el nulo, remuestreé los casos con reemplazo y conduje una prueba de hipótesis en cada conjunto de datos remuestreado. Estas pruebas de hipótesis rechazan el nulo más del 5% de las veces.

En la siguiente simulación, muy simple, genero conjuntos de datos con $X \sim N(0,1) \amalg Y \sim N(0,1)$ y le pongo un simple modelo de OLS a cada uno. Luego, para cada conjunto de datos, genero 1000 nuevos conjuntos de datos remuestreando filas del conjunto de datos original con reemplazo (un algoritmo específicamente descrito en el texto clásico de Davison & Hinkley como apropiado para la regresión lineal). Para cada uno de ellos, me ajusto al mismo modelo de OLS. En última instancia, alrededor del 16% de las pruebas de hipótesis dentro de las muestras de bootstrap rechazan el nulo mientras que deberíamos obtener un 5% (como en los conjuntos de datos originales).

Sospeché que tiene algo que ver con las observaciones repetidas que causan asociaciones infladas, así que para comparar, probé otros dos enfoques en el código que sigue (comentado). En el método 2, fijé $X$ y luego reemplaza $Y$ con residuos remuestreados del modelo OLS en el conjunto de datos original. En el método 3, tomo una submuestra aleatoria sin reemplazo. Ambas alternativas funcionan, es decir, sus pruebas de hipótesis rechazan el nulo 5% de las veces.

Mi pregunta: ¿Tengo razón en que las observaciones repetidas son las culpables? Si es así, dado que es un enfoque estándar para el bootstrapping, ¿dónde exactamente estamos violando la teoría estándar de bootstrap?

Actualización

Intenté un escenario aún más simple, un modelo de regresión de sólo interceptación para $Y$ . El mismo problema ocurre.

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

Answer 1

Si toma una muestra con reemplazo de su muestra normal original, la muestra de bootstrap resultante no es normal . La distribución conjunta de la muestra bootstrap sigue una distribución de mezcla nudosa que es muy probable que contenga registros duplicados, mientras que los valores duplicados tienen probabilidad cero cuando se toman muestras iid de una distribución normal.

Como ejemplo simple, si su muestra original son dos observaciones de una distribución normal univariante, entonces una muestra bootstrap con reemplazo consistirá la mitad del tiempo en la muestra original, y la mitad del tiempo consistirá en uno de los valores originales, duplicados. Está claro que la varianza de la muestra de bootstrap será en promedio menor que la del original, de hecho será la mitad del original.

La principal consecuencia es que la inferencia que estás haciendo basada en la teoría normal devuelve el mal $p$ -Valores cuando se aplican a la muestra de bootstrap. En particular, la teoría normal da lugar a reglas de decisión anticonservadoras, porque su muestra de bootstrap producirá $t$ estadísticas cuyos denominadores son más pequeños de lo que cabría esperar en la teoría normal, debido a la presencia de duplicados. Como resultado, la prueba de hipótesis de la teoría normal termina rechazando la hipótesis nula más de lo esperado.