Cómo realizar un ANOVA y pruebas post hoc utilizando R
Un tutorial detallado sobre el análisis ANOVA en R: supuestos, pruebas, funciones básicas de R, pruebas post hoc y paquete rstatix.
Imagen creada por mí con el paquete {aRtsy}Esta semana hice un tutorial sobre cómo realizar un análisis ANOVA utilizando R para el Psych Rstats Club, una plataforma de aprendizaje de estadísticas y R basada en ciencia abierta y colaborativa. Aquí está el guion, ¡que lo disfruten!
ANOVA
ANOVA, que significa “análisis de varianza” en inglés, es una extensión de la prueba t independiente de dos muestras. Utilizamos ANOVA cuando queremos determinar si las medias de más de dos grupos son significativamente diferentes. En este caso, tenemos varios grupos definidos por una variable de agrupación única, que a veces se denomina variable de factor.
Para esta demostración, utilizaremos el dataset iris. Como siempre, un buen primer paso es inspeccionar los datos:
# inspeccionar los datos
iris## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## 11 5.4 3.7 1.5 0.2 setosa
## 12 4.8 3.4 1.6 0.2 setosa
## 13 4.8 3.0 1.4 0.1 setosa
## 14 4.3 3.0 1.1 0.1 setosa
## 15 5.8 4.0 1.2 0.2 setosa
## 16 5.7 4.4 1.5 0.4 setosa
## 17 5.4 3.9 1.3 0.4 setosa
## 18 5.1 3.5 1.4 0.3 setosa
## 19 5.7 3.8 1.7 0.3 setosa
## 20 5.1 3.8 1.5 0.3 setosa
## 21 5.4 3.4 1.7 0.2 setosa
## 22 5.1 3.7 1.5 0.4 setosa
## 23 4.6 3.6 1.0 0.2 setosa
## 24 5.1 3.3 1.7 0.5 setosa
## 25 4.8 3.4 1.9 0.2 setosa
## 26 5.0 3.0 1.6 0.2 setosa
## 27 5.0 3.4 1.6 0.4 setosa
## 28 5.2 3.5 1.5 0.2 setosa
## 29 5.2 3.4 1.4 0.2 setosa
## 30 4.7 3.2 1.6 0.2 setosa
## 31 4.8 3.1 1.6 0.2 setosa
## 32 5.4 3.4 1.5 0.4 setosa
## 33 5.2 4.1 1.5 0.1 setosa
## 34 5.5 4.2 1.4 0.2 setosa
## 35 4.9 3.1 1.5 0.2 setosa
## 36 5.0 3.2 1.2 0.2 setosa
## 37 5.5 3.5 1.3 0.2 setosa
## 38 4.9 3.6 1.4 0.1 setosa
## 39 4.4 3.0 1.3 0.2 setosa
## 40 5.1 3.4 1.5 0.2 setosa
## 41 5.0 3.5 1.3 0.3 setosa
## 42 4.5 2.3 1.3 0.3 setosa
## 43 4.4 3.2 1.3 0.2 setosa
## 44 5.0 3.5 1.6 0.6 setosa
## 45 5.1 3.8 1.9 0.4 setosa
## 46 4.8 3.0 1.4 0.3 setosa
## 47 5.1 3.8 1.6 0.2 setosa
## 48 4.6 3.2 1.4 0.2 setosa
## 49 5.3 3.7 1.5 0.2 setosa
## 50 5.0 3.3 1.4 0.2 setosa
## 51 7.0 3.2 4.7 1.4 versicolor
## 52 6.4 3.2 4.5 1.5 versicolor
## 53 6.9 3.1 4.9 1.5 versicolor
## 54 5.5 2.3 4.0 1.3 versicolor
## 55 6.5 2.8 4.6 1.5 versicolor
## 56 5.7 2.8 4.5 1.3 versicolor
## 57 6.3 3.3 4.7 1.6 versicolor
## 58 4.9 2.4 3.3 1.0 versicolor
## 59 6.6 2.9 4.6 1.3 versicolor
## 60 5.2 2.7 3.9 1.4 versicolor
## 61 5.0 2.0 3.5 1.0 versicolor
## 62 5.9 3.0 4.2 1.5 versicolor
## 63 6.0 2.2 4.0 1.0 versicolor
## 64 6.1 2.9 4.7 1.4 versicolor
## 65 5.6 2.9 3.6 1.3 versicolor
## 66 6.7 3.1 4.4 1.4 versicolor
## 67 5.6 3.0 4.5 1.5 versicolor
## 68 5.8 2.7 4.1 1.0 versicolor
## 69 6.2 2.2 4.5 1.5 versicolor
## 70 5.6 2.5 3.9 1.1 versicolor
## 71 5.9 3.2 4.8 1.8 versicolor
## 72 6.1 2.8 4.0 1.3 versicolor
## 73 6.3 2.5 4.9 1.5 versicolor
## 74 6.1 2.8 4.7 1.2 versicolor
## 75 6.4 2.9 4.3 1.3 versicolor
## 76 6.6 3.0 4.4 1.4 versicolor
## 77 6.8 2.8 4.8 1.4 versicolor
## 78 6.7 3.0 5.0 1.7 versicolor
## 79 6.0 2.9 4.5 1.5 versicolor
## 80 5.7 2.6 3.5 1.0 versicolor
## 81 5.5 2.4 3.8 1.1 versicolor
## 82 5.5 2.4 3.7 1.0 versicolor
## 83 5.8 2.7 3.9 1.2 versicolor
## 84 6.0 2.7 5.1 1.6 versicolor
## 85 5.4 3.0 4.5 1.5 versicolor
## 86 6.0 3.4 4.5 1.6 versicolor
## 87 6.7 3.1 4.7 1.5 versicolor
## 88 6.3 2.3 4.4 1.3 versicolor
## 89 5.6 3.0 4.1 1.3 versicolor
## 90 5.5 2.5 4.0 1.3 versicolor
## 91 5.5 2.6 4.4 1.2 versicolor
## 92 6.1 3.0 4.6 1.4 versicolor
## 93 5.8 2.6 4.0 1.2 versicolor
## 94 5.0 2.3 3.3 1.0 versicolor
## 95 5.6 2.7 4.2 1.3 versicolor
## 96 5.7 3.0 4.2 1.2 versicolor
## 97 5.7 2.9 4.2 1.3 versicolor
## 98 6.2 2.9 4.3 1.3 versicolor
## 99 5.1 2.5 3.0 1.1 versicolor
## 100 5.7 2.8 4.1 1.3 versicolor
## 101 6.3 3.3 6.0 2.5 virginica
## 102 5.8 2.7 5.1 1.9 virginica
## 103 7.1 3.0 5.9 2.1 virginica
## 104 6.3 2.9 5.6 1.8 virginica
## 105 6.5 3.0 5.8 2.2 virginica
## 106 7.6 3.0 6.6 2.1 virginica
## 107 4.9 2.5 4.5 1.7 virginica
## 108 7.3 2.9 6.3 1.8 virginica
## 109 6.7 2.5 5.8 1.8 virginica
## 110 7.2 3.6 6.1 2.5 virginica
## 111 6.5 3.2 5.1 2.0 virginica
## 112 6.4 2.7 5.3 1.9 virginica
## 113 6.8 3.0 5.5 2.1 virginica
## 114 5.7 2.5 5.0 2.0 virginica
## 115 5.8 2.8 5.1 2.4 virginica
## 116 6.4 3.2 5.3 2.3 virginica
## 117 6.5 3.0 5.5 1.8 virginica
## 118 7.7 3.8 6.7 2.2 virginica
## 119 7.7 2.6 6.9 2.3 virginica
## 120 6.0 2.2 5.0 1.5 virginica
## 121 6.9 3.2 5.7 2.3 virginica
## 122 5.6 2.8 4.9 2.0 virginica
## 123 7.7 2.8 6.7 2.0 virginica
## 124 6.3 2.7 4.9 1.8 virginica
## 125 6.7 3.3 5.7 2.1 virginica
## 126 7.2 3.2 6.0 1.8 virginica
## 127 6.2 2.8 4.8 1.8 virginica
## 128 6.1 3.0 4.9 1.8 virginica
## 129 6.4 2.8 5.6 2.1 virginica
## 130 7.2 3.0 5.8 1.6 virginica
## 131 7.4 2.8 6.1 1.9 virginica
## 132 7.9 3.8 6.4 2.0 virginica
## 133 6.4 2.8 5.6 2.2 virginica
## 134 6.3 2.8 5.1 1.5 virginica
## 135 6.1 2.6 5.6 1.4 virginica
## 136 7.7 3.0 6.1 2.3 virginica
## 137 6.3 3.4 5.6 2.4 virginica
## 138 6.4 3.1 5.5 1.8 virginica
## 139 6.0 3.0 4.8 1.8 virginica
## 140 6.9 3.1 5.4 2.1 virginica
## 141 6.7 3.1 5.6 2.4 virginica
## 142 6.9 3.1 5.1 2.3 virginica
## 143 5.8 2.7 5.1 1.9 virginica
## 144 6.8 3.2 5.9 2.3 virginica
## 145 6.7 3.3 5.7 2.5 virginica
## 146 6.7 3.0 5.2 2.3 virginica
## 147 6.3 2.5 5.0 1.9 virginica
## 148 6.5 3.0 5.2 2.0 virginica
## 149 6.2 3.4 5.4 2.3 virginica
## 150 5.9 3.0 5.1 1.8 virginicaRecordá que podés ejecutar una línea de código haciendo clic en “Run” cuando estés en esa línea, o podés usar el atajo de teclado control+enter o command+enter en Mac.
En este conjunto de datos, tenemos la longitud y el ancho del sépalo y los pétalos de las flores de cada especie. En este ejercicio, nos gustaría saber si las diferentes especies difieren en cuanto al ancho de sus sépalos.
Antes de adentrarnos en los aspectos estadísticos de ANOVA, discutamos brevemente sus supuestos. ANOVA asume que los datos están distribuidos de forma normal, las varianzas son desconocidas pero iguales entre los grupos, y las muestras son independientes.
Para verificar el supuesto de normalidad, podemos utilizar la prueba de Shapiro-Wilk. Tenemos que ejecutar la función shapiro.test() para la variable dependiente (el ancho del sépalo) de cada especie, por lo que usamos una sintaxis específica para seleccionar los datos del conjunto de datos.
# prueba de Shapiro-Wilk (supuesto de normalidad)
shapiro.test(iris$Sepal.Width[iris$Species == "setosa"])##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Width[iris$Species == "setosa"]
## W = 0.97172, p-value = 0.2715# primero obtenemos los valores del ancho del sépalo y luego seleccionamos esos valores para la especie setosa
# luego repetimos esto para cada especie, copiando y pegando
shapiro.test(iris$Sepal.Width[iris$Species == "virginica"])##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Width[iris$Species == "virginica"]
## W = 0.96739, p-value = 0.1809shapiro.test(iris$Sepal.Width[iris$Species == "versicolor"])##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Width[iris$Species == "versicolor"]
## W = 0.97413, p-value = 0.338Según los resultados de estas pruebas, ninguno de ellos arrojó resultados estadísticamente significativos, lo que indica que los datos están distribuidos de forma normal.
Para evaluar el supuesto de homogeneidad de varianzas, podemos utilizar la prueba de Levene del paquete car. El primer argumento es la fórmula: recuerda que siempre usamos primero la variable dependiente, luego la tilde o virgulilla, y por último la variable de agrupación. El segundo argumento es el conjunto de datos.
# Prueba de Levene (supuesto de homoscedasticidad)
car::leveneTest(Sepal.Width ~ Species, data = iris)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5902 0.5555
## 147La prueba de Levene evalúa si la varianza del ancho del sépalo es igual en diferentes grupos de especies. En este caso, dado que el valor p no es estadísticamente significativo, podemos concluir que se cumple el supuesto de homogeneidad de varianzas.
Ahora que hemos verificado los supuestos, podemos proceder con el análisis ANOVA. Sin embargo, es importante tener en cuenta que si tus datos no siguen una distribución normal, podés considerar el uso de pruebas no paramétricas como las pruebas de Kruskal-Wallis o Friedman. Además, si se viola el supuesto de homogeneidad de varianzas, se puede utilizar un ANOVA de Welch.
Realicemos el análisis ANOVA utilizando funciones básicas de R:
# ANOVA con funciones básicas de R
iris_anova <- aov(Sepal.Width ~ Species, data = iris)
# para ver el resultado
summary(iris_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.35 5.672 49.16 <2e-16 ***
## Residuals 147 16.96 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El resumen de la salida ANOVA muestra que el análisis ANOVA es estadísticamente significativo. Sin embargo, no proporciona información sobre qué grupos específicos difieren entre sí. Para determinar las diferencias por pares, podemos utilizar pruebas post hoc como la prueba de Diferencia Significativa Honesta (HSD) de Tukey, también conocida como corrección de Tukey:
# Corrección de Tukey
TukeyHSD(iris_anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802#tukeyhsd(iris_anova) # error!Recordá que R distingue entre mayúsculas y minúsculas, así que es importante prestar atención a la capitalización.
La prueba de Tukey revela que los tres grupos son estadísticamente diferentes entre sí. Alternativamente, también podés utilizar otras pruebas post hoc, como la corrección de Bonferroni o el método de Holm, para ajustar los valores p para comparaciones múltiples.
# Corrección de Bonferroni
# utilizamos una función que calcula una prueba t y especificamos que queremos una corrección de Bonferroni
# el primer argumento es la variable dependiente y el segundo es la variable de agrupación o factor.
pairwise.t.test(iris$Sepal.Width, iris$Species, p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: iris$Sepal.Width and iris$Species
##
## setosa versicolor
## versicolor < 2e-16 -
## virginica 1.4e-09 0.0094
##
## P value adjustment method: bonferroniTambién podés utilizar el método de Holm, con la misma sintaxis:
# Corrección de Holm
pairwise.t.test(iris$Sepal.Width, iris$Species, p.adjust.method = "holm")##
## Pairwise comparisons using t tests with pooled SD
##
## data: iris$Sepal.Width and iris$Species
##
## setosa versicolor
## versicolor < 2e-16 -
## virginica 9.1e-10 0.0031
##
## P value adjustment method: holmAunque los valores p pueden variar ligeramente entre estas correcciones, la conclusión general sigue siendo la misma.
El paquete rstatix
Ahora, exploremos un paquete llamado rstatix, que simplifica el análisis ANOVA y también es compatible con las funciones de tidyverse:
# para instalar el paquete
install.packages("rstatix")
# para cargar el paquete
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterCon rstatix, podemos realizar el análisis ANOVA:
# ANOVA con rstatix
anova_test(formula = Sepal.Width ~ Species, data = iris)## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 Species 2 147 49.16 4.49e-17 * 0.401Los resultados resumidos proporcionados por la función anova_test incluyen no solo el valor p, sino también el tamaño del efecto. De forma predeterminada, la función anova_test calcula un eta cuadrad generalizado, pero también podés calcular un eta cuadrado parcial si es necesario:
# Eta cuadrado parcial
anova_test(data = iris, formula = Sepal.Width ~ Species, effect.size = "pes")## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 pes
## 1 Species 2 147 49.16 4.49e-17 * 0.401# en este caso es lo mismorstatix también proporciona su propia función para realizar la corrección de Tukey:
# Corrección de Tukey con rstatix
tukey_hsd(x = iris, formula = Sepal.Width ~ Species)## # A tibble: 3 × 9
## term group1 group2 null.value estimate conf.low conf.high p.adj
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Species setosa versicolor 0 -0.658 -0.819 -0.497 3.10e-14
## 2 Species setosa virginica 0 -0.454 -0.615 -0.293 1.36e- 9
## 3 Species versicolor virginica 0 0.204 0.0431 0.365 8.78e- 3
## # ℹ 1 more variable: p.adj.signif <chr>En general, estas son algunas de las formas en que puedes realizar el análisis ANOVA en R. Espero que esta explicación te ayude a aclarar el proceso. Si tienes alguna pregunta o necesitas más aclaraciones, ¡escribime!
Como siempre, recordá que podés suscribirte a mi blog para no perderte ninguna actualización, y si tenés alguna pregunta, no dudes en contactarme. Y si te gusta lo que hago, puedes invitarme a tomar un cafecito desde Argentina o un kofi desde otros países.
Macarena Quiroga
Lingüista/Becaria doctoral
Investigo la adquisición del lenguaje. Estudio estadística y ciencia de datos con R/Rstudio. Si te gusta lo que hago, podés invitarme un cafecito desde Argentina, o un kofi desde otros países. Suscribite a mi blog aquí.