Valores p iguales a 0 en la prueba de permutación

Tengo dos conjuntos de datos y me gustaría saber si son significativamente diferentes o no (esto proviene de "¿ Dos grupos son significativamente diferentes? Prueba de uso ").

Decidí usar una prueba de permutación, haciendo lo siguiente en R:

permutation.test <- function(coding, lncrna) {
    coding <- coding[,1] # dataset1
    lncrna <- lncrna[,1] # dataset2

    ### Under null hyphotesis, both datasets would be the same. So:
    d <- c(coding, lncrna)

    # Observed difference
    diff.observed = mean(coding) - mean(lncrna)
    number_of_permutations = 5000
    diff.random = NULL

    for (i in 1:number_of_permutations) {
        # Sample from the combined dataset
        a.random = sample (d, length(coding), TRUE)
        b.random = sample (d, length(lncrna), TRUE)
        # Null (permuated) difference
        diff.random[i] = mean(b.random) - mean(a.random)
    }

    # P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    pvalue
}

Sin embargo, los valores p no deberían ser 0 según este documento: http://www.statsci.org/smyth/pubs/permp.pdf

¿Qué me recomiendan hacer? Es esta forma de calcular el valor p:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

¿una buena manera? ¿O es mejor hacer lo siguiente?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1

Discusión

Una prueba de permutación genera todas las permutaciones relevantes de un conjunto de datos, calcula un estadístico de prueba designado para cada permutación y evalúa el estadístico de prueba real en el contexto de la distribución de permutación resultante de las estadísticas. Una forma común de evaluarlo es informar la proporción de estadísticas que son (en cierto sentido) "tan o más extremas" que las estadísticas reales. Esto a menudo se llama un "valor p".

Debido a que el conjunto de datos real es una de esas permutaciones, su estadística necesariamente estará entre las que se encuentran dentro de la distribución de permutación. Por lo tanto, el valor p nunca puede ser cero.

A menos que el conjunto de datos sea muy pequeño (menos de aproximadamente 20-30 números totales, por lo general) o que el estadístico de prueba tenga una forma matemática particularmente agradable, no es factible generar todas las permutaciones. (Un ejemplo donde todas las permutaciones se generan aparece en prueba de permutación en R .) Por lo tanto implementaciones de ordenador de pruebas de permutación típicamente de muestra de la distribución de permutación. Lo hacen generando algunas permutaciones aleatorias independientes y esperan que los resultados sean una muestra representativa de todas las permutaciones.

Por lo tanto, cualquier número (como un "valor p") derivado de dicha muestra es solo estimadores de las propiedades de la distribución de permutación. Es muy posible, y a menudo ocurre cuando los efectos son grandes, que el valor p estimado sea ​​cero. No hay nada de malo en eso, pero inmediatamente plantea el problema hasta ahora descuidado de cuánto podría diferir el valor p estimado del correcto. Debido a que la distribución muestral de una proporción (como un valor p estimado) es binomial, esta incertidumbre se puede abordar con un intervalo de confianza binomial .

Arquitectura

Una implementación bien construida seguirá de cerca la discusión en todos los aspectos. Comenzaría con una rutina para calcular la estadística de prueba, como esta para comparar las medias de dos grupos:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

Escriba otra rutina para generar una permutación aleatoria del conjunto de datos y aplique la estadística de prueba. La interfaz de este permite al llamante proporcionar la estadística de prueba como argumento. Comparará los primeros melementos de una matriz (presuntamente un grupo de referencia) con los elementos restantes (el grupo "tratamiento").

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

La prueba de permutación se lleva a cabo primero encontrando la estadística de los datos reales (se supone que se almacenan en dos matrices) control y treatment) y luego buscando estadísticas para muchas permutaciones aleatorias independientes de los mismos:

z <- stat(control, treatment) # Test statistic for the observed data
sim<- sapply(1:1e4, f, sample=c(control,treatment), m=length(control), statistic=diff.means)

Ahora calcule la estimación binomial del valor p y un intervalo de confianza para él. Un método utiliza el binconfprocedimiento incorporado en el HMiscpaquete:

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

No es una mala idea comparar el resultado con otra prueba, incluso si se sabe que no es del todo aplicable: al menos podría obtener un sentido de orden de magnitud de dónde debería estar el resultado. En este ejemplo (de medias de comparación), una prueba t de Student generalmente da un buen resultado de todos modos:

t.test(treatment, control)

Esta arquitectura se ilustra en una situación más compleja, con Rcódigo de trabajo , en Probar si las variables siguen la misma distribución .

Ejemplo

$10$$0$$20$$1.5$ .

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

Después de usar el código anterior para ejecutar una prueba de permutación, tracé la muestra de la distribución de permutación junto con una línea roja vertical para marcar la estadística real:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

El cálculo del límite de confianza binomial resultó en

 PointEst Lower        Upper
        0     0 0.0003688199

$0$$0.00037$3.16e-05$0.00037$$0.00037$$0.05$$0.01$$0.001$ ).

Comentarios

$k$$N$ $k/N$$(k+1)/(N+1)$$N$

$10$$10^2=100$$0.000005$$1.6$$11.7$partes por millón: un poco más pequeño que la prueba t de Student informada. Aunque los datos se generaron con generadores de números aleatorios normales, lo que justificaría el uso de la prueba t de Student, los resultados de la prueba de permutación difieren de los resultados de la prueba t de Student porque las distribuciones dentro de cada grupo de observaciones no son perfectamente normales.

$\frac{B}{M}$$\frac{B+1}{M+1}$ ) es un estimador de valor p válido (pero conservador), no conduce a un rechazo excesivo del valor nulo.

(B es el número de permutaciones aleatorias en las que se obtiene una estadística mayor o igual que la observada y M es el número total de permutaciones aleatorias muestreadas).

$\frac{B}{M}$