Utilisation du bootstrap sous H0 pour effectuer un test de la différence de deux moyennes: remplacement au sein des groupes ou au sein de l'échantillon groupé

Supposons que j'ai des données avec deux groupes indépendants:

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)

Il est évident que la taille de l'échantillon par groupe est biaisée lorsque g1 a 6 observations et g2 en a 22 . L'ANOVA traditionnelle suggère que les groupes ont des moyens différents lorsque la valeur critique est définie sur 0,05 (la valeur p est de 0,0044 ).

summary (aov (lengths~group, data = lengths))  

Étant donné que mon objectif est de comparer la différence moyenne, de telles données non équilibrées et de petits échantillons pourraient donner des résultats inappropriés avec l'approche traditionnelle. Par conséquent, je veux effectuer un test de permutation et un bootstrap.

ESSAI DE PERMUTATION

L'hypothèse nulle (H0) indique que les moyennes du groupe sont les mêmes. Cette hypothèse dans le test de permutation est justifiée par la mise en commun des groupes en un seul échantillon. Cela garantit que les échantillons de deux groupes ont été tirés de la même distribution. Par échantillonnage répété (ou plus précisément - remaniement) à partir des données regroupées, les observations sont réallouées (mélangées) aux échantillons d'une nouvelle manière, et la statistique de test est calculée. L'exécution de cette n fois, donnera la distribution d'échantillonnage des statistiques de test dans l'hypothèse où H0 est VRAI. À la fin, sous la valeur H0, p est la probabilité que la statistique de test soit égale ou supérieure à la valeur observée.

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)

La valeur p rapportée du test de permutation est de 0,0053 . OK, si je l'ai fait correctement, les permutations et l'ANOVA paramétrique donnent des résultats presque identiques.

AMORCER

Tout d'abord, je suis conscient que le bootstrap ne peut pas aider lorsque la taille des échantillons est trop petite. Ce message a montré que cela peut être encore pire et trompeur . De plus, le second a souligné que le test de permutation est généralement meilleur que le bootstrap lorsque le test d'hypothèse est l'objectif principal. Néanmoins, ce grand article traite des différences importantes entre les méthodes à forte intensité informatique. Cependant, ici, je veux soulever (je crois) une question différente.

Permettez-moi de présenter d'abord l'approche de bootstrap la plus courante (Bootstrap1: rééchantillonnage au sein de l'échantillon regroupé ):

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)

La valeur P du bootstrap effectué de cette manière est de 0,005 . Même si cela semble raisonnable et presque identique à l'ANOVA paramétrique et au test de permutation, est-il approprié de justifier H0 dans ce bootstrap sur la base que nous venons de regrouper les échantillons dont nous avons tiré les échantillons suivants?

Approche différente que j'ai trouvée dans plusieurs articles scientifiques. Plus précisément, j'ai vu que les chercheurs modifient les données afin de rencontrer H0 avant le bootstrap. En cherchant autour, j'ai trouvé un post très intéressant dans CV où @ jan.s a expliqué des résultats inhabituels de bootstrap dans la question du post où le but était de comparer deux moyennes. Cependant, dans cet article, il n'est pas expliqué comment effectuer un amorçage lorsque les données sont modifiées avant l'amorçage. L'approche où les données sont modifiées avant le bootstrap ressemble à ceci:

1. H0 indique que les moyennes de deux groupes sont les mêmes
2. H0 est vrai si l'on soustrait les observations individuelles de la moyenne de l'échantillon groupé

Dans ce cas, la modification des données devrait affecter les moyennes des groupes, et donc sa différence, mais pas la variation au sein (et entre) les groupes.

3. Des données modifiées serviront de base à d'autres bootstrap, avec des avertissements que l'échantillonnage est effectué séparément dans chaque groupe .
4. La différence entre la moyenne bootstrappée de g1 et g2 est calculée et comparée à la différence observée (non modifiée) entre les groupes.
5. La proportion de valeurs égales ou plus extrêmes que celle observée, divisée par le nombre d'itérations, donnera la valeur p.

Voici le code (Bootstrap2: rééchantillonnage au sein des groupes après modification que H0 est VRAI ):

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)

Un tel bootstrap effectué donnera une valeur de p de 0,514, ce qui est extrêmement différent des tests précédents. Je crois que cela doit traiter l' explication de @ jan.s , mais je ne peux pas comprendre où est la clé ...

Voici mon point de vue, basé sur le chapitre 16 de An Introduction to the bootstrap d'Efron et Tibshirani (page 220-224). En bref, votre deuxième algorithme d'amorçage a été mal mis en œuvre, mais l'idée générale est correcte.

Lors des tests de bootstrap, il faut s'assurer que la méthode de ré-échantillonnage génère des données qui correspondent à l'hypothèse nulle. Je vais utiliser les données de sommeil dans R pour illustrer ce post. Notez que j'utilise la statistique de test studentisée plutôt que simplement la différence de moyenne, qui est recommandée par le manuel.

Le test t classique, qui utilise un résultat analytique pour obtenir des informations sur la distribution d'échantillonnage de la statistique t, donne le résultat suivant:

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

$n1$$n2$

# 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

$H_0$$H_0$$H_0$$\bar{z}$

$\tilde{x}/\tilde{y}$ H$\bar{z}$$H_0$

# 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

Cette fois-ci, nous nous sommes retrouvés avec des valeurs de p similaires pour les trois approches. J'espère que cela t'aides!