Tutorial R: ANOVA en R

ANOVA de un factor y multifactorial

El Análisis de Varianza (ANOVA) es una técnica estadística que nos permite comparar las medias de tres o más grupos simultáneamente. A diferencia de la prueba t que solo compara dos grupos, ANOVA nos ayuda a determinar si existen diferencias significativas entre múltiples grupos.

ANOVA de un factor

El ANOVA de un factor (también llamado ANOVA unidireccional o one-way ANOVA) examina el efecto de una única variable independiente (factor) sobre una variable dependiente continua. Por ejemplo, podríamos analizar si diferentes métodos de enseñanza afectan las calificaciones de los estudiantes.

Veamos un ejemplo sencillo. Imaginemos que queremos comparar el rendimiento de tres variedades de plantas (A, B y C) en términos de altura:

# Creamos un conjunto de datos de ejemplo
set.seed(123)  # Para reproducibilidad
variedad_A <- rnorm(10, mean = 25, sd = 5)
variedad_B <- rnorm(10, mean = 30, sd = 5)
variedad_C <- rnorm(10, mean = 22, sd = 5)

# Organizamos los datos en un data frame
datos <- data.frame(
  altura = c(variedad_A, variedad_B, variedad_C),
  variedad = factor(rep(c("A", "B", "C"), each = 10))
)

# Visualizamos los primeros registros
head(datos)

Para realizar un ANOVA de un factor en R, utilizamos la función aov():

# Realizamos el ANOVA
modelo_anova <- aov(altura ~ variedad, data = datos)

# Vemos el resumen de resultados
summary(modelo_anova)

La salida nos mostrará:

• La suma de cuadrados entre grupos e intra-grupos
• Los grados de libertad
• La media cuadrática (suma de cuadrados dividida por grados de libertad)
• El valor F (cociente entre la varianza entre grupos y la varianza intra-grupos)
• El valor p (probabilidad de obtener un valor F igual o más extremo si la hipótesis nula fuera cierta)

Interpretación básica

La hipótesis nula en ANOVA es que todas las medias de los grupos son iguales. Si el valor p es menor que nuestro nivel de significancia (típicamente 0.05), rechazamos esta hipótesis y concluimos que al menos un grupo difiere significativamente de los demás.

# Obtenemos estadísticas descriptivas por grupo
aggregate(altura ~ variedad, data = datos, FUN = function(x) c(
  media = mean(x),
  desv_est = sd(x)
))

ANOVA multifactorial

El ANOVA multifactorial examina el efecto de dos o más factores sobre una variable dependiente. Esto nos permite no solo estudiar los efectos principales de cada factor, sino también sus interacciones.

Veamos un ejemplo con dos factores: variedad de planta (A, B, C) y tipo de fertilizante (1, 2):

# Creamos datos para un ANOVA de dos factores
set.seed(456)

# Creamos un diseño factorial completo
datos_multi <- expand.grid(
  variedad = factor(c("A", "B", "C")),
  fertilizante = factor(c("1", "2"))
)

# Añadimos 5 repeticiones por combinación
datos_multi <- datos_multi[rep(1:nrow(datos_multi), each = 5), ]

# Generamos alturas con efectos de ambos factores
datos_multi$altura <- with(datos_multi, {
  base <- ifelse(variedad == "A", 25, 
                ifelse(variedad == "B", 30, 22))
  fert_efecto <- ifelse(fertilizante == "1", 0, 5)
  # Añadimos interacción: el fertilizante 2 es especialmente efectivo con la variedad B
  interaccion <- ifelse(variedad == "B" & fertilizante == "2", 3, 0)
  
  rnorm(nrow(datos_multi), mean = base + fert_efecto + interaccion, sd = 3)
})

# Visualizamos los primeros registros
head(datos_multi)

Para realizar un ANOVA multifactorial, incluimos todos los factores y sus interacciones en la fórmula:

# ANOVA con dos factores e interacción
modelo_multi <- aov(altura ~ variedad * fertilizante, data = datos_multi)

# Vemos el resumen
summary(modelo_multi)

