Utilizo el paquete "boot" para calcular un valor p de arranque de 2 lados aproximado , pero el resultado está demasiado lejos del valor p de usar t.test. No puedo entender qué hice mal en mi código R. ¿Puede alguien darme una pista para esto?

time = c(14,18,11,13,18,17,21,9,16,17,14,15,
         12,12,14,13,6,18,14,16,10,7,15,10)
group=c(rep(1:2, each=12))
sleep = data.frame(time, group)

require(boot)
diff = function(d1,i){
    d = d1[i,]
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

set.seed(1234)
b3 = boot(data = sleep, statistic = diff, R = 5000, strata=sleep$group)

pvalue = mean(abs(b3$t) > abs(b3$t0))
pvalue 

El valor p de arranque de 2 lados (pvalue) = 0.4804 pero el valor p de 2 lados de t.test es 0.04342. Ambos valores p tienen una diferencia de alrededor de 11 veces. ¿Cómo puede pasar esto?

Está utilizando bootstrap para generar datos bajo la distribución empírica de los datos observados. Esto puede ser útil para dar un intervalo de confianza en la diferencia entre las dos medias:

> quantile(b3$t,c(0.025,0.975))
     2.5%     97.5% 
0.4166667 5.5833333 

Para obtener un valor , debe generar permutaciones bajo la hipótesis nula. Esto se puede hacer, por ejemplo, así:$p$

diff2 = function(d1,i){
    d = d1; 
    d$group <- d$group[i];  # randomly re-assign groups
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

> set.seed(1234)
> b4 = boot(data = sleep, statistic = diff2, R = 5000)
> mean(abs(b4$t) > abs(b4$t0))
[1] 0.046

En esta solución, el tamaño de los grupos no es fijo, usted reasigna aleatoriamente un grupo a cada individuo iniciando desde el conjunto de grupos inicial. Me parece legítimo, sin embargo, una solución más clásica es fijar el número de individuos de cada grupo, por lo que simplemente permutas los grupos en lugar de arrancar (esto generalmente está motivado por el diseño del experimento, donde los tamaños de los grupos se fijan de antemano ):

> R <- 10000; d <- sleep
> b5 <- numeric(R); for(i in 1:R) { 
+    d$group <- sample(d$group, length(d$group)); 
+    b5[i] <- mean(d$time[d$group==1])-mean(d$time[d$group==2]); 
+ }
> mean(abs(b5) > 3)
[1] 0.0372

La respuesta de Elvis se basa en permutaciones, pero en mi opinión no aclara qué está mal con el enfoque de arranque original. Permítanme discutir una solución basada únicamente en bootstrap.

El problema crucial de su simulación original es que bootstrap siempre le proporciona la distribución VERDADERA de la estadística de prueba. Sin embargo, al calcular el valor p, debe comparar el valor obtenido del estadístico de prueba con su distribución BAJO H0, es decir, ¡no con la distribución verdadera!

[Dejémoslo claro. Por ejemplo, se sabe que el estadístico de prueba T de la prueba t clásica tiene la distribución t "central" clásica bajo H0 y una distribución no central en general. Sin embargo, todos están familiarizados con el hecho de que el valor observado de T se compara con la distribución t "central" clásica, es decir, no se intenta obtener la distribución t [no central] verdadera para hacer la comparación con T.]

Su valor p 0.4804 es muy grande, porque el valor observado "t0" del estadístico de prueba Media [1] -Mean [2] se encuentra muy cerca del centro de la muestra "t". Es natural y, por lo general, siempre lo es [es decir, independientemente de la validez de H0], porque la muestra de arranque "t" emula la distribución ACTUAL de la media [1] - Media [2]. Pero, como se señaló anteriormente [y también por Elvis], lo que realmente necesita es la distribución de Mean [1] -Mean [2] BAJO H0. Eso es obvio

1) bajo H0, la distribución de Media [1] -Mean [2] se centrará alrededor de 0,

2) su forma no depende de la validez de H0.

Estos dos puntos implican que la distribución de la media [1] -Medio [2] bajo H0 puede ser emulada por la muestra inicializada "t" SHIFTED para que se centre alrededor de 0. En R:

b3.under.H0 <- b3$t - mean(b3$t)

y el valor p correspondiente será:

mean(abs(b3.under.H0) > abs(b3$t0))

