Conocí un comportamiento paradójico de las llamadas "pruebas exactas" o "pruebas de permutación", cuyo prototipo es la prueba de Fisher. Aquí está.

Imagine que tiene dos grupos de 400 individuos (p. Ej., 400 casos de control frente a 400) y una covariable con dos modalidades (p. Ej., Expuestos / no expuestos). Solo hay 5 individuos expuestos, todos en el segundo grupo. La prueba de Fisher es así:

> x <- matrix( c(400, 395, 0, 5) , ncol = 2)
> x
     [,1] [,2]
[1,]  400    0
[2,]  395    5
> fisher.test(x)

    Fisher's Exact Test for Count Data

data:  x
p-value = 0.06172
(...)

Pero ahora, hay cierta heterogeneidad en el segundo grupo (los casos), por ejemplo, la forma de la enfermedad o el centro de reclutamiento. Se puede dividir en 4 grupos de 100 individuos. Es probable que ocurra algo como esto:

> x <- matrix( c(400, 99, 99 , 99, 98, 0, 1, 1, 1, 2) , ncol = 2)
> x
     [,1] [,2]
[1,]  400    0
[2,]   99    1
[3,]   99    1
[4,]   99    1
[5,]   98    2
> fisher.test(x)

    Fisher's Exact Test for Count Data

data:  x 
p-value = 0.03319
alternative hypothesis: two.sided
(...)

Ahora, tenemos ...$p < 0.05$

Este es solo un ejemplo. Pero podemos simular el poder de las dos estrategias de análisis, suponiendo que en los primeros 400 individuos, la frecuencia de exposición es 0, y que es 0.0125 en los 400 individuos restantes.

Podemos estimar el poder del análisis con dos grupos de 400 individuos:

> p1 <- replicate(1000, { n <- rbinom(1, 400, 0.0125); 
                          x <- matrix( c(400, 400 - n, 0, n), ncol = 2); 
                          fisher.test(x)$p.value} )
> mean(p1 < 0.05)
[1] 0.372

Y con un grupo de 400 y 4 grupos de 100 individuos:

> p2 <- replicate(1000, { n <- rbinom(4, 100, 0.0125); 
                          x <- matrix( c(400, 100 - n, 0, n), ncol = 2);
                          fisher.test(x)$p.value} )
> mean(p2 < 0.05)
[1] 0.629

Hay una gran diferencia de poder. Dividir los casos en 4 subgrupos da una prueba más poderosa, incluso si no hay diferencia de distribución entre estos subgrupos. Por supuesto, esta ganancia de poder no es atribuible a una mayor tasa de error tipo I.

¿Es este fenómeno bien conocido? ¿Significa eso que la primera estrategia tiene poca potencia? ¿Sería un valor p bootstrapped una mejor solución? Todos sus comentarios son bienvenidos.

Post Scriptum

Como señaló @MartijnWeterings, una gran parte de la razón de este comportamiento (¡que no es exactamente mi pregunta!) Radica en el hecho de que los verdaderos errores de tipo I de las estrategias de análisis de remolque no son los mismos. Sin embargo, esto no parece explicar todo. Traté de comparar las curvas ROC para vs .H 1 : p 0 = 0.05 ≠ p 1 = $H_0 : p_0 = p_1 = 0.005$$H1 : p_0 = 0.05 \ne p1 = 0.0125$

Aquí está mi código.

B <- 1e5
p0 <- 0.005
p1 <- 0.0125

# simulation under H0 with p = p0 = 0.005 in all groups
# a = 2 groups 400:400, b = 5 groupe 400:100:100:100:100