La notación variedad * fertilizante es equivalente a variedad + fertilizante + variedad:fertilizante, donde:

• variedad y fertilizante representan los efectos principales
• variedad:fertilizante representa la interacción entre ambos factores

Interpretación de efectos principales e interacciones

En un ANOVA multifactorial, interpretamos:

1. Efectos principales: Si cada factor por separado tiene un efecto significativo sobre la variable dependiente.
2. Interacciones: Si el efecto de un factor depende del nivel del otro factor.

# Estadísticas descriptivas por combinación de factores
aggregate(altura ~ variedad + fertilizante, data = datos_multi, 
          FUN = function(x) c(media = mean(x), desv_est = sd(x)))

Cuando existe una interacción significativa, debemos ser cautelosos al interpretar los efectos principales, ya que estos pueden variar según los niveles de otros factores.

Verificación de supuestos

Para que los resultados de ANOVA sean válidos, debemos verificar ciertos supuestos:

1. Normalidad de los residuos
2. Homogeneidad de varianzas entre grupos
3. Independencia de las observaciones

Podemos verificar estos supuestos con gráficos de diagnóstico:

# Gráficos de diagnóstico para el modelo ANOVA
par(mfrow = c(2, 2))
plot(modelo_anova)

Para verificar específicamente la homogeneidad de varianzas, podemos usar la prueba de Levene:

# Prueba de homogeneidad de varianzas
# Requiere el paquete car
if (!require(car)) {
  install.packages("car")
  library(car)
}

leveneTest(altura ~ variedad, data = datos)

Si los supuestos no se cumplen, podemos considerar:

• Transformaciones de la variable dependiente
• Pruebas no paramétricas como Kruskal-Wallis
• ANOVA con correcciones de heteroscedasticidad

ANOVA con medidas repetidas

Cuando las observaciones no son independientes (por ejemplo, mediciones repetidas en los mismos sujetos), necesitamos un ANOVA de medidas repetidas:

# Ejemplo simplificado de datos con medidas repetidas
datos_rep <- data.frame(
  sujeto = factor(rep(1:10, each = 3)),
  tiempo = factor(rep(c("pre", "durante", "post"), times = 10)),
  rendimiento = c(
    # Sujeto 1-10, tres mediciones cada uno
    rnorm(10, mean = 20, sd = 5),  # pre
    rnorm(10, mean = 25, sd = 5),  # durante
    rnorm(10, mean = 22, sd = 5)   # post
  )
)

# ANOVA de medidas repetidas
modelo_rep <- aov(rendimiento ~ tiempo + Error(sujeto/tiempo), data = datos_rep)
summary(modelo_rep)

En este modelo:

• tiempo es nuestro factor de interés
• Error(sujeto/tiempo) especifica que las mediciones están anidadas dentro de cada sujeto

El ANOVA de medidas repetidas tiene en cuenta la correlación entre mediciones del mismo sujeto, lo que aumenta la potencia estadística del análisis.

Comparaciones post-hoc con TukeyHSD

Cuando realizamos un ANOVA y encontramos diferencias significativas entre grupos, sabemos que al menos un grupo difiere de los demás, pero no sabemos exactamente cuáles grupos son diferentes entre sí. Para responder a esta pregunta necesitamos realizar pruebas post-hoc, siendo una de las más utilizadas la prueba de Tukey de diferencia honestamente significativa (TukeyHSD).

La prueba TukeyHSD nos permite realizar comparaciones por pares entre todos los grupos mientras controla la tasa de error global, evitando el problema de las comparaciones múltiples que ocurriría si realizáramos múltiples pruebas t independientes.

Aplicación de TukeyHSD en R

Aplicar la prueba de Tukey en R es muy sencillo. Primero necesitamos un modelo ANOVA y luego aplicamos la función TukeyHSD():

# Continuando con el ejemplo anterior de variedades de plantas
# Primero creamos el modelo ANOVA
modelo_anova <- aov(altura ~ variedad, data = datos)

# Aplicamos la prueba de Tukey
resultados_tukey <- TukeyHSD(modelo_anova)

# Visualizamos los resultados
print(resultados_tukey)

