(Nota del autor: El objetivo de esta píldora es la de mostrar boxplots con los p-valores de comparaciones post-hoc, por eso no se detalla toda la teoría que hay detrás del ANOVA y de los post-hoc, ni se detalla si los contrastes son paramétricos o no paramétricos, ni se habla de corrección de p-valores por tests múltiples, etc)

Cuando queremos estudiar el efecto que producen las categorías de un factor \(X\) en una variable numérica \(Y\) realizamos un Análisis de la Varianza (ANOVA) de un factor (+info: (https://en.wikipedia.org/wiki/Analysis_of_variance)). Es un contraste que se basa en la distribución \(F\) de Fisher-Snedecor cuya hipótesis nula es \(H_0\): El factor X no influye en la variable Y. El resultado suele proporcionarse en una tabla con información sobre la variabilidad y el p-valor de la \(F\).

En el siguiente ejemplo vamos a estudiar la influencia del factor \(X\) = “tipo de medicamento para el colesterol” (con categorías {‘M1’, ‘M2’, ‘M3’, ‘M4’}) sobre la variable \(Y\) = “variación del nivel de colesterol”. Simulamos unos datos, observamos la media de \(Y\) en cada tipo de medicamento y aplicamos el test ANOVA aov() para ver si hay diferencias significativas entre dichas medias.

library(dplyr)    # install.packages('dplyr') si no está instalado
set.seed(1968)                                 
data <- data.frame(X = factor(rep(c('M1', 'M2', 'M3', 'M4'), each = 100)),
                   Y = c(rnorm(100), rnorm(100, 0.75, 2.5), rnorm(100),
                         rnorm(100, -0.5)))
data %>% group_by(X) %>% summarise(mean_Y = mean(Y))                               
# A tibble: 4 × 2
  X      mean_Y
  <fct>   <dbl>
1 M1     0.0202
2 M2     0.765 
3 M3     0.0369
4 M4    -0.495 
summary(aov(Y ~ X, data = data))
             Df Sum Sq Mean Sq F value   Pr(>F)    
X             3   80.5  26.826   11.83 1.94e-07 ***
Residuals   396  897.6   2.267                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Observamos en la tabla ANOVA un p-valor = 1.9399169^{-7} < 0.05 por lo que tenemos evidencias suficientes para rechazar la hipótesis nula \(H_0\), y por tanto concluir que el tipo de medicamento \(X\) influye en la variación del colesterol \(Y\).

El test ANOVA nos dice únicamente si la variable \(Y\) presenta (o no) diferencias significativas entre la categorías de \(X\), pero no nos dice entre qué categorías son esas diferencias. Para aclarar este punto se realizan los llamados contrastes ‘post-hoc’ que son comparaciones 2 a 2 de las categorías del factor \(X\) (+info: (https://en.wikipedia.org/wiki/Post_hoc_analysis)).

En nuestro ejemplo hay evidencias de que el tipo de medicamento influye en el colesterol así que realizamos los contrastes post-hoc con la función pairwise.t.test() que realiza todos los contrastes simultáneamente:

pairwise.t.test(data$Y, data$X, p.adjust.method = 'none')

    Pairwise comparisons using t tests with pooled SD 

data:  data$Y and data$X 

   M1      M2      M3     
M2 0.00052 -       -      
M3 0.93747 0.00069 -      
M4 0.01603 7.1e-09 0.01293

P value adjustment method: none

Observamos que para una significación de \(\alpha\) = 0.05 los únicos medicamentos para los que el colesterol no presenta diferencias significativas son “M1” vs “M3”.

Ahora llegamos al objetivo de esta píldora estadística, que es representar la información de los contrastes post-hoc en un gráfico de cajas. Cada caja muestra la distribución de la variable \(Y\) en las categorías de la variable \(X\). Mediante segmentos podemos unir estas cajas de 2 en 2 e indicar el p-valor del correspondiente contraste post-hoc.

Utilizamos la librería ggplot2 para construir el gráfico boxplot y la librería ggsignif (https://github.com/const-ae/ggsignif) para añadir la información de los contrastes al gráfico, usando la geometría geom_signif().

library(ggplot2)      # install.packages('ggplot2') si no está instalada
library(ggsignif)     # install.packages('ggsignif') si no está instalada
# Boxplot para visualizar cambios entre categorías
g <- ggplot(data = data, aes(x = X, y = Y)) + geom_boxplot()
g + ggtitle('Boxplot de Y en cada categoría de X')

# Calculamos todas las combinaciones de 2 en 2 de las categorías con combn()
comparisons <- as.data.frame(combn(levels(data$X), 2))
g + geom_signif(comparisons = comparisons, test = 't.test', vjust = 0.3,
                y_position = seq(max(data$Y)+0.5, length.out = ncol(comparisons)),
                textsize = 4, size = 0.5, col = 'darkblue') +
  ggtitle('Pairwise t-test p-values')

(Nota: los p-valores del gráfico no coinciden exactamente con los del post-hoc debido al tipo de corrección por test múltiples. Se muestran los p-valores de un T test pero se podrían mostrar los de un Wilcoxon test)

En lugar de añadir a cada segmento el p-valor se puede utilizar una notación con asteriscos que significan: \(***\) = p-valor ≤ 0.001, \(**\) = p-valor ≤ 0.01, \(*\) = p-valor ≤ 0.05 y \(NS\) = p-valor > 0.05. Para ello se añade el parámetro map_signif_level = TRUE.

# Todas las combinaciones de 2 en 2 con notación de asteriscos
g + geom_signif(comparisons = comparisons, test = 't.test', vjust = 0.2,
                y_position = seq(max(data$Y)+0.5, length.out = ncol(comparisons)),
                textsize = 6, size = 0.5, col = 'darkred', map_signif_level = TRUE) +
  ggtitle('Pairwise t-test significance')