p.H0.a <- replicate(B, { n <- rbinom( 2, c(400,400), p0);
                           x <- matrix( c( c(400,400) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

p.H0.b <- replicate(B, { n <- rbinom( 5, c(400,rep(100,4)), p0);
                           x <- matrix( c( c(400,rep(100,4)) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

# simulation under H1 with p0 = 0.005 (controls) and p1 = 0.0125 (cases)

p.H1.a <- replicate(B, { n <- rbinom( 2, c(400,400), c(p0,p1) );
                           x <- matrix( c( c(400,400) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

p.H1.b <- replicate(B, { n <- rbinom( 5, c(400,rep(100,4)), c(p0,rep(p1,4)) );
                           x <- matrix( c( c(400,rep(100,4)) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

# roc curve 

ROC <- function(p.H0, p.H1) {
  p.threshold <- seq(0, 1.001, length=501)
  alpha <- sapply(p.threshold, function(th) mean(p.H0 <= th) )
  power <- sapply(p.threshold, function(th) mean(p.H1 <= th) )
  list(x = alpha, y = power)
}

par(mfrow=c(1,2))
plot( ROC(p.H0.a, p.H1.a) , type="b", xlab = "alpha", ylab = "1-beta" , xlim=c(0,1), ylim=c(0,1), asp = 1)
lines( ROC(p.H0.b, p.H1.b) , col="red", type="b" )
abline(0,1)

plot( ROC(p.H0.a, p.H1.a) , type="b", xlab = "alpha", ylab = "1-beta" , xlim=c(0,.1) )
lines( ROC(p.H0.b, p.H1.b) , col="red", type="b" )
abline(0,1)

Aquí está el resultado:

Entonces vemos que una comparación con el mismo error verdadero de tipo I todavía conduce a diferencias (de hecho mucho más pequeñas).

Por qué los valores p son diferentes

Hay dos efectos en curso:

• $\chi^2$

  $\chi^2$

  Es por eso que el valor p difiere en casi un factor 2. (no exactamente por el siguiente punto)

• Mientras pierde el 5 0 0 0 0 como un caso igualmente extremo, gana el 1 4 0 0 0 como un caso más extremo que 0 2 1 1 1.

$\chi^2$

ejemplo de código:

# probability of distribution a and b exposures among 2 groups of 400
draw2 <- function(a,b) {
  choose(400,a)*choose(400,b)/choose(800,5)
}

# probability of distribution a, b, c, d and e exposures among 5 groups of resp 400, 100, 100, 100, 100
draw5 <- function(a,b,c,d,e) {
choose(400,a)*choose(100,b)*choose(100,c)*choose(100,d)*choose(100,e)/choose(800,5)
}

# looping all possible distributions of 5 exposers among 5 groups
# summing the probability when it's p-value is smaller or equal to the observed value 0 2 1 1 1
sumx <- 0
for (f in c(0:5)) {
  for(g in c(0:(5-f))) {
    for(h in c(0:(5-f-g))) {
      for(i in c(0:(5-f-g-h))) {
        j = 5-f-g-h-i
        if (draw5(f, g, h, i, j) <= draw5(0, 2, 1, 1, 1)) {
          sumx <- sumx + draw5(f, g, h, i, j)
        }
      }
    }
  } 
}
sumx  #output is 0.3318617

# the split up case (5 groups, 400 100 100 100 100) can be calculated manually
# as a sum of probabilities for cases 0 5 and 1 4 0 0 0 (0 5 includes all cases 1 a b c d with the sum of the latter four equal to 5)
fisher.test(matrix( c(400, 98, 99 , 99, 99, 0, 2, 1, 1, 1) , ncol = 2))[1]
draw2(0,5) + 4*draw(1,4,0,0,0)

# the original case of 2 groups (400 400) can be calculated manually
# as a sum of probabilities for the cases 0 5 and 5 0 
fisher.test(matrix( c(400, 395, 0, 5) , ncol = 2))[1]
draw2(0,5) + draw2(5,0)

salida de ese último bit

> fisher.test(matrix( c(400, 98, 99 , 99, 99, 0, 2, 1, 1, 1) , ncol = 2))[1]
$p.value
[1] 0.03318617

> draw2(0,5) + 4*draw(1,4,0,0,0)
[1] 0.03318617

> fisher.test(matrix( c(400, 395, 0, 5) , ncol = 2))[1]
$p.value
[1] 0.06171924

> draw2(0,5) + draw2(5,0)
[1] 0.06171924

Cómo afecta el poder al dividir grupos

• Existen algunas diferencias debido a los pasos discretos en los niveles 'disponibles' de los valores p y la conservaduría de la prueba exacta de Fishers (y estas diferencias pueden llegar a ser bastante grandes).

• también la prueba de Fisher se ajusta al modelo (desconocido) basado en los datos y luego usa este modelo para calcular los valores de p. El modelo en el ejemplo es que hay exactamente 5 individuos expuestos. Si modela los datos con un binomio para los diferentes grupos, ocasionalmente obtendrá más o menos de 5 individuos. Cuando aplique la prueba de Fisher a esto, se ajustará parte del error y los residuos serán más pequeños en comparación con las pruebas con márgenes fijos. El resultado es que la prueba es demasiado conservadora, no exacta.

$\alpha$

$H_0$$H_0$$H_a$

La pregunta sigue siendo si esto es válido para todas las situaciones posibles.

Ajuste de código 3 veces de su análisis de potencia (y 3 imágenes):

usando binomio restringiendo al caso de 5 individuos expuestos

$H_0$

Es interesante ver que el efecto es mucho más fuerte para el caso 400-400 (rojo) frente al caso 400-100-100-100-100 (azul). Por lo tanto, es posible que usemos esta división para aumentar la potencia y aumentar la probabilidad de rechazar el H_0. (aunque no nos importa tanto hacer que el error tipo I sea más probable, entonces el punto de hacer esta división para aumentar la potencia puede no ser siempre tan fuerte)

usando binomial no restringiendo a 5 individuos expuestos

Si usamos un binomio como usted, ninguno de los dos casos 400-400 (rojo) o 400-100-100-100-100 (azul) da un valor p exacto. Esto se debe a que la prueba exacta de Fisher supone totales fijos de fila y columna, pero el modelo binomial permite que estos sean libres. La prueba de Fisher 'ajustará' los totales de fila y columna, haciendo que el término residual sea más pequeño que el término de error verdadero.

¿El aumento de poder tiene un costo?

$H_0$$H_a$$H_a$

# using binomial distribution for 400, 100, 100, 100, 100
# x uses separate cases
# y uses the sum of the 100 groups
p <- replicate(4000, { n <- rbinom(4, 100, 0.006125); m <- rbinom(1, 400, 0.006125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:1000)/1000
m1 <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

plot(ps,ps,type="l",
     xlab = "chosen alpha level",
     ylab = "p rejection")
lines(ps,m1,col=4)
lines(ps,m2,col=2)

title("due to concervative test p-value will be smaller\n leading to differences")

# using all samples also when the sum exposed individuals is not 5
ps <- c(1:1000)/1000
m1 <- sapply(ps,FUN = function(x) mean(p[2,] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,] < x))

plot(ps,ps,type="l", 
     xlab = "chosen alpha level",
     ylab = "p rejection")
lines(ps,m1,col=4)
lines(ps,m2,col=2)

title("overly conservative, low effective p-values \n fitting marginals makes residuals smaller than real error")


#   
# Third graph comparing H_0 and H_a
#
# using binomial distribution for 400, 100, 100, 100, 100
# x uses separate cases
# y uses the sum of the 100 groups
offset <- 0.5
p <- replicate(10000, { n <- rbinom(4, 100, offset*0.0125); m <- rbinom(1, 400, (1-offset)*0.0125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:10000)/10000
m1 <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

offset <- 0.6
p <- replicate(10000, { n <- rbinom(4, 100, offset*0.0125); m <- rbinom(1, 400, (1-offset)*0.0125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:10000)/10000
m1a <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2a <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

plot(ps,ps,type="l",
     xlab = "p rejecting if H_0 true",
     ylab = "p rejecting if H_a true",log="xy")
points(m1,m1a,col=4)
points(m2,m2a,col=2)

legend(0.01,0.001,c("400-400","400-100-100-100-100"),pch=c(1,1),col=c(2,4))

title("comparing H_0:p=0.5 \n with H_a:p=0.6")

¿Por qué afecta el poder?

Creo que la clave del problema está en la diferencia de los valores de resultado que se eligen como "significativos". La situación es de cinco individuos expuestos provenientes de 5 grupos de tamaño 400, 100, 100, 100 y 100. Se pueden hacer diferentes selecciones que se consideran "extremas". aparentemente el poder aumenta (incluso cuando el error efectivo de tipo I es el mismo) cuando elegimos la segunda estrategia.

Si bosquejáramos la diferencia entre la primera y la segunda estrategia gráficamente. Luego imagino un sistema de coordenadas con 5 ejes (para los grupos de 400 100 100 100 y 100) con un punto para los valores de hipótesis y la superficie que representa una distancia de desviación más allá de la cual la probabilidad está por debajo de cierto nivel. Con la primera estrategia, esta superficie es un cilindro, con la segunda estrategia, esta superficie es una esfera. Lo mismo es cierto para los valores verdaderos y una superficie a su alrededor para el error. Lo que queremos es que la superposición sea lo más pequeña posible.

Podemos hacer un gráfico real cuando consideramos un problema ligeramente diferente (con menor dimensionalidad).

$H_0: p=0.5$

El gráfico muestra cómo se distribuyen los grupos de 500 y 500 (en lugar de un solo grupo de 1000).

La prueba de hipótesis estándar evaluaría (para un nivel alfa del 95%) si la suma de X e Y es mayor que 531 o menor que 469.

Pero esto incluye una distribución desigual muy poco probable de X e Y.

$H_0$$H_a$

Sin embargo, esto no es (necesariamente) cierto cuando no seleccionamos la división de los grupos al azar y cuando puede haber un significado para los grupos.