Pourquoi les tests d'hypothèses sur les jeux de données rééchantillonnés rejettent-ils trop souvent le null?

tl; dr: En commençant par un ensemble de données généré sous la valeur nulle, j'ai rééchantillonné les cas avec remplacement et effectué un test d'hypothèse sur chaque ensemble de données rééchantillonné. Ces tests d'hypothèse rejettent le nul plus de 5% du temps.

Dans la simulation ci-dessous, très simple, je génère des ensembles de données avec , et un modèle OLS simple à chacun. Ensuite, pour chaque ensemble de données, je génère 1000 nouveaux ensembles de données en rééchantillonnant les lignes de l'ensemble de données d'origine avec remplacement (un algorithme spécifiquement décrit dans le texte classique de Davison & Hinkley comme étant approprié pour la régression linéaire). Pour chacun d'eux, je correspond au même modèle OLS. En fin de compte, environ 16% des tests d'hypothèse dans les échantillons de bootstrap rejettent le zéro , alors que nous devrions obtenir 5% (comme nous le faisons dans les ensembles de données d'origine).$X \sim N(0,1) \amalg Y \sim N(0,1)$

Je soupçonnais que cela avait quelque chose à voir avec des observations répétées provoquant des associations gonflées, donc à titre de comparaison, j'ai essayé deux autres approches dans le code ci-dessous (commenté). Dans la méthode 2, je corrige , puis je remplace par des résidus rééchantillonnés du modèle OLS sur le jeu de données d'origine. Dans la méthode 3, je dessine un sous-échantillon aléatoire sans remplacement. Ces deux alternatives fonctionnent, c'est-à-dire que leurs tests d'hypothèse rejettent les 5% nuls du temps.$X$$Y$

Ma question: ai-je raison de dire que les observations répétées sont le coupable? Si oui, étant donné qu'il s'agit d'une approche standard du bootstrap, où violons-nous exactement la théorie du bootstrap standard?

Mise à jour # 1: plus de simulations

J'ai essayé un scénario encore plus simple, un modèle de régression ordonnée à l' origine uniquement pour . Le même problème se produit.$Y$

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

Mise à jour # 2: la réponse

Plusieurs possibilités ont été proposées dans les commentaires et réponses, et j'ai fait plus de simulations pour les tester empiriquement. Il s'avère que JWalker a raison que le problème est que nous avions besoin de centrer les statistiques de bootstrap par l'estimation des données d'origine afin d'obtenir la distribution d'échantillonnage correcte sous . Cependant, je pense également que le commentaire de Whuber sur la violation des hypothèses de test paramétrique est également correct, bien que dans ce cas, nous obtenions en fait des faux positifs nominaux lorsque nous résolvons le problème de JWalker.$H_0$

Lorsque vous rééchantillonnez le null, la valeur attendue du coefficient de régression est zéro. Lorsque vous rééchantillonnez un ensemble de données observé, la valeur attendue est le coefficient observé pour ces données. Ce n'est pas une erreur de type I si P <= 0,05 lorsque vous rééchantillonnez les données observées. En fait, c'est une erreur de type II si vous avez P> 0,05.

Vous pouvez acquérir une certaine intuition en calculant la corrélation entre les abs (b) et la moyenne (P). Voici un code plus simple pour reproduire ce que vous avez fait et calculer la corrélation entre b et l'erreur "type I" sur l'ensemble des simulations

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)

mettre à jour la réponse par grand_chat n'est pas la raison pour laquelle la fréquence de P <= 0,05 est> 0,05. La réponse est très simple et ce que j'ai dit ci-dessus - la valeur attendue de la moyenne de chaque rééchantillonnage est la moyenne originale observée. C'est toute la base du bootstrap, qui a été développé pour générer des erreurs standard / limites de confiance sur une moyenne observée et non pas comme un test d'hypothèse. Étant donné que l'attente n'est pas nulle, bien sûr, "l'erreur de type I" sera supérieure à alpha. Et c'est pourquoi il y aura une corrélation entre la magnitude du coefficient (à quelle distance de zéro) et la magnitude de l'écart de "l'erreur de type I" par rapport à alpha.

Si vous échantillonnez avec remplacement à partir de votre échantillon normal d'origine, l'échantillon d'amorçage résultant n'est pas normal . La distribution conjointe de l'échantillon bootstrap suit une distribution de mélange noueux qui est très susceptible de contenir des enregistrements en double, tandis que les valeurs en double ont une probabilité nulle lorsque vous prenez des échantillons iid d'une distribution normale.