La salida de la prueba TukeyHSD incluye para cada par de grupos:

• La diferencia entre las medias de los grupos
• Los límites inferior y superior del intervalo de confianza (por defecto al 95%)
• El valor p ajustado para esa comparación específica

Interpretación de los resultados

Veamos cómo interpretar los resultados de la prueba TukeyHSD:

# Ejemplo de salida (los valores exactos dependerán de tus datos)
#   Tukey multiple comparisons of means
#     95% family-wise confidence level
#
# Fit: aov(formula = altura ~ variedad, data = datos)
#
# $variedad
#     diff        lwr       upr     p adj
# B-A  5.12  0.7234567  9.516543 0.0214563
# C-A -2.87 -7.2665433  1.526543 0.2543210
# C-B -7.99 -12.3865433 -3.593457 0.0008765

Para interpretar estos resultados:

1. Observamos la columna diff que muestra la diferencia entre medias
2. Revisamos la columna p adj que contiene el valor p ajustado
3. Si el valor p es menor que nuestro nivel de significancia (típicamente 0.05), concluimos que existe una diferencia significativa entre ese par de grupos

En el ejemplo anterior, podríamos concluir que:

• Las variedades B y A difieren significativamente (p = 0.021)
• Las variedades C y A no muestran diferencias significativas (p = 0.254)
• Las variedades C y B difieren significativamente (p < 0.001)

TukeyHSD para ANOVA multifactorial

Cuando trabajamos con un ANOVA multifactorial, podemos aplicar TukeyHSD para examinar las diferencias entre los niveles de cada factor:

# Continuando con el ejemplo multifactorial
modelo_multi <- aov(altura ~ variedad * fertilizante, data = datos_multi)

# Aplicamos TukeyHSD
tukey_multi <- TukeyHSD(modelo_multi)

# Visualizamos los resultados para cada factor y la interacción
print(tukey_multi)

En este caso, obtendremos comparaciones para:

• Los niveles del factor variedad
• Los niveles del factor fertilizante
• Las combinaciones de la interacción variedad:fertilizante

Manejo de interacciones significativas

Cuando existe una interacción significativa entre factores, la interpretación de las comparaciones post-hoc se vuelve más compleja. En estos casos, es recomendable:

1. Analizar primero la interacción
2. Si la interacción es significativa, realizar comparaciones dentro de cada nivel de un factor

Podemos hacer esto creando subconjuntos de datos:

# Si la interacción es significativa, podemos analizar el efecto de la variedad
# para cada tipo de fertilizante por separado

# Para fertilizante tipo 1
datos_fert1 <- subset(datos_multi, fertilizante == "1")
modelo_fert1 <- aov(altura ~ variedad, data = datos_fert1)
TukeyHSD(modelo_fert1)

# Para fertilizante tipo 2
datos_fert2 <- subset(datos_multi, fertilizante == "2")
modelo_fert2 <- aov(altura ~ variedad, data = datos_fert2)
TukeyHSD(modelo_fert2)

Ajuste del nivel de confianza

Por defecto, TukeyHSD utiliza un nivel de confianza del 95%, pero podemos ajustarlo según nuestras necesidades:

# Cambiar el nivel de confianza al 99%
TukeyHSD(modelo_anova, conf.level = 0.99)

Un nivel de confianza más alto (99%) resultará en intervalos más amplios y un enfoque más conservador, mientras que un nivel más bajo (90%) producirá intervalos más estrechos pero aumentará la probabilidad de falsos positivos.

Ejemplo práctico: Comparación de tratamientos médicos

Veamos un ejemplo más aplicado. Supongamos que estamos comparando tres tratamientos para reducir la presión arterial:

# Creamos datos de ejemplo
set.seed(789)
tratamiento_A <- rnorm(15, mean = -10, sd = 3)  # Reducción media de 10 mmHg
tratamiento_B <- rnorm(15, mean = -15, sd = 4)  # Reducción media de 15 mmHg
tratamiento_C <- rnorm(15, mean = -8, sd = 3)   # Reducción media de 8 mmHg

