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

Tengo datos de un experimento que analicé usando pruebas t. La variable dependiente se escala a intervalos y los datos no están emparejados (es decir, 2 grupos) o están emparejados (es decir, dentro de los sujetos). Por ejemplo (dentro de las asignaturas):

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 usemos algo más que la prueba t. Sin embargo, como se puede ver fácilmente, los datos no solo no se distribuyen normalmente, sino que las distribuciones no son iguales entre las condiciones:

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

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

Sé que hay algunos paquetes R que pueden hacer esto por mí (por ejemplo, coin, perm, exactoRankTest, etc.), pero no sé cuál elegir. Por lo tanto, si alguien con algo de experiencia en el uso de estas pruebas me diera un impulso inicial, sería excelente.

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

No debería importar mucho ya que la estadística de prueba siempre será la diferencia en las medias (o algo equivalente). Pequeñas diferencias pueden venir de la implementación de los métodos de Monte-Carlo. Probar 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 verificar el valor p exacto con un cálculo manual de todas las permutaciones, restringiré los datos a los primeros 9 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

coiny exactRankTestsambos son del mismo autor, pero coinparece ser más general y extenso, también en términos de documentación. exactRankTestsya no se desarrolla activamente. Por lo tanto, elegiría coin(también debido a funciones informativas como support()), a menos que no le guste tratar con objetos S4.

EDITAR: 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

Algunos comentarios son, creo, en orden.

1) Le animo a que pruebe varias pantallas visuales de sus datos, ya que pueden capturar cosas que se pierden con los histogramas (como gráficos), y también le recomiendo encarecidamente que dibuje en ejes uno al lado del otro. En este caso, no creo que los histogramas hagan un muy buen trabajo al comunicar las características más destacadas de sus datos. Por ejemplo, eche un vistazo a los diagramas de caja de lado a lado:

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

O incluso diagramas de tira de lado a lado

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

¡Mira los centros, los spreads y las formas de estos! Alrededor de las tres cuartas partes de los datos están muy por encima de la mediana de los datos . La propagación de es pequeña, mientras que la propagación de es enorme. Tanto como están muy sesgadas a la izquierda, pero de diferentes maneras. Por ejemplo, tiene cinco (!) Valores repetidos de cero.y 1 x 1 y 1 x 1 y 1 y $x1$$y1$$x1$$y1$$x1$$y1$$y1$

2) No explicó con mucho detalle de dónde provienen sus datos, ni cómo se midieron, pero esta información es muy importante cuando llega el momento de seleccionar un procedimiento estadístico. ¿Sus dos muestras anteriores son 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 previas al estudio que lo llevaron a buscar evidencia 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, puede recurrir al CLT (debido a su muestra de tamaño moderado) para usar una prueba (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 una apelación tan convincente. Claro, puede usarlo de todos modos para calcular un valor , pero ¿qué hace eso por usted? Si los supuestos no se cumplen, entonces un valor es solo una estadística; no dice lo que (supuestamente) quiere saber: si hay evidencia de que las dos muestras provienen de diferentes distribuciones.p $z$$p$$p$

4) Una prueba de permutación también sería inapropiada para estos datos. La suposición única y a menudo pasada por alto para las pruebas de permutación es que las dos muestras son intercambiables bajo la hipótesis nula. Eso significaría que tienen distribuciones marginales idénticas (bajo nulo). Pero está 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 puede probar (válidamente) una diferencia de ubicación porque las escalas son diferentes, y no puede probar (válidamente) una diferencia de escala porque las ubicaciones son diferentes. Ups De nuevo, puedes hacer la prueba de todos modos y obtener un valor , pero ¿y qué? ¿Qué has logrado realmente?$p$

5) En mi opinión, estos datos son un ejemplo perfecto (?) De que una imagen bien elegida vale 1000 pruebas de hipótesis. No necesitamos estadísticas para diferenciar entre un lápiz y un granero. La declaración apropiada en mi opinión para estos datos sería "Estos datos exhiben marcadas diferencias con respecto a la ubicación, escala y forma". Puede hacer un seguimiento con estadísticas descriptivas (robustas) para cada una de ellas para cuantificar las diferencias y explicar qué significan las diferencias en el contexto de su estudio original.

6) Su revisor probablemente (y tristemente) va a insistir en algún tipo de valor como condición previa para la publicación. ¡Suspiro! Si fuera yo, dadas las diferencias con respecto a todo , probablemente usaría una prueba no paramétrica de Kolmogorov-Smirnov para escupir un valor que demuestre que las distribuciones son diferentes, y luego proceder con estadísticas descriptivas como arriba. Debería agregar algo de ruido a las dos muestras para deshacerse de los lazos. (Y, por supuesto, todo esto supone que sus muestras son independientes, lo que no indicó explícitamente).$p$$p$

Esta respuesta es mucho más larga de lo que originalmente pretendía que fuera. Lo siento por eso.

Mis comentarios no se refieren a la implementación de la prueba de permutación, sino a los problemas más generales planteados por estos datos y su discusión, en particular la publicación de G. Jay Kerns.

Las dos distribuciones en realidad se parecen bastante a mí, EXCEPTO para el grupo de 0s en Y1, que son muy diferentes de las otras observaciones en esa muestra (la más pequeña es aproximadamente 50 en la escala 0-100), así como todas las de X1. Primero investigaría si había algo diferente en esas observaciones.

En segundo lugar, suponiendo que esos 0s pertenecen al análisis, decir que la prueba de permutación no es válida porque las distribuciones parecen diferir plantea la pregunta. Si el valor nulo fuera verdadero (las distribuciones son idénticas), ¿podría (con una probabilidad razonable) hacer que las distribuciones se vean tan diferentes como estas dos? Respondiendo que ese es el objetivo de la prueba, ¿no? Quizás en este caso algunos consideren obvia la respuesta sin ejecutar la prueba, pero con estas distribuciones pequeñas y peculiares, no creo que lo haga.

Como esta pregunta apareció nuevamente, puedo agregar otra respuesta inspirada en una reciente publicación de blog a través de R-Bloggers de Robert Kabacoff, el autor de Quick-R y R in Action usando el lmPermpaquete.

Sin embargo, este método produce resultados muy contrastantes (y muy inestables) con el producido por el coinpaquete en la respuesta de @caracakl (el valor p del análisis dentro de los sujetos 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 ejecuta esto varias veces, los valores p saltan entre ~ .05 y ~ .1.

Aunque es una respuesta a la pregunta, permítame plantear una pregunta al final (puedo moverla a una nueva pregunta si lo desea):
cualquier idea de por qué este análisis es tan inestable y produce valores p tan divergentes para el analisis de monedas? ¿Hice algo mal?

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

Simplemente agregando otro enfoque, ezPermde ezpaquete:

> # preparing the data
> DV <- c(x1, y1)
> IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
> id <- factor(rep(1:length(x1), 2))
> df <- data.frame(id=id,DV=DV,IV=IV)
>
> library(ez)
> ezPerm( data = df, dv = DV, wid = id, within = IV, perms = 1000)
|=========================|100%              Completed after 17 s 
  Effect     p p<.05
1     IV 0.016     *

Esto parece ser consistente con el oneway_testdel coinpaquete:

> library(coin)
> pvalue(oneway_test(DV ~ IV | id,  distribution=approximate(B=999999)))
[1] 0.01608002
99 percent confidence interval:
 0.01575782 0.01640682

Sin embargo, tenga en cuenta que este no es el mismo ejemplo proporcionado por @caracal . En su ejemplo, que incluye alternative="greater", por lo tanto, la diferencia en los valores de p ~0.008vs ~0.016.

El aovppaquete sugirió en una de las respuestas producir recelo disminuir los valores de p, y se ejecuta sospechoso rápido, incluso cuando intento valores altos para el Iter, Cay maxIterargumentos:

library(lmPerm)
summary(aovp(DV ~ IV + Error(id/IV), data=df,  maxIter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Iter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Ca = 0.00000000001))

Dicho esto, los argumentos parecen estar reduciendo ligeramente las variaciones de los valores p de ~.03y ~.1(obtuve un rango mayor que el informado por @Henrik), a 0.03y 0.07.