lo que le da un valor "muy agradable" de 0.0232. :-)

Permítanme señalar que el punto "2)" mencionado anteriormente se llama "equivalencia de traducción" del estadístico de prueba y NO tiene que ser válido en general. Es decir, para algunas estadísticas de prueba, el desplazamiento de la "t" de arranque no le proporciona una estimación válida de la distribución de la estadística de prueba bajo HO! Eche un vistazo a esta discusión y especialmente a la respuesta de P. Dalgaard: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Su problema de prueba produce una distribución perfectamente simétrica del estadístico de prueba, pero tenga en cuenta que existen algunos problemas con la obtención de valores p de DOS CARAS en caso de una distribución de arranque sesgada del estadístico de prueba. Nuevamente, lea el enlace de arriba.

[Y finalmente, usaría la prueba de permutación "pura" en su situación; es decir, la segunda mitad de la respuesta de Elvis. :-)]

Existen numerosas formas de calcular los CI de arranque y los valores p. El problema principal es que es imposible que el bootstrap genere datos bajo una hipótesis nula. La prueba de permutación es una alternativa viable basada en remuestreo para esto. Para usar una rutina de arranque adecuada, debe hacer algunas suposiciones sobre la distribución de muestreo de la estadística de prueba.

$\beta^*_0 = \hat{\beta} - \hat{\beta}^*$$\beta^*_0 = \hat{\beta}^* - \hat{\beta}$

bootstrap normal

Un enfoque es un bootstrap normal donde se toma la media y la desviación estándar de la distribución de bootstrap, se calcula la distribución de muestreo debajo de la nula cambiando la distribución y utilizando los percentiles normales de la distribución nula en el punto de la estimación en la muestra de arranque original . Este es un enfoque razonable cuando la distribución de bootstrap es normal, la inspección visual generalmente es suficiente aquí. Los resultados que utilizan este enfoque suelen ser muy similares a la estimación de errores robusta o basada en sándwich, que es robusta frente a los supuestos de heterocedasticidad y / o varianza de muestra finita. La suposición de una estadística de prueba normal es una condición más fuerte de los supuestos en la próxima prueba de arranque que discutiré.

bootstrap percentil

$\mathcal{F}_0^*$$2 \times \min (\mathcal{F}_0^*(\hat{\beta}), 1-\mathcal{F}_0^*(\hat{\beta}))$

Bootstrap Studentizado

$p$

Ejemplo de programación

Como ejemplo, usaré los citydatos en el paquete bootstrap. Los intervalos de confianza de bootstrap se calculan con este código:

ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.ci(city.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

y producir esta salida:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm", 
    "basic", "perc", "bca"))

Intervals : 
Level      Normal              Basic         
90%   ( 1.111,  1.837 )   ( 1.030,  1.750 )   
95%   ( 1.042,  1.906 )   ( 0.895,  1.790 )  

Level     Percentile            BCa          
90%   ( 1.291,  2.011 )   ( 1.292,  2.023 )   
95%   ( 1.251,  2.146 )   ( 1.255,  2.155 )  
Calculations and Intervals on Original Scale

El IC del 95% para el bootstrap normal se obtiene calculando:

with(city.boot, 2*t0 - mean(t) + qnorm(c(0.025, 0.975)) %o% sqrt(var(t)[1,1]))

El valor p se obtiene así:

> with(city.boot, pnorm(abs((2*t0 - mean(t) - 1) / sqrt(var(t)[1,1])), lower.tail=F)*2)
[1] 0.0315

Lo que concuerda en que el IC normal del 95% no incluye el valor de la relación nula de 1.

Se obtiene el IC percentil (con algunas diferencias debido a los métodos para los lazos):

quantile(city.boot$t, c(0.025, 0.975))

Y el valor p para el percentil bootstrap es:

cvs <- quantile(city.boot$t0 - city.boot$t + 1, c(0.025, 0.975))
mean(city.boot$t > cvs[1] & city.boot$t < cvs[2])

Da una p de 0.035 que también concuerda con el intervalo de confianza en términos de exclusión de 1 del valor. En general, no podemos observar que, mientras que el ancho del IC del percentil es casi tan ancho como el IC normal y que el IC del percentil está más lejos de ser nulo, el IC del percentil debería proporcionar valores p más bajos. Esto se debe a que la forma de la distribución de muestreo subyacente al IC para el método de percentil no es normal.