# Organizamos los datos
datos_presion <- data.frame(
  reduccion = c(tratamiento_A, tratamiento_B, tratamiento_C),
  tratamiento = factor(rep(c("A", "B", "C"), each = 15))
)

# ANOVA
modelo_presion <- aov(reduccion ~ tratamiento, data = datos_presion)
summary(modelo_presion)

# Si el ANOVA es significativo, procedemos con TukeyHSD
tukey_presion <- TukeyHSD(modelo_presion)
print(tukey_presion)

En este caso, la interpretación sería en términos de eficacia clínica:

• Si el tratamiento B muestra una diferencia significativa respecto a A y C, podríamos concluir que es el más efectivo
• Si no hay diferencia significativa entre A y C, podríamos considerar otros factores como costo o efectos secundarios para elegir entre ellos

Representación tabular de resultados

Para presentar los resultados de forma más ordenada, podemos convertir la salida de TukeyHSD en un data frame:

# Convertir los resultados de Tukey a un data frame
tukey_df <- as.data.frame(tukey_presion$tratamiento)

# Añadir una columna con los pares de comparación
tukey_df$comparacion <- rownames(tukey_df)

# Reordenar las columnas
tukey_df <- tukey_df[, c("comparacion", "diff", "lwr", "upr", "p adj")]

# Añadir una columna que indique si la diferencia es significativa
tukey_df$significativo <- ifelse(tukey_df$`p adj` < 0.05, "Sí", "No")

# Mostrar la tabla
print(tukey_df)

Esta representación facilita la interpretación y comunicación de los resultados, especialmente cuando trabajamos con muchos grupos o necesitamos incluir los resultados en un informe.

Limitaciones de TukeyHSD

Aunque TukeyHSD es una herramienta muy útil, tiene algunas limitaciones:

• Asume que los tamaños de muestra son iguales o similares entre grupos
• Requiere que se cumplan los supuestos del ANOVA (normalidad, homogeneidad de varianzas)
• No es la mejor opción cuando hay una gran disparidad en las varianzas entre grupos

En casos donde no se cumplen estos supuestos, podemos considerar alternativas como:

• La prueba de Games-Howell (disponible en paquetes como PMCMRplus)
• Comparaciones con corrección de Bonferroni o Holm
• Métodos no paramétricos como la prueba de Dunn después de Kruskal-Wallis

# Ejemplo de uso de corrección de Bonferroni con pairwise.t.test
pairwise.t.test(datos_presion$reduccion, datos_presion$tratamiento, 
                p.adjust.method = "bonferroni")

Visualización de resultados ANOVA

La visualización de los resultados de un análisis ANOVA es fundamental para comunicar eficazmente los hallazgos estadísticos y facilitar su interpretación. R ofrece diversas herramientas gráficas que nos permiten representar tanto los datos originales como los resultados del análisis de varianza de manera clara y efectiva.

Diagramas de caja (Boxplots)

Los diagramas de caja son una excelente manera de visualizar la distribución de los datos por grupos y detectar visualmente posibles diferencias entre ellos:

# Continuando con el ejemplo de variedades de plantas
# Creamos un boxplot básico
boxplot(altura ~ variedad, data = datos, 
        col = c("lightblue", "lightgreen", "lightpink"),
        main = "Altura por variedad de planta",
        ylab = "Altura (cm)",
        xlab = "Variedad")

# Añadimos una línea horizontal con la media global
abline(h = mean(datos$altura), col = "red", lty = 2)

Para ANOVA multifactorial, podemos crear boxplots agrupados:

# Boxplot para datos multifactoriales
boxplot(altura ~ variedad:fertilizante, data = datos_multi,
        col = rep(c("lightblue", "lightgreen", "lightpink"), each = 2),
        main = "Altura por variedad y tipo de fertilizante",
        ylab = "Altura (cm)",
        xlab = "Combinación variedad:fertilizante")

# Mejoramos las etiquetas del eje x
axis(1, at = 1:6, 
     labels = c("A:1", "A:2", "B:1", "B:2", "C:1", "C:2"), 
     cex.axis = 0.8)

Gráficos de medias con barras de error

Los gráficos de medias con barras de error son ideales para visualizar las diferencias entre grupos y su significancia estadística:

