Utilizar bootstrap bajo H0 para realizar una prueba de la diferencia de dos medias: sustitución dentro de los grupos o dentro de la muestra conjunta.

Question

Supongamos que tengo unos datos con dos grupos independientes:

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)

Es evidente que el tamaño de la muestra por grupo está sesgado cuando g1 tiene 6 observaciones y g2 tiene 22 . El ANOVA tradicional sugiere que los grupos tienen medias diferentes cuando el valor crítico se fija en 0,05 (el valor p es 0.0044 ).

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

Dado que mi objetivo es comparar la diferencia de medias, unos datos tan desequilibrados y con una muestra tan pequeña podrían dar resultados inadecuados con el enfoque tradicional. Por lo tanto, quiero realizar una prueba de permutación y bootstrap.

PRUEBA DE PERMUTACIÓN

La hipótesis nula (H0) afirma que las medias de los grupos son iguales. Esta hipótesis en la prueba de permutación se justifica agrupando los grupos en una sola muestra. Esto garantiza que las muestras de los dos grupos se extrajeron de la misma distribución. Mediante el muestreo repetido (o, más exactamente, la reorganización) de los datos agrupados, las observaciones se reasignan (barajan) a las muestras de una nueva forma y se calcula el estadístico de la prueba. Si se realiza esto n veces, se obtendrá la distribución de muestreo de los estadísticos de la prueba en el supuesto de que H0 sea VERDADERO. Al final, en la hipótesis H0, el valor p es la probabilidad de que el estadístico de prueba sea igual o superior al valor observado.

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)

El valor p de la prueba de permutación es 0.0053 . OK, si lo he hecho correctamente, las permutaciones y el ANOVA paramétrico dan resultados casi idénticos.

BOOTSTRAP

En primer lugar, soy consciente de que bootstrap no puede ayudar cuando el tamaño de la muestra es demasiado pequeño. Este post demostró que puede ser aún peor y engañoso . También, el segundo puso de relieve que la prueba de permutación suele ser mejor que la bootstrap cuando la comprobación de hipótesis es el objetivo principal. No obstante, este gran post aborda importantes diferencias entre los métodos informáticos intensivos. Sin embargo, aquí quiero plantear (creo) una cuestión diferente.

Permítanme presentarles primero el enfoque bootstrap más común (Bootstrap1: remuestreo dentro de la muestra conjunta ):

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)

El valor P del bootstrap realizado de esta forma es 0.005 . Aunque esto suene razonable y casi idéntico al ANOVA paramétrico y a la prueba de permutación, ¿es apropiado justificar H0 en este bootstrap sobre la base de que sólo muestras agrupadas de la que extrajimos las muestras posteriores?

Enfoque diferente que he encontrado en varios artículos científicos. En concreto, vi que los investigadores modifican los datos para cumplir H0 antes del bootstrap. Buscando por ahí, he encontrado muy interesante post en CV donde @jan.s explicó resultados inusuales de bootstrap en la pregunta posterior en la que el objetivo era comparar dos medias. Sin embargo, en ese post no se explica cómo realizar el bootstrap cuando los datos se modifican antes del bootstrap. El enfoque en el que los datos se modifican antes del bootstrap es el siguiente:

1. H0 afirma que las medias de dos grupos son iguales
2. H0 se cumple si restamos las observaciones individuales de la media de la muestra conjunta

En este caso, la modificación de los datos debería afectar a las medias de los grupos, y por tanto a su diferencia, pero no a la variación dentro de (y entre) los grupos.

3. Los datos modificados servirán de base para otros bootstrap, con la advertencia de que el muestreo se lleva a cabo dentro de cada grupo por separado .
4. Se calcula la diferencia entre la media bootstrapped de g1 y g2 y se compara con la diferencia observada (no modificada) entre grupos.
5. La proporción de valores extremos iguales o superiores al observado dividida por el número de iteraciones dará el valor p.

Aquí está el código (Bootstrap2: remuestreo dentro de los grupos tras modificar que H0 es VERDADERO ):

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)

Tal bootstrap realizado dará un valor p de 0.514 que es tremendamente diferente en comparación con las pruebas anteriores. Creo que esto tiene que ver con @jan.s explicación, pero no puedo averiguar dónde está la clave ...

Answer 1

Esta es mi opinión al respecto, basada en el capítulo 16 del libro de Efron y Tibshirani Introducción al bootstrap (página 220-224). En resumidas cuentas, tu segundo algoritmo bootstrap está mal implementado, pero la idea general es correcta.

Al realizar pruebas bootstrap, hay que asegurarse de que el método de remuestreo genera datos que corresponden a la hipótesis nula. Utilizaré los datos del sueño en R para ilustrar este post. Tenga en cuenta que estoy usando la estadística de prueba studentized en lugar de sólo la diferencia de medias, que es recomendado por el libro de texto.

La prueba t clásica, que utiliza un resultado analítico para obtener información sobre la distribución muestral del estadístico t, arroja el siguiente resultado:

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

Un enfoque es similar en espíritu a la prueba de permutación más conocida: se toman muestras de todo el conjunto de observaciones ignorando las etiquetas de agrupación. A continuación, la primera $n1$ se asignan al primer grupo y los restantes $n2$ al segundo grupo.

# 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

Sin embargo, este algoritmo en realidad está probando si la distribución de x e y son idénticas. Si simplemente nos interesa saber si sus medias poblacionales son iguales o no, sin hacer ninguna suposición sobre su varianza, deberíamos generar datos bajo $H_0$ de una manera ligeramente diferente. Ibas por buen camino con tu planteamiento, pero tu traducción a $H_0$ es un poco diferente de la propuesta en el libro de texto. Para generar $H_0$ tenemos que restar la media del primer grupo de las observaciones del primer grupo y luego sumar la media común o agrupada $\bar{z}$ . Para el segundo grupo hacemos lo mismo.

$$ \tilde{x}_i = x_i - \bar{x} + \bar{z} $$ $$ \tilde{y}_i = y_i - \bar{y} + \bar{z}$$

Esto resulta más intuitivo cuando se calculan las medias de las nuevas variables $\tilde{x}/\tilde{y}$ . Al restar primero las medias de sus respectivos grupos, las variables se centran en torno a cero. Sumando la media global $\bar{z}$ obtenemos una muestra de observaciones centrada en torno a la media global. En otras palabras, transformamos las observaciones para que tengan la misma media, que es también la media general de ambos grupos juntos, que es exactamente $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

Esta vez obtuvimos valores p similares para los tres enfoques.