À titre d'exemple simple, si votre échantillon d'origine est constitué de deux observations d'une distribution normale univariée, un échantillon bootstrap avec remplacement sera composé pour moitié de l'échantillon d'origine et pour moitié de l'une des valeurs d'origine, dupliquées. Il est clair que la variance de l'échantillon de l'échantillon de bootstrap sera en moyenne inférieure à celle de l'original - en fait, elle sera la moitié de l'original.

$p$$t$

Je suis totalement d'accord avec la réponse de @ JWalker.

Il y a un autre aspect de ce problème. C'est dans votre processus de rééchantillonnage. Vous vous attendez à ce que le coefficient de régression soit centré autour de zéro parce que vous supposez Xet Yêtes indépendant. Cependant, dans votre rééchantillonnage, vous

ids = sample( 1:nrow(d), replace=TRUE )
  b = d[ ids, ]

ce qui crée une corrélation parce que vous échantillonnez Xet Yensemble. Par exemple, supposons que la première ligne de l'ensemble de données dsoit (x1, y1), Dans l'ensemble de données rééchantillonné P(Y = y1|X = x1) = 1, alors que si Xet Ysont indépendants, P(Y|X = x1)suit une distribution normale.

Donc, une autre façon de résoudre ce problème est d'utiliser

b = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

le même code que vous utilisez pour générer d, afin de rendre X et Y indépendants l'un de l'autre.

La même raison explique pourquoi il fonctionne avec un rééchantillonnage résiduel (car il Xest indépendant du nouveau Y).

Le plus gros problème ici est que les résultats du modèle sont parasites et donc très instables, car le modèle correspond juste au bruit. Dans un sens très littéral. Y1 n'est pas une variable dépendante en raison de la façon dont les données d'échantillon ont été générées.

Modifier, en réponse aux commentaires. Permettez-moi d'essayer à nouveau d'expliquer ma pensée.

Avec un OLS, l'intention générale est de découvrir et de quantifier les relations sous-jacentes dans les données. Avec des données réelles, nous ne les connaissons généralement pas exactement.

Mais c'est une situation de test artificielle. Nous connaissons le mécanisme de génération de données EXACT, il est juste là dans le code affiché par l'OP C'est

X1 = rnorm (n = 1000), Y1 = rnorm (n = 1000)

Si nous exprimons cela sous la forme familière d'une régression OLS, c'est-à-dire

Y1 = interception + Beta1 * X1 + Erreur
qui devient
Y1 = moyenne (X1) + 0 (X1) + Erreur

Donc, dans mon esprit, c'est un modèle exprimé sous forme linéaire, mais ce n'est PAS en fait une relation / modèle linéaire, car il n'y a pas de pente. Beta1 = 0,000000.

Lorsque nous générons les 1000 points de données aléatoires, le nuage de points va ressembler au spray circulaire de fusil de chasse classique. Il pourrait y avoir une certaine corrélation entre X1 et Y1 dans l'échantillon spécifique de 1000 points aléatoires qui a été généré, mais si c'est le cas, c'est un hasard. Si l'OLS trouve une corrélation, c'est-à-dire qu'il rejette l'hypothèse nulle selon laquelle il n'y a pas de pente, puisque nous savons définitivement qu'il n'y a vraiment aucune relation entre ces deux variables, l'OLS a littéralement trouvé un modèle dans le composant d'erreur. J'ai qualifié cela de «convenable au bruit» et de «faux».

De plus, l'une des hypothèses / exigences std d'un OLS est que (+/-) "le modèle de régression linéaire est" linéaire en paramètres ". Compte tenu des données, je pense que nous ne satisfaisons pas à cette hypothèse. Par conséquent, les statistiques sous-jacentes du test de signification sont inexactes. Je pense que la violation de l'hypothèse de linéarité est la cause directe des résultats non intuitifs du bootstrap.

Lorsque j'ai lu ce problème pour la première fois, il n'a pas sombré dans la mesure où l'OP avait l'intention de tester sous l'hypothèse nulle.

Mais les mêmes résultats non intuitifs se produiraient-ils si l'ensemble de données avait été généré

X1 = rnorm (n = 1000), Y1 = X1 * .4 + rnorm (n = 1000)?