ANOVA en R: aov(), la tabla ANOVA y Tukey HSD

Por qué ANOVA y no un montón de tests t

Un test t compara dos medias. Con tres o más grupos, lo tentador es hacer un test t con cada par - y esa es exactamente la trampa. Cada test ejecutado al nivel 0.05 conlleva un 5% de riesgo de falso positivo, y los riesgos se acumulan: con 4 grupos son 6 tests por pares y aproximadamente un 26% de probabilidad de al menos un resultado "significativo" espurio; con 5 grupos (10 tests), aproximadamente un 40%. Fabricarías descubrimientos a partir de puro ruido.

El ANOVA (ANalysis Of VAriance, análisis de varianza) lo soluciona planteando una pregunta con un valor p: ¿son todas las medias de grupo iguales, o al menos una difiere? A pesar del nombre, compara medias - solo que lo hace analizando la varianza: si los promedios de los grupos están más dispersos de lo que el ruido dentro de los grupos puede explicar, algo real está ocurriendo.

ANOVA de una vía con aov()

PlantGrowth está hecho para esto: pesos de plantas bajo un control y dos condiciones de tratamiento. Ajusta con aov() y luego - importante - imprime la tabla con summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Lee la fórmula como "¿depende weight de group?". La variable de agrupación debe ser un factor - PlantGrowth$group ya lo es, pero si tus grupos están codificados como números (dosis, IDs de lote), envuélvelos: aov(y ~ factor(dose), ...). De lo contrario aov() ajusta en silencio una recta de regresión a través de los códigos de grupo en lugar de comparar las medias de los grupos - modelo equivocado, sin mensaje de error.

Cómo leer la tabla ANOVA

La tabla tiene dos filas - el factor y los residuos - y cinco columnas. Celda a celda:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df - grados de libertad. La fila del factor recibe grupos − 1 (3 grupos → 2); la fila de los residuos recibe observaciones − grupos (30 − 3 = 27). Una comprobación rápida de que R vio el diseño que pretendías.
• Sum Sq - la variación, repartida en dos montones. La fila group es el montón entre grupos: cuán lejos se sitúan las medias de los grupos respecto a la media global. Residuals es el montón dentro de los grupos: cuánto varían las plantas alrededor de la media de su propio grupo. Juntos suman la variación total de los datos.
• Mean Sq - cada Sum Sq dividido entre sus Df, convirtiendo montones de variación en tasas comparables por grado de libertad. El Mean Sq residual (0.389) es el nivel de ruido.
• F value - la razón: el Mean Sq de group sobre el Mean Sq de Residuals (1.8832 / 0.3886 ≈ 4.85). Si todas las medias de grupo fueran realmente iguales, esta razón ronda el 1. Cuanto mayor es F, más superan las diferencias entre grupos lo que el ruido puede explicar.
• Pr(>F) - el valor p: la probabilidad de un F así de grande si las tres medias verdaderas fueran idénticas. Aquí 0.016 - lo bastante pequeño al nivel convencional de 0.05 para rechazar "todas iguales". Como siempre, no es la probabilidad de que la nula sea cierta, y no dice nada sobre qué grupos difieren ni en cuánto.

Ese último punto es la limitación crucial: un F significativo dice "alguna diferencia, en alguna parte". Nada más.

¿Qué grupos difieren? TukeyHSD()

La pregunta de seguimiento necesita un test post-hoc. La Diferencia Honestamente Significativa de Tukey contrasta cada par manteniendo la tasa de error a nivel de familia en el 5% sobre todas las comparaciones a la vez:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

Una fila por par, cuatro números por fila:

• diff - la diferencia estimada de medias (el segundo grupo nombrado menos el primero).
• lwr, upr - el intervalo de confianza del 95% a nivel de familia para esa diferencia.
• p adj - el valor p, ya ajustado por hacer tres comparaciones.

