Usando bootstrap bajo H0 para realizar una prueba para la diferencia de dos medios: reemplazo dentro de los grupos o dentro de la muestra agrupada

Supongamos que tengo 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 donde g1 tiene 6 observaciones y g2 tiene 22 . ANOVA tradicional sugiere que los grupos tienen medios diferentes cuando el valor crítico se establece en 0.05 (el valor de p es 0.0044 ).

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

Dado que mi objetivo es comparar la diferencia de medias, estos datos muestreados pequeños y no balanceados pueden dar resultados inapropiados con el enfoque tradicional. Por lo tanto, quiero realizar pruebas de permutación y bootstrap.

PRUEBA DE PERMUTACION

La hipótesis nula (H0) establece que las medias del grupo son las mismas. Esta suposición en la prueba de permutación se justifica agrupando grupos en una muestra. Esto asegura que las muestras para dos grupos se extrajeron de la distribución idéntica. Mediante el muestreo repetido (o más precisamente, la reorganización) de los datos agrupados, las observaciones se reasignan (barajan) a las muestras de una manera nueva, y se calcula el estadístico de prueba. Realizando esto n veces, dará una distribución de muestreo de las estadísticas de prueba bajo el supuesto de que H0 es VERDADERO. Al final, bajo el valor H0, 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 informado de la prueba de permutación es 0.0053 . Bien, si lo hice correctamente, las permutaciones y el ANOVA paramétrico dan resultados casi idénticos.

OREJA

En primer lugar, soy consciente de que bootstrap no puede ayudar cuando los tamaños de muestra son demasiado pequeños. Esta publicación mostró que puede ser aún peor y engañoso . Además, el segundo destacó que la prueba de permutación es generalmente mejor que bootstrap cuando la prueba de hipótesis es el objetivo principal. Sin embargo, esta gran publicación aborda diferencias importantes entre los métodos intensivos en computadora. Sin embargo, aquí quiero plantear (creo) una pregunta diferente.

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

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 de bootstrap realizado de esta manera es 0.005 . Incluso si esto suena razonable y casi idéntico al ANOVA paramétrico y la prueba de permutación, ¿es apropiado justificar H0 en este arranque sobre la base de que acabamos de agrupar las muestras de las que extrajimos las muestras posteriores?

Enfoque diferente que encontré en varios artículos científicos. Específicamente, vi que los investigadores modifican los datos para cumplir con H0 antes del arranque. Al buscar, encontré publicaciones muy interesantes en CV donde @ jan.s explicó resultados inusuales de bootstrap en la pregunta de la publicación donde el objetivo era comparar dos medios. Sin embargo, en esa publicación no se explica cómo realizar el arranque cuando se modifican los datos antes del arranque. El enfoque en el que se modifican los datos antes del arranque se ve así:

1. H0 afirma que las medias de dos grupos son iguales
2. H0 es cierto si restamos observaciones individuales de la media de la muestra agrupada

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

3. Los datos modificados serán la base para un mayor arranque, con advertencias de que el muestreo se lleva a cabo dentro de cada grupo por separado .
4. La diferencia entre la media de bootstrapped de g1 y g2 se calcula y se compara con la diferencia observada (no modificada) entre los grupos.
5. La proporción de valores iguales o más extremos que los observados dividido por el número de iteraciones dará el valor p.

Aquí está el código (Bootstrap2: remuestreo dentro de los grupos después de la modificación de 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 arranque 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 la explicación de @ jan.s , pero no puedo entender dónde está la clave ...

Aquí está mi opinión al respecto, basada en el capítulo 16 de Una introducción a la rutina de Efron y Tibshirani (página 220-224). En resumen, su segundo algoritmo de arranque se implementó incorrectamente, pero la idea general es correcta.

Al realizar pruebas de bootstrap, uno debe asegurarse de que el método de remuestreo genera datos que corresponden a la hipótesis nula. Usaré los datos del sueño en R para ilustrar esta publicación. Tenga en cuenta que estoy usando la estadística de prueba estudiada en lugar de solo la diferencia de medias, que recomienda 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

$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}$$\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

Esta vez terminamos con valores p similares para los tres enfoques. ¡Espero que esto ayude!