# Calculamos medias y errores estándar por grupo
library(dplyr)

resumen <- datos %>%
  group_by(variedad) %>%
  summarise(
    media = mean(altura),
    ee = sd(altura)/sqrt(n()),
    ic_inf = media - qt(0.975, n()-1) * ee,
    ic_sup = media + qt(0.975, n()-1) * ee
  )

# Creamos el gráfico de barras con barras de error
barplot_pos <- barplot(resumen$media, 
                      names.arg = resumen$variedad,
                      col = c("lightblue", "lightgreen", "lightpink"),
                      ylim = c(0, max(resumen$ic_sup) * 1.1),
                      main = "Altura media por variedad",
                      ylab = "Altura (cm)",
                      xlab = "Variedad")

# Añadimos barras de error (intervalos de confianza 95%)
arrows(x0 = barplot_pos, 
       y0 = resumen$ic_inf, 
       y1 = resumen$ic_sup, 
       angle = 90, code = 3, length = 0.1)

Gráficos de interacción

Para ANOVA multifactorial, los gráficos de interacción son esenciales para visualizar cómo un factor modifica el efecto de otro:

# Calculamos medias por combinación de factores
medias_interaccion <- with(datos_multi, 
                          tapply(altura, list(variedad, fertilizante), mean))

# Creamos el gráfico de interacción
interaction.plot(x.factor = datos_multi$fertilizante,
                trace.factor = datos_multi$variedad, 
                response = datos_multi$altura,
                type = "b",
                col = c("blue", "red", "green"),
                pch = c(16, 17, 18),
                legend = TRUE,
                xlab = "Tipo de fertilizante",
                ylab = "Altura media (cm)",
                main = "Interacción entre variedad y fertilizante")

Un gráfico de interacción con líneas no paralelas sugiere la presencia de interacción entre factores. Cuanto más se crucen las líneas, más fuerte es la interacción.

Visualización de resultados post-hoc

Podemos visualizar los resultados de las pruebas post-hoc de Tukey para facilitar su interpretación:

# Visualizamos gráficamente los resultados de TukeyHSD
plot(TukeyHSD(modelo_anova), las = 1)

Para una visualización más personalizada:

# Extraemos los resultados de Tukey
tukey_result <- TukeyHSD(modelo_anova)$variedad

# Creamos un gráfico de intervalos de confianza
plot(1:nrow(tukey_result), tukey_result[,"diff"],
     ylim = range(tukey_result[,c("lwr","upr")]),
     pch = 16, xaxt = "n",
     xlab = "Comparación", ylab = "Diferencia",
     main = "Intervalos de confianza de Tukey")

# Añadimos los nombres de las comparaciones
axis(1, at = 1:nrow(tukey_result), labels = rownames(tukey_result), las = 2)

# Añadimos los intervalos de confianza
arrows(1:nrow(tukey_result), tukey_result[,"lwr"], 
       1:nrow(tukey_result), tukey_result[,"upr"], 
       angle = 90, code = 3, length = 0.1)

# Añadimos una línea horizontal en cero
abline(h = 0, lty = 2, col = "red")

# Coloreamos según significancia
points(1:nrow(tukey_result), tukey_result[,"diff"],
       pch = 16,
       col = ifelse(tukey_result[,"p adj"] < 0.05, "red", "black"))

En este gráfico, los intervalos que no cruzan la línea de cero (en rojo) representan diferencias estadísticamente significativas.

Visualización con ggplot2

El paquete ggplot2 ofrece opciones más avanzadas y estéticamente agradables para visualizar resultados ANOVA:

# Cargamos ggplot2
library(ggplot2)

# Boxplot con ggplot2
ggplot(datos, aes(x = variedad, y = altura, fill = variedad)) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.5) +  # Añade puntos individuales
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "black") +
  labs(title = "Altura por variedad de planta",
       x = "Variedad",
       y = "Altura (cm)") +
  theme_minimal() +
  scale_fill_brewer(palette = "Pastel1")

Para visualizar interacciones:

# Gráfico de interacción con ggplot2
ggplot(datos_multi, aes(x = fertilizante, y = altura, color = variedad, group = variedad)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Interacción entre variedad y fertilizante",
       x = "Tipo de fertilizante",
       y = "Altura media (cm)") +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

Visualización de supuestos ANOVA

Es importante visualizar también los supuestos del ANOVA para validar nuestros resultados:

# Gráficos de diagnóstico para verificar supuestos
par(mfrow = c(2, 2))
plot(modelo_anova)

Para una visualización más específica de la normalidad:

# Gráfico QQ para verificar normalidad de residuos
residuos <- residuals(modelo_anova)
qqnorm(residuos)
qqline(residuos)

# Histograma de residuos
hist(residuos, breaks = 10, main = "Histograma de residuos",
     xlab = "Residuos", col = "lightblue")

Visualización de la homogeneidad de varianzas

Podemos visualizar la homogeneidad de varianzas entre grupos:

# Gráfico de dispersión de residuos vs. valores ajustados
plot(fitted(modelo_anova), residuos,
     xlab = "Valores ajustados", ylab = "Residuos",
     main = "Residuos vs. Valores ajustados")
abline(h = 0, lty = 2, col = "red")

# Boxplot de residuos por grupo
boxplot(residuos ~ datos$variedad,
        main = "Residuos por variedad",
        xlab = "Variedad", ylab = "Residuos")
abline(h = 0, lty = 2, col = "red")

Visualización de efectos principales

Para ANOVA multifactorial, es útil visualizar los efectos principales de cada factor:

# Efectos principales con ggplot2
# Para el factor variedad
ggplot(datos_multi, aes(x = variedad, y = altura)) +
  stat_summary(fun = mean, geom = "bar", fill = "lightblue") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Efecto principal: Variedad",
       x = "Variedad",
       y = "Altura media (cm)") +
  theme_minimal()

# Para el factor fertilizante
ggplot(datos_multi, aes(x = fertilizante, y = altura)) +
  stat_summary(fun = mean, geom = "bar", fill = "lightgreen") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  labs(title = "Efecto principal: Fertilizante",
       x = "Tipo de fertilizante",
       y = "Altura media (cm)") +
  theme_minimal()

Visualización integrada de resultados

Finalmente, podemos crear un panel de visualización que integre diferentes aspectos de nuestro análisis ANOVA:

# Configuramos un panel de múltiples gráficos
par(mfrow = c(2, 2))

# 1. Boxplot de los datos originales
boxplot(altura ~ variedad, data = datos, 
        col = c("lightblue", "lightgreen", "lightpink"),
        main = "Distribución por grupo",
        ylab = "Altura (cm)")

# 2. Gráfico de medias con intervalos de confianza
barplot_pos <- barplot(resumen$media, 
                      names.arg = resumen$variedad,
                      col = c("lightblue", "lightgreen", "lightpink"),
                      ylim = c(0, max(resumen$ic_sup) * 1.1),
                      main = "Medias por grupo",
                      ylab = "Altura (cm)")
arrows(x0 = barplot_pos, y0 = resumen$ic_inf, y1 = resumen$ic_sup, 
       angle = 90, code = 3, length = 0.1)

# 3. Gráfico QQ de residuos
qqnorm(residuos, main = "Normalidad de residuos")
qqline(residuos)

# 4. Resultados de Tukey
plot(tukey_result[,"diff"], 
     ylim = range(tukey_result[,c("lwr","upr")]),
     pch = 16, xaxt = "n",
     main = "Comparaciones de Tukey",
     xlab = "", ylab = "Diferencia")
axis(1, at = 1:nrow(tukey_result), labels = rownames(tukey_result), las = 2, cex.axis = 0.7)
arrows(1:nrow(tukey_result), tukey_result[,"lwr"], 
       1:nrow(tukey_result), tukey_result[,"upr"], 
       angle = 90, code = 3, length = 0.1)
abline(h = 0, lty = 2, col = "red")

Esta visualización integrada proporciona una visión completa del análisis, desde los datos originales hasta los resultados de las comparaciones post-hoc, facilitando la interpretación y comunicación de los hallazgos estadísticos.