Para PlantGrowth: trt2-trt1 muestra una diferencia de aproximadamente 0.87 con p adj ≈ 0.012 y un intervalo alejado del cero - el tratamiento 2 crece más que el tratamiento 1. Las dos filas de tratamiento contra control (trt1-ctrl, trt2-ctrl) tienen intervalos que abarcan el cero y p adj muy por encima de 0.05 - ninguno de los tratamientos se distingue del control en esta muestra. La regla general refleja la de los intervalos de confianza en todas partes: intervalo que excluye el cero ⇔ p adj por debajo de 0.05.

Fíjate en cómo el F global dijo "existe una diferencia" mientras que Tukey la localiza exactamente en un par - la estructura en dos pasos es todo el diseño: un test global honesto, y luego trabajo detectivesco por pares debidamente corregido.

ANOVA de dos vías: dos factores y su interacción

Con dos variables de agrupación, una sola llamada contrasta ambas - además de si interactúan. ToothGrowth cruza el tipo de suplemento (supp) con la dosis (0.5, 1, 2 mg - numérica, así que necesita factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) se expande en tres efectos, una fila de la tabla cada uno:

• supp - promediando sobre las dosis, ¿importa el tipo de suplemento? (Sí: p ≈ 0.0002.)
• factor(dose) - promediando sobre los suplementos, ¿importa la dosis? (Rotundamente: p es minúsculo.)
• supp:factor(dose) - la interacción: ¿cambia el efecto del suplemento dependiendo de la dosis? Aquí p ≈ 0.022 - sí cambia. El zumo de naranja supera al ácido ascórbico en dosis bajas, pero la brecha se cierra en la dosis de 2 mg.

Una interacción significativa es una etiqueta de advertencia sobre los efectos principales: "el tipo de suplemento importa" solo es cierto en promedio, y el promedio esconde una historia dependiente de la dosis. Cuando la interacción es significativa, describe las combinaciones (un TukeyHSD(fit2, "supp:factor(dose)") o una tabla de medias de grupo mediante resúmenes por grupos) en lugar de reportar los efectos principales por sí solos. Usa + en lugar de * solo cuando quieras deliberadamente un modelo sin interacción.

Supuestos - y el plan B basado en rangos

Las matemáticas del ANOVA se apoyan en tres supuestos, en orden decreciente de negociabilidad:

• Independencia - las observaciones no se influyen entre sí. Nada arregla una violación a posteriori; es una propiedad del diseño del estudio.
• Varianzas iguales entre grupos - el Mean Sq residual es una única estimación agrupada del ruido, así que los grupos deberían ser más o menos igual de ruidosos. Comprobación de una línea: bartlett.test(weight ~ group, data = PlantGrowth) (un valor p grande significa que no hay evidencia de varianzas desiguales).
• Residuos aproximadamente normales - es lo que menos importa con diseños balanceados y tamaños de grupo decentes; échale un vistazo a un histograma o boxplot de los residuos mediante residuals(fit).

Cuando los datos están muy sesgados, son ordinales o están plagados de valores atípicos, la alternativa basada en rangos al ANOVA de una vía es kruskal.test(weight ~ group, data = PlantGrowth) - el hermano para muchos grupos de wilcox.test(). Es una línea, y cambia algo de potencia por robustez.

Lo que te llevas

• El ANOVA plantea una pregunta sobre las medias de 3+ grupos con un solo valor p - el remedio a la multiplicidad de tests t.
• fit <- aov(y ~ group, data = df); summary(fit) - y la variable de agrupación debe ser un factor.
• El valor F es la señal entre grupos sobre el ruido dentro de los grupos; un Pr(>F) pequeño significa "alguna diferencia en alguna parte", y nada más.
• TukeyHSD(fit) encuentra qué pares difieren, con la corrección por comparaciones múltiples incorporada.
• a * b ajusta dos factores más su interacción; una interacción significativa significa que los efectos principales no cuentan la historia por sí solos.
• Comprueba las varianzas iguales (bartlett.test) y la normalidad de los residuos; kruskal.test() es el plan B basado en rangos.

A continuación: de comparar medias de grupos a modelar una relación - regresión lineal con lm().