ANOVA

Este material es parte de la Unidad 2 del Curso de Epidemiología - Nivel Avanzado del Instituto Nacional de Epidemiología “Dr. Juan H. Jara” - ANLIS v2024

ANOVA por Christian Ballejo bajo licencia CC BY-NC 4.0

Introducción

La técnica de análisis de varianza (ANOVA) es el test estadístico a emplear cuando se desea comparar las medias de más de dos grupos.

La hipótesis nula (\(H_0\)) de la prueba es que la media de la variable de interés es la misma en los diferentes grupos, en contraposición a la hipótesis alternativa de que al menos dos medias difieren de forma significativa.

Por lo tanto, ANOVA permite comparar múltiples medias, pero lo hace mediante el estudio de sus varianzas.

El funcionamiento básico consiste comparar la varianza de las medias por grupo (varianza explicada por la variable grupo, intervarianza) frente a la varianza promedio dentro de los grupos (la no explicada por la variable grupo, intravarianza).

El estadístico estudiado en el ANOVA, que sigue una distribución “F de Fisher/Snedecor”, es el cociente entre la varianza de las medias de los grupos y el promedio de la varianza dentro de los grupos.

Si se cumple la \(H_0\), el estadístico \(F\) adquiere el valor de 1 ya que la intervarianza será igual a la intravarianza.

Cuanto más diferentes son las medias de los grupos, mayor será la varianza entre esas medias en comparación al promedio de la varianza dentro de los grupos, obteniéndose valores de F superiores a 1.

ANOVA de un factor o de una vía

Esta prueba es una extensión de los t-test independientes para más de dos grupos.

Entendemos por factor a la propiedad o característica que nos permite distinguir entre sí a distintas poblaciones o grupos.

Los requisitos necesarios son:

• Observaciones aleatorias
• Independencia de las observaciones entre grupos
• Variable respuesta: cuantitativa continua
• Variable explicatoria o factor: cualitativa (> 2 grupos)
• Normalidad para cada grupo
• Homocedasticidad (homogeneidad de varianzas, asume que las varianzas dentro de cada grupo a comparar son homogéneas.)

Las hipótesis contrastadas en un ANOVA de una vía son:

• \(H_0\): No hay diferencias entre las medias de los diferentes grupos

• \(H_1\): Al menos un par de medias son significativamente distintas entre ellas.

Ejemplo de ANOVA con R

La concentración plasmática elevada de lipoproteínas de alta densidad (HDL) se acompaña de un menor riesgo de padecer cardiopatía coronaria.

Varios estudios sugieren que el ejercicio vigoroso eleva la concentración de HDL. Con el fin de investigar si el trote incrementa la concentración plasmática de HDL, G. Harley Hartung et al. cuantificaron la concentración de HDL en corredores de maratón, trotadores y varones sedentarios (35 a 66 años de edad).

La concentración promedio de HDL en estos últimos fue de 43,3 mg/100 ml con una desviación estándar de 14,2 mg/100 ml.

La media y desvío estándar de la concentración de HDL en los trotadores y maratonistas fueron de 58,1 y 17,2 mg/100 ml y de 64,8 y 14,1 mg/100 ml, respectivamente.

Si cada grupo constaba de 70 varones, comprobemos si existen o no existen diferencias significativas en la concentración promedio de HDL en los diversos grupos.

Vayamos paso a paso con el código en R:

Lectura y exploración de datos

El archivo Anova-HDL.txt es de texto plano y utiliza el punto y coma (;) como separador.

library(tidyverse)
library(skimr)
library(dlookr)

HDL <- read_csv2("Anova-HDL.txt")

Exploramos con skim de skimr

skim(HDL)

Data summary

Name	HDL
Number of rows	210
Number of columns	2
_______________________
Column type frequency:
character	1
numeric	1
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
grupo	0	1	8	11	0	3	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
hdl	0	1	55.39	17.65	14.7	43.62	55.95	66.65	96.5	▃▅▇▆▂

Observamos que la tabla de datos contiene dos variables y 210 observaciones.

Una llamada grupo de tipo caracter con tres categorías diferentes (n_unique = 3) y sin valores NA (n_missing = 0).

La otra con nombre hdl, numérica con valores entre 14,7 y 96,5 y sin NA´s.

Antes de continuar vamos a convertir la variable grupo de caracter a factor dándole cierto orden a sus categorías.

HDL <- HDL %>%
  mutate(grupo = factor(grupo, 
                    levels = c("Sedentario", "Trotador", "Maratonista")))

