Destacando resultados significativos de comparaciones múltiples no paramétricas en diagramas de caja

9

Tengo parcelas de 13 grupos que muestro en una parcela. Los grupos tienen poblaciones desequilibradas y normalmente no están distribuidos. Quiero mostrar qué pares son estadísticamente similares (es decir, tienen kruskal.test p-value <0.05) colocando a, b, c, etc. en la parte superior de los cuadros que coinciden. Aquí hay un pseudocódigo para mostrar lo que tengo:

A = c(1, 5, 8, 17, 16, 3, 24, 19, 6) 
B = c(2, 16, 5, 7, 4, 7, 3) 
C = c(1, 1, 3, 7, 9, 6, 10, 13) 
D = c(2, 15, 2, 9, 7) 
junk = list(g1=A, g2=B, g3=C, g4=D) 
boxplot(junk)

Aquí hay una trama que encontré que hace lo que quiero (excepto que tengo 13 grupos en una fila):

r statistical-significance boxplot Ilik
fuente

Respuestas:

7

El código más simple que viene a mi mente se muestra a continuación. Estoy bastante seguro de que hay algunas funciones ya existentes para hacer eso en CRAN, pero soy demasiado vago para buscarlas, incluso en R-seek .

dd <- data.frame(y=as.vector(unlist(junk)), 
                 g=rep(paste("g", 1:4, sep=""), unlist(lapply(junk, length))))

aov.res <- kruskal.test(y ~ g, data=dd)
alpha.level <- .05/nlevels(dd$g)  # Bonferroni correction, but use 
                                  # whatever you want using p.adjust()

# generate all pairwise comparisons
idx <- combn(nlevels(dd$g), 2)

# compute p-values from Wilcoxon test for all comparisons
pval.res <- numeric(ncol(idx))
for (i in 1:ncol(idx))
  # test all group, pairwise
  pval.res[i] <- with(dd, wilcox.test(y[as.numeric(g)==idx[1,i]], 
                                      y[as.numeric(g)==idx[2,i]]))$p.value

# which groups are significantly different (arranged by column)
signif.pairs <- idx[,which(pval.res<alpha.level)]

boxplot(y ~ g, data=dd, ylim=c(min(dd$y)-1, max(dd$y)+1))
# use offset= to increment space between labels, thanks to vectorization
for (i in 1:ncol(signif.pairs))
    text(signif.pairs[,i], max(dd$y)+1, letters[i], pos=4, offset=i*.8-1)

Aquí hay un ejemplo de lo que produciría el código anterior (con diferencias significativas entre los cuatro grupos):

ingrese la descripción de la imagen aquí

En lugar de la prueba de Wilcoxon, uno podría confiar en el procedimiento implementado en la kruskalmc()función del paquete pgirmess (vea una descripción del procedimiento utilizado aquí ).

Además, asegúrese de consultar los consejos R de Rudolf Cardinal sobre R: gráficos básicos 2 (ver en particular, Otro gráfico de barras, con anotaciones ).

chl
fuente
Gracias chl. No soy estadístico, así que tendré que aprender más tu respuesta, pero parece que hace lo que necesito.
Ilik
Si hay algo que no está claro, no dude en preguntar. No hay muchas cosas estadísticas aquí: utilicé las pruebas de Wilcoxon para las comparaciones grupales por pares y corregí los valores p individuales para limitar el riesgo general de falsos positivos al 5%. El paquete npmc incluye funciones adicionales para manejar comparaciones múltiples no paramétricas, pero también existe el marco de monedas . El resto es puramente código R usando gráficos base R; también se puede hacer con celosía o ggplot2.
chl
Pido disculpas por mi ignorancia en las estadísticas. Así que probé su código y primero noté que calculaba el kruskal.test pero no usaba el resultado (aov.res). Ahora entiendo que wilcox.test es un caso especial para kruskal para dos muestras. Pero luego intenté cambiar los valores en mis grupos para hacerlos (intuitivamente) diferentes y ver qué surge. A = c (10, 50, 18, 17, 16, 31, 24, 19, 6) B = c (10, 50, 18, 17, 16, 30, 25, 18, 7) C = c (1, 1, 3, 7, 9, 6, 10, 13) D = c (200, 158, 206, 119, 77). ¿g1 y g2 ahora son significativamente diferentes? y g3 y g4 no lo son? No estoy seguro de entender por qué.
Ilik
@Ilik Algo salió mal en mi código (la with(dd,instrucción estaba en el lugar equivocado, lo que resultó en una prueba extraña). Gracias por atrapar eso! Sí, no utilizo los resultados de la prueba KW, pero siempre es una buena idea verificarlo primero, de lo contrario, múltiples pruebas post-hoc no tendrían sentido (o al menos, sugerirían espionaje de datos). Tenga en cuenta que he corregido el código y nada resultó ser significativo, pero dejé la imagen original en aras de la claridad con resultados significativos.
chl
Para cualquier persona interesada, he cambiado un poco el código para escribir los resultados en el diagrama de caja: MyText = rep ('', nlevels (ddg))for(iin1:ncol(signif.pairs))MyText[signif.pairs[,i]]=paste(MyText[signif.pairs[,i]],letters[i],sep=″)text(c(1:nlevels(ddg)), rep (max (dd g)), MyText)y)+1,nlevels(dd
Ilik