¿Qué implementación de la prueba de permutación en R debe utilizarse en lugar de las pruebas t (emparejadas y no emparejadas)?

Question

Tengo datos de un experimento que he analizado mediante pruebas t. La variable dependiente es de escala de intervalo y los datos son no apareados (es decir, 2 grupos) o apareados (es decir, dentro de los sujetos). Por ejemplo, (dentro de los sujetos):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

Sin embargo, los datos no son normales, por lo que un revisor nos pidió que utilizáramos algo distinto a la prueba t. Sin embargo, como se puede ver fácilmente, los datos no sólo no están distribuidos normalmente, sino que las distribuciones no son iguales entre las condiciones:

Por lo tanto, no se pueden utilizar las pruebas no paramétricas habituales, la prueba U de Mann-Whitney (no emparejada) y la prueba de Wilcoxon (emparejada), ya que requieren distribuciones iguales entre las condiciones. Por lo tanto, decidí que lo mejor sería un remuestreo o una prueba de permutación.

Ahora, estoy buscando una implementación en R de un equivalente basado en permutaciones de la prueba t, o cualquier otro consejo sobre qué hacer con los datos.

Sé que hay algunos paquetes de R que pueden hacer esto por mí (por ejemplo, coin, perm, exactRankTest, etc.), pero no sé cuál elegir. Así que, si alguien con experiencia en el uso de estas pruebas puede darme un empujón, sería genial.

ACTUALIZACIÓN: Sería ideal si pudiera proporcionar un ejemplo de cómo informar de los resultados de esta prueba.

Answer 1

No debería importar mucho, ya que la estadística de la prueba siempre será la diferencia de medias (o algo equivalente). Las pequeñas diferencias pueden provenir de la implementación de los métodos de Monte-Carlo. Pruebe los tres paquetes con sus datos con una prueba unilateral para dos variables independientes:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

Para comprobar el valor p exacto con un cálculo manual de todas las permutaciones, restringiré los datos a los 9 primeros valores.

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coin y exactRankTests son ambos del mismo autor, pero coin parece ser más general y extensa, también en cuanto a la documentación. exactRankTests ya no se desarrolla activamente. Por lo tanto, yo elegiría coin (también por funciones informativas como support() ), a menos que no le guste tratar con objetos S4.

EDIT: para dos variables dependientes, la sintaxis es

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081

Answer 2

Creo que es necesario hacer algunos comentarios.

1) Te animo a que pruebes con múltiples representaciones visuales de tus datos, porque pueden captar cosas que se pierden con (gráficos como) los histogramas, y también te recomiendo encarecidamente que los representes en ejes paralelos. En este caso, no creo que los histogramas hagan un buen trabajo para comunicar las características más destacadas de sus datos. Por ejemplo, eche un vistazo a los boxplots de lado a lado:

boxplot(x1, y1, names = c("x1", "y1"))

O incluso gráficos de bandas de lado a lado:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

¡Mira los centros, las extensiones y las formas de estos! Alrededor de tres cuartas partes de los $x1$ datos se sitúan muy por encima de la mediana del $y1$ datos. La difusión de $x1$ es minúsculo, mientras que la propagación de $y1$ es enorme. Ambos $x1$ y $y1$ están muy sesgados a la izquierda, pero de forma diferente. Por ejemplo, $y1$ tiene cinco (!) valores repetidos de cero.

2) No has explicado con mucho detalle de dónde proceden tus datos, ni cómo se han medido, pero esta información es muy importante a la hora de seleccionar un procedimiento estadístico. ¿Son sus dos muestras anteriores independientes? ¿Hay alguna razón para creer que las distribuciones marginales de las dos muestras deberían ser las mismas (excepto por una diferencia de ubicación, por ejemplo)? ¿Cuáles fueron las consideraciones antes al estudio que le llevó a buscar pruebas de una diferencia entre los dos grupos?

3) La prueba t no es apropiada para estos datos porque las distribuciones marginales son marcadamente no normales, con valores extremos en ambas muestras. Si lo desea, podría apelar al CLT (debido a su muestra de tamaño moderado) para utilizar una $z$ -(que sería similar a una prueba z para muestras grandes), pero dada la asimetría (en ambas variables) de sus datos no juzgaría muy convincente tal recurso. Claro, de todos modos se puede utilizar para calcular un $p$ -valor, pero ¿qué hace eso por ti? Si los supuestos no se cumplen, entonces un $p$ -El valor es sólo una estadística; no dice lo que usted (presumiblemente) quiere saber: si hay evidencia de que las dos muestras provienen de distribuciones diferentes.

4) Una prueba de permutación tampoco sería apropiada para estos datos. La única y a menudo olvidada suposición para las pruebas de permutación es que las dos muestras son intercambiable bajo la hipótesis nula. Eso significaría que tienen distribuciones marginales idénticas (bajo la nula). Pero estás en problemas, porque los gráficos sugieren que las distribuciones difieren tanto en ubicación como en escala (y también en forma). Por lo tanto, no se puede probar (válidamente) una diferencia en la ubicación porque las escalas son diferentes, y no se puede probar (válidamente) una diferencia en la escala porque las ubicaciones son diferentes. Uy. De nuevo, puedes hacer la prueba de todos modos y obtener un $p$ -valor, pero ¿y qué? ¿Qué has conseguido realmente?

5) En mi opinión, estos datos son un ejemplo perfecto (?) de que una imagen bien elegida vale más que 1000 pruebas de hipótesis. No necesitamos la estadística para distinguir entre un lápiz y un granero. La afirmación adecuada, en mi opinión, para estos datos sería "Estos datos muestran marcadas diferencias con respecto a la ubicación, la escala y la forma." Podrías seguir con estadísticas descriptivas (robustas) para cada una de ellas para cuantificar las diferencias, y explicar qué significan las diferencias en el contexto de tu estudio original.

6) Su revisor probablemente (y lamentablemente) va a insistir en algún tipo de $p$ -como condición previa a la publicación. Suspiro. Si fuera yo, dadas las diferencias con respecto a todo Probablemente utilizaría una prueba no paramétrica de Kolmogorov-Smirnov para escupir un $p$ -valor que demuestre que las distribuciones son diferentes, y luego proceder con las estadísticas descriptivas como arriba. Tendrías que añadir algo de ruido a las dos muestras para eliminar los empates. (Y, por supuesto, todo esto supone que tus muestras son independientes, lo cual no has declarado explícitamente).

Esta respuesta es mucho más larga de lo que pretendía en un principio. Lo siento.

Answer 3

Como esta pregunta ha vuelto a surgir, puedo añadir otra respuesta inspirada en una reciente entrada del blog vía R-Bloggers de Robert Kabacoff, el autor de Quick-R y R en acción utilizando el lmPerm paquete.

Sin embargo, este método produce resultados muy contrastados (y muy inestables) con el producido por el coin en la respuesta de @caracakl (el valor p del análisis intra-sujeto es 0.008 ). El análisis toma la preparación de datos de la respuesta de @caracal también:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

produce:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839

Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Si se ejecuta esto varias veces, los valores p saltan entre ~.05 y ~.1.

Aunque se trata de una respuesta a la pregunta, permítame plantear una pregunta al final (puedo trasladar esto a una nueva pregunta si lo desea):
¿Alguna idea de por qué este análisis es tan inestable y produce unos valores p tan divergentes a los del análisis de la moneda? ¿He hecho algo mal?

Answer 4

Mis comentarios no se refieren a la aplicación de la prueba de permutación, sino a las cuestiones más generales que plantean estos datos y el debate al respecto, en particular el post de G. Jay Kerns.

En realidad, las dos distribuciones me parecen bastante similares, EXCEPTO el grupo de 0 en Y1, que es muy diferente de las demás observaciones de esa muestra (la siguiente más pequeña es de unos 50 en la escala 0-100), así como todas las de X1. Yo investigaría primero si hay algo diferente en esas observaciones.

En segundo lugar, suponiendo que esos 0s pertenezcan al análisis, decir que la prueba de permutación no es válida porque las distribuciones parecen ser diferentes es una pregunta. Si la nula fuera cierta (las distribuciones son idénticas), ¿podrías (con una probabilidad razonable) obtener distribuciones que parezcan tan diferentes como estas dos? Responder a eso es el objetivo de la prueba, ¿no? Tal vez en este caso algunos consideren que la respuesta es obvia sin hacer la prueba, pero con estas distribuciones tan pequeñas y peculiares, no creo que yo lo haga.

Answer 5

¿Son estas puntuaciones proporciones? Si es así, ciertamente no deberías utilizar una prueba paramétrica gaussiana, y aunque podrías seguir adelante con un enfoque no paramétrico como una prueba de permutación o un bootstrap de las medias, sugeriría que obtendrías más poder estadístico empleando un enfoque paramétrico no gaussiano adecuado. Específicamente, cada vez que pueda calcular una medida de proporción dentro de una unidad de interés (por ejemplo, un participante en un experimento), puede y probablemente debería utilizar un modelo de efectos mixtos que especifique las observaciones con un error distribuido binomialmente. Véase Dixon 2004 .