# el orden de los niveles se otorga al escribir sus niveles (en este ejemplo van en aumento de actividad fisica)

## verificamos con class() y levels()

class(HDL$grupo)

[1] "factor"

levels(HDL$grupo)

[1] "Sedentario"  "Trotador"    "Maratonista"

Graficamente podemos ver la distribución con un histograma y un boxplot.

HDL %>% 
  ggplot(aes(x = hdl)) + 
  geom_histogram(fill="mediumpurple1", binwidth = 10, color = "white")

HDL %>% 
  ggplot(aes(y = hdl)) + 
  geom_boxplot(fill="indianred")

La distribución global es simétrica y no hay valores atípicos.

Ahora veamos como es la distribución de los valores de HDL por grupo.

HDL %>% 
  ggplot(aes(x = hdl, fill = grupo)) + 
  geom_histogram(binwidth = 10, color = "white") +
  facet_grid(grupo ~ .) # variable grupo en filas

HDL %>% 
  ggplot(aes(y = hdl, x = grupo, fill = grupo)) + 
  geom_boxplot() 

Observamos que las distribuciones en cada grupo parecen simétricas y que a simple vista hay diferencias entre las medianas. La pregunta es si podemos afirmar que esas diferencias son estadísticamente significativas o pueden deberse al azar.

Análisis de supuestos

El análisis de los supuestos puede realizarse previo al test de ANOVA, dado que si no se cumplen no tiene mucho sentido seguir adelante. Sin embargo la forma habitual de comprobar que se satisfacen las condiciones necesarias es estudiando los residuos del modelo una vez generado el objeto ANOVA en R.

Normalidad

Dado que los grupos tienen mas de 50 observaciones empleamos el test de Kolmogorov-Smirnov con la corrección de Lilliefors. La función en R se llama lillie.test() y se encuentra en el paquete nortest. Si fuesen menos de 50 eventos por grupo podríamos utilizar el test Shapiro-Wilk (shapiro.test() en R base).

library(nortest)

HDL %>% 
  filter(grupo == "Sedentario") %>%  # filtramos grupo Sedentario
  pull(hdl) %>%  # con pull() extraemos datos de hdl como vector
  lillie.test()  # aplicamos test de bondad de ajuste

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  .
D = 0.096922, p-value = 0.1063

HDL %>% 
  filter(grupo == "Trotador") %>% 
  pull(hdl) %>% 
  lillie.test()

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  .
D = 0.061224, p-value = 0.7434

HDL %>% 
  filter(grupo == "Maratonista") %>% 
  pull(hdl) %>% 
  lillie.test()

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  .
D = 0.075981, p-value = 0.4053

En los tres grupos los test de hipótesis no muestran evidencias de falta de normalidad.

También podemos graficar con Q-Q Plots

HDL %>%
  filter(grupo == "Sedentario") %>% 
  plot_normality(hdl)

HDL %>%
  filter(grupo == "Trotador") %>% 
  plot_normality(hdl, col = "orange")

HDL %>%
  filter(grupo == "Maratonista") %>% 
  plot_normality(hdl, col = "purple")

La gráfica muestra que los puntos se distribuyen a lo largo de la línea, incluso para el grupo de sedentarios que se aleja un poco en sus colas.

Como conclusión podemos asumir el ajuste normal de las distribuciones.

Varianza constante entre grupos (homocedasticidad)

Usaremos el test de bondad de ajuste implementado para probar homocedasticidad de Bartlett.

bartlett.test(hdl ~ grupo, # utiliza sintaxis fórmula
              data = HDL) 

    Bartlett test of homogeneity of variances

data:  hdl by grupo
Bartlett's K-squared = 3.5696, df = 2, p-value = 0.1678

No hay evidencias significativas de falta de homocedasticidad (p valor > 0,05).

Análisis de varianza ANOVA

La función del lenguaje para el ANOVA es aov() y se encuentra en R base.

Utiliza sintaxis fórmula del tipo \(var\_cuanti \sim factor\) y su resultado es un objeto lista (modelo “aov-lm”).

Siempre que trabajemos con modelos, sea para el ANOVA como para las regresiones, conviene que asignemos los resultados a objetos que luego nos servirán para aplicar otras funciones.

En las siguientes dos líneas corremos el análisis ANOVA para el conjunto de datos y luego ejecutamos un resumen del análisis con summary() para obtener los resultados completos, incluido el p-valor.

resultado <- aov(formula = hdl ~ grupo, 
                 data = HDL)

summary(resultado)

             Df Sum Sq Mean Sq F value   Pr(>F)    
grupo         2  16934    8467   36.38 2.88e-14 ***
Residuals   207  48172     233                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La tabla de resultados muestra un p-valor muy menor a 0,05 sugiriendo que hay evidencias suficientes para considerar que al menos dos medias son distintas.

Comparaciones múltiples

Con el resultado anterior solo podemos decir que hay diferencias significativas entre al menos dos grupos de los tres que tenemos pero no podemos decir entre cuales sucede.

Existen diferentes algotirmos de comparaciones múltiples y correcciones implementados en R. Aquí vamos a aplicar la prueba de rango de Tukey cuya función es TukeyHSD().

La HSD de Tukey proviene de Diferencia Honestamente Significativa de Tukey es una técnica de comparaciones múltiples y de rangos a la vez. Se aplica para grupos equilibrados (mismo tamaño) y varianzas similares (homogeneidad de varianzas).

Compara los grupos entre ellos de dos en dos y calcula el intervalo de confianza de esas diferencias.

Es una prueba conservadora, dado que mantiene bajo el error de tipo I, sacrificando la capacidad de detectar diferencias existentes.

Si lo aplicamos sobre el objeto resultado con el modelo ANOVA realizado, obtenemos:

TukeyHSD(resultado)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = hdl ~ grupo, data = HDL)

$grupo
                           diff        lwr      upr     p adj
Trotador-Sedentario    14.77857  8.6913607 20.86578 0.0000001
Maratonista-Sedentario 21.49857 15.4113607 27.58578 0.0000000
Maratonista-Trotador    6.72000  0.6327893 12.80721 0.0264602

Hay mayor diferencia entre los grupos Sedentario-Maratonista y Trotador-Sedentario que en Trotador-Maratonista, aunque todos tiene valor-p ajustado significativo.

Una mejor forma de ver este resultado es en un gráfico:

plot(TukeyHSD(resultado))

La línea punteada vertical es el 0 de “no hay diferencias” y cada segmento de comparación de dos grupos está conformado por los intervalos de confianza. Se observa que ninguno de estos tocan el valor 0 y claramente se encuentran separados unos de otros.

Tamaño de efecto y potencia de la prueba

Los test de potencia permiten determinar la probabilidad de encontrar diferencias significativas entre las medias para un determinado nivel de significancia indicando el tamaño de los grupos.

El paquete pwr contiene la función pwr.anova.test() que necesita como argumentos el número de grupos, el número de observaciones por cada grupo, el tamaño del efecto y el nivel de significancia (habitualmente 0,05).

El tamaño de efecto comúnmente utilizado para el caso de ANOVA es eta cuadrado, que se calcula haciendo el cociente de la suma de cuadrados del efecto sobre la suma de cuadrados total. (los datos se pueden leer en la salida del summary(modelo))

library(pwr)

eta_cuadrado <- 16934/(16934 + 48172)

eta_cuadrado

[1] 0.2600989

## tamaño de efecto mediano convencional de Cohen 
cohen.ES(test = c("anov"), size = c("medium"))

     Conventional effect size from Cohen (1982) 

           test = anov
           size = medium
    effect.size = 0.25

pwr.anova.test(k = 3, n = 70, f = eta_cuadrado, sig.level = 0.05)

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 70
              f = 0.2600989
      sig.level = 0.05
          power = 0.9280577

NOTE: n is number in each group

Para el tamaño de efecto calculado (cercano al medio convencional propuesto por Cohen) la potencia del este de ANOVA efectuado es 92,8 %.

Análisis de residuos

Decíamos anteriormente que para poder dar por validos los resultados del ANOVA es necesario verificar que se satisfacen las condiciones analizando los residuos del modelo.

Describiremos este procedimiento en la Unidad III cuando trabajemos con modelos de regreión lineal.

Conclusión

Se ha realizado un ANOVA de 1 vía para comparar el valor de HDL en personas según tres grupos/niveles de actividad física. Se ha encontrado una diferencia estadísticamente significativa entre las medias de los grupos (F = 36,38, p < 0,01). El eta cuadrado fue de 0,26, lo que indica que la variable independiente (la actividad física) explica el 26% de la varianza total en los valores de HDL. La potencia de la prueba estadística para un nivel de significancia de 0,05 % fue de 92,8%.

Se realizó la prueba post hoc de Tukey HSD y se encontró que las diferencias entre la media de cada grupo fue estadísticamente significativa (p < 0,05). (parece que todas las medias de HDL de los grupos difieren entre ellas.)