El Análisis de la Varianza o ANOVA es una de las técnicas que se utilizan para comparar los valores medios de varias poblaciones estadísticas. Consiste en descomponer toda la variabilidad que interviene en un muestreo o experimento, en las fuentes de variación que han sido controladas o consideradas por el investigador. La manera de obtener los datos, es lo que condiciona el tipo de descomposición.

ANOVA DE UN FACTOR (one-way ANOVA)

Cuando se toman muestras aleatorias en distintas poblaciones claramente definidas (machos y hembras; 3 estadios larvales de un crustáceo; las 4 estaciones del año), los datos obtenidos tienen 2 posibles fuentes de variación: la que existe entre las poblaciones, y la que naturalmente hay dentro de cada una de ellas. La primera es un efecto fijo y sistemático (que puede o no existir), provocado por el hecho de que un individuo está en una o en otra población. La segunda es aleatoria, y se genera por todos los factores que hacen que los individuos de una misma población no sean iguales.

Ejemplo I. Se tienen 3 hospedadores diferentes (Lenguado, Pejerrey y Pescadilla) para el tercer estadio larval de un nematodo parásito. La variable Y es la longitud (mm) de la larva. Para cada hospedador se puede definir una población con media \(\mu_i\) y varianza \(\sigma_i^2\). Las larvas que pertenecen a una misma población no tienen por qué tener igual longitud. Esto se debe a una serie de factores aleatorios que determinan la magnitud de \(\sigma_i^2\) (varianza dentro de la especie hospedadora). Cuando pasamos de una población a otra, podría existir una variación originada por el cambio de la especie a parasitar; lo que provocaría una modificación sistemática de los valores medios (\(\mu_i\)). Se toma una muestra aleatoria de tamaño \(n_i\) en cada población, y se estiman los parámetros poblacionales correspondientes con las medias y varianzas muestrales.

Este diseño puede representarse a través de un modelo que describe cada uno de los elementos que conforman la variable: \[\LARGE y_{ij} = \mu_i + \epsilon_{ij}\]

donde:
\(y_{ij}\) es la observación “j” de la población “i”

\(\mu_i\) es la media de la población “i”

\(\epsilon_{ij}\) es una variable aleatoria debida a los factores no controlados en el muestreo

Vamos a analizar sólo el caso en que se toma igual número de observaciones por muestra (muestreo balanceado), aunque en este Modelo el desbalance sólo complica un poco los cálculos y no es un obstáculo para su aplicación.

Nomenclatura y fórmulas:

\(I\): número de poblaciones

\(n\): tamaño de la muestra en cada población

\(N\): número total de observaciones

\(y_{ij}\): observación “j” de la población “i”

\(\bar{y_i}\): promedio de las observaciones de la población “i”

\(S_i^2\): varianza de las observaciones de la población “i”

\(y\) : promedio de todas las observaciones

En el ANOVA, las varianzas se denominan cuadrados medios. Mientras que sus numeradores son llamados sumas de cuadrados, los denominadores constituyen los grados de libertad.

CM: cuadrado medio   SC: suma de cuadrados   gl: grados de libertad       \(CM = \frac{{SC}}{gl}\)

\(\mathrm{SC_{Total}}\) = mide la variación entre cada dato y el promedio de todas las observaciones
\(\mathrm{SC_{Entre}}\) = mide la variación entre los promedios obtenidos en cada población y el promedio total
\(\mathrm{SC_{Dentro}}\) = mide la variación entre cada dato y el promedio obtenido de la población a la que pertenecen

Aditividad de las Sumas de Cuadrados y los grados de libertad

La \(\mathrm{SC_{Total}}\) se puede descomponer de manera aditiva en las \(\mathrm{SC_{Entre}}\) y \(\mathrm{SC_{Dentro}}\); del mismo modo que sus grados de libertad:

\[\underbrace{\sum_{i=1}^{I}\sum_{j=1}^{n}(y_{ij} - \bar{y})^2}_{\mathrm{SC_{Total}}} = \underbrace{n\times\sum_{i=1}^{I}(\bar{y}_{i} - \bar{y})^2}_{\mathrm{SC_{Entre}}} + \underbrace{\sum_{i=1}^{I}\sum_{j=1}^{n}(y_{ij} - \bar{y_{i}})^2}_{\mathrm{SC_{Dentro}}}\]                                              \(\mathrm{\textbf{gl}:\hspace{2.25cm} N - 1 \hspace{1.1cm} = \hspace{1.2cm} I - 1 \hspace{1.2cm} + \hspace{1.1cm} N - I}\)

Los grados de libertad son la diferencia entre los términos independientes a comparar y las restricciones con las que se comparan.

Cuando se realiza un experimento, se aplican al azar distintos tratamientos a una serie de unidades experimentales (u.e.) que, inicialmente, pertenecían a una misma población. Luego, cada una de ellas pasa a pertenecer a una población diferente, según el tratamiento que le tocó en suerte. Aquí, las 2 posibles fuentes de variación son: la que existe entre los tratamientos, y la que hay dentro de cada uno de ellos.

Como antes, la primera es un efecto fijo, provocado por el hecho de que una u.e. recibe un determinado tratamiento, y no otro. La segunda es aleatoria, y se genera por todos los factores que hacen que las u.e. no sean iguales, a pesar de estar bajo un mismo tratamiento.

Ejemplo II. Se tienen “N” ratones (u.e.) infectados con un virus, a los que se les aplicarán por sorteo 3 dosis diferentes de un suero: 0, 10 y 20 ml (tratamientos). De este modo quedan \(\mathrm{n_i}\) u.e. en cada tratamiento, que constituyen muestras aleatorias de las correspondientes poblaciones. Sus parámetros se estiman a partir de las medias y las varianzas muestrales.

La variable Y es carga viral al cabo de 1 semana de la aplicación de los sueros. Para cada dosis se puede definir una población con media \(\mu_i\) y varianza \(\sigma_i^2\). Los ratones que recibieron una misma dosis de suero no presentan la misma carga viral, debido a una serie de factores aleatorios que determinan la magnitud de \(\mu_i\) (varianza dentro de cada dosis). Cuando pasamos de un tratamiento a otro, podría existir una variación originada por el cambio de dosis; lo que provocaría una variación sistemática en los valores medios poblacionales (\(\mu_i\)).

Este diseño experimental puede representarse a través de un modelo que describe cada uno de los elementos que conforman la variable: \[\LARGE y_{ij} = \mu + \alpha_i + \epsilon_{ij}\]

donde:
\(y_{ij}\) es la observación “j” del tratamiento “i”

\(\mu\) es la media de todo el experimento

\(\alpha_i\) es el efecto que el tratamiento “i” produce respecto de \(\mu\)

\(\epsilon_{ij}\) es una variable aleatoria debida a los factores no controlados del experimento

Si bien la nomenclatura es otra, el modelo es el mismo que el expresado para el muestreo de poblaciones (basta reemplazar ”\(\mu_i\)“ por “\(\mu\) + \(\alpha_i\)“). Por lo tanto, la descomposición aditiva de las sumas de cuadrados y de los grados de libertad es exactamente igual.

Para que la prueba estadística que se realiza en el ANOVA tenga validez, se deben cumplir determinados supuestos:

  1. Cada población (o tratamiento) debe tener una distribución Normal
  2. Todas las poblaciones (o tratamientos) deben tener la misma varianza (homocedasticidad):

\(\hspace{15cm}\sigma_1^2 = \sigma_2^2 =...= \sigma_I^2 = \sigma^2\)

Esta varianza común es la que estima el \(\mathrm{CM_{Dentro}}\).

El \(\mathrm{CM_{Entre}}\) estima, además, los cambios que se producen en el valor medio por aplicar distintos tratamientos o pertenecer a diferentes poblaciones (es una estimación de \(\sigma^2\), más los efectos debidos al factor Controlado).

Esto se puede ver explícitamente a través de la Esperanza de los CM (lo que se espera que estén estimando estas expresiones). \[E(CM_{Entre}) = \sigma^2 + n\times\frac{\sum{\alpha_i^2}}{I - 1}\] \[\hspace{-2.4cm}E(CM_{Dentro}) = \sigma^2\]

La Hipótesis Nula plantea que no hay diferencias entre las medias poblacionales: \[\mathrm{H_0}:\ \mu_1 = \mu_2 =...= \mu_I = \mu\]

versus la Hipótesis Alternativa de que por lo menos 1 de estas igualdades no se cumple.

Si utilizáramos la otra nomenclatura, plantearíamos: \(\mathrm{H_0}:\ \alpha_1 = \alpha_2 =...= \alpha_I = 0\), es decir: que los tratamientos no producen cambios en la media total del experimento.

Cuando \(\mathrm{H_0}\) es cierta, el \(\mathrm{CM_{Entre}}\) es sólo otro estimador de \(\sigma^2\) (el término de la derecha en la expresión de la esperanza se hace “cero”). Pero cuando \(\mathrm{H_0}\) es falsa, el \(\mathrm{CM_{Entre}}\) contiene, además, los efectos que provocan los tratamientos. Por lo tanto, será MAYOR que el \(\mathrm{CM_{Dentro}}\).

Una forma de comparar ambos CM es a través de su cociente:   \(\Large{F = \frac{CM_{Entre}}{CM_{Dentro}}}\)

“F” es un estadístico que, bajo \(\mathrm{H_0}\), tiene una distribución conocida:

\(\mathscr{F}\) de Snedecor. Es positiva y asimétrica. Su forma depende de 2 parámetros: los grados de libertad del numerador y los grados de libertad del denominador (la cantidad de información que tiene cada CM sobre \(\sigma^2\)).

En este caso:

\[ \mathrm{F} \sim \mathscr{F} (I-1; N-I)\]

Si este cociente es muy grande, tendríamos razones para sospechar que el numerador contiene algo más que \(\sigma^2\): las diferencias entre las medias poblacionales (o entre los tratamientos). Por lo tanto, deberíamos rechazar \(\mathrm{H_0}\).

En gráfico siguiente se muestra la distribución del estadístico F para (6;30) grados de libertad, junto con los valores de F para un nivel de significancia del 95% y 99% y las zonas de rechazo (áreas sombreadas).

  • Si el estadístico F supera al valor crítico de tabla al 5%, rechazamos \(\mathrm{H_0}\) y afirmamos que existe por lo menos alguna diferencia entre las medias poblacionales (entre los tratamientos), con una probabilidad de error menor del 5% (p < 0.05). También se suele decir que se hallaron diferencias significativas, y se agrega un asterisco (*) al lado del estadístico.

  • Si el estadístico F supera al valor crítico de tabla al 1%, rechazamos \(\mathrm{H_0}\) y afirmamos que existe por lo menos alguna diferencia entre las medias poblacionales (entre los tratamientos), con una probabilidad de error menor del 1% (p < 0.01). También se suele decir que se hallaron diferencias altamente significativas, y se agregan dos asteriscos (**) al lado del estadístico.

  • Si el estadístico F no supera al valor crítico de tabla al 5%, no rechazamos \(\mathrm{H_0}\) y, por lo tanto, no hallamos diferencias entre las medias poblacionales (entre los tratamientos), trabajando con un nivel de significación mayor del 5% (p > 0.05). También se suele decir que las diferencias fueron no significativas, y se agregan las siglas (ns) al lado del estadístico. Por supuesto que si, al aumentar el \(\alpha\) seguimos aceptando \(\mathrm{H_0}\), podemos decir que no se hallaron diferencias, trabajando con un nivel de significación mayor (p > \(\alpha\)).

Ejemplo de parásitos de tres peces

Veamos el ejemplo I, donde teniamos 3 peces que hospedan un estadio larval de Nematoda. Se tomaron 9 larvas al azar para cada una de las especies. Se midió la longitud (mm) de cada larva.

Resultados Lenguado Pejerrey Pescadilla
4.60 4.75 4.99
4.75 5.50 5.72
5.00 6.17 6.12
5.22 6.27 6.91
Datos 5.55 6.36 7.06
5.69 6.80 7.52
5.81 6.87 7.65
6.40 7.00 8.27
6.66 7.07 8.58
Media 5.52 6.31 6.98
Varianza 0.50 0.59 1.41

\[\mathrm{I = 3 \hspace{0.4cm} n = 9 \hspace{0.4cm} N = 27 \hspace{0.4cm} \bar{y} = 6.27}\]

Resolvamos este ejemplo aplicando R. En primer lugar carguemos los datos. Podes descargar el data frame parasitos.csv del siguiente link.

## cargamos los datos
data = read.csv("parasitos.csv", header = TRUE, sep = ";") 
head(data) ## mostramos las primeras 6 filas
##   Longitud  Pescado
## 1     4.60 Lenguado
## 2     4.75 Lenguado
## 3     5.00 Lenguado
## 4     5.22 Lenguado
## 5     5.55 Lenguado
## 6     5.69 Lenguado
str(data)  ## vemos la estructura de los datos
## 'data.frame':    27 obs. of  2 variables:
##  $ Longitud: num  4.6 4.75 5 5.22 5.55 5.69 5.81 6.4 6.66 4.75 ...
##  $ Pescado : Factor w/ 3 levels "Lenguado","Pejerrey",..: 1 1 1 1 1 1 1 1 1 2 ...

Como se puede ver, nuestros datos están formados por 27 observaciones (9 larvas por cada uno de los peces) y dos variables Longitud (tamaño del III estadio larval del parásito) y Pescado, cada uno de los tres pescados donde se extrajo la larva (lenguado, pejerrey y pescadilla). Tener presenten que la última variable es del tipo factor (categórica) con 3 niveles (cada uno de los pescados). A lo que nosotros nos interesa saber es si la longitud de los parásitos difiere entre los tres pescados. Por lo tanto Y = longitud de la larva (mm) y X = especie de pescado. Nuestra \(\mathrm{H_0}\) es: \[\mathrm{H_0}:\ \mu_{Lenguado} = \mu_{Pejerrey} = \mu_{Pescadilla}\]

Como tenemos una sola variable podemos aplicar un ANOVA simple para ver si existen diferencias entre las medias de las longitudes de las larvas entre los 3 pescados.

Un buen primer paso para empezar sería realizar algún gráfico exploratorio para corroborar como se comportan los datos.

## gráfico exploratorio de f(x) = Longitud(Pescado)
library(ggplot2)

data = read.csv("parasitos.csv", header = T, sep = ";")

p <- ggplot(data, aes(factor(Pescado), Longitud, fill = Pescado)) + 
  stat_boxplot(geom ='errorbar')
  
p + geom_boxplot() + theme_bw() + 
  labs(
    title = "Boxplot de las tres poblaciones de pescados",
    x = ""
    ) + 
  theme(legend.position = "none")

Como se observa en la figura anterior, parecería ser que existe una diferencia entre las medias de longitud de las larvas entre los diferentes pescados. Con ANOVA podemos ver si esa diferencia es estadísticamente significativa. Pues hagamoslo.

## usamos la función lm() para realizar un anova (recordar que un anova se puede interpretar 
## como una regresión lineal donde las variables X son categóricas)
datos = read.csv("parasitos.csv", header = TRUE, sep = ";") # cargamos los datos
fit.larvas = lm(Longitud ~ Pescado, data = datos) #hacemos la regresión
anova(fit.larvas) ## una vez aplicamos lm(), la función anova() nos devuelve la TABLA DE ANOVA
## Analysis of Variance Table
## 
## Response: Longitud
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## Pescado    2  9.6138  4.8069  5.7765 0.008954 **
## Residuals 24 19.9716  0.8322                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## En la Tabla de ANOVA, Pescado corresponde a la variación "ENTRE" y Residuals a la variación "DENTRO"

De la Tabla deducimos que como el esadístico F (5.7765) supera el \(\mu\) del 99%, entonces rechazamos \(\mathrm{H_0}\). Podemos afirmar que hay diferencias altamente significativas entre los hospedadores (p < 0.01). Pero…¿son todas las poblaciones diferentes? No lo sabemos, podemos afirmar que por lo menos una de las poblaciones es diferente a alguna de las otras. Para poder saber quién difiere con quien o si todas difieren entre sí debemos hacer comparaciones de a pares (contrastes).

La forma más rápida de realizar contrastes (alguno) es pidiendole a R que nos devuelva el summary() de la función lm().

## usamos el summary() para realizar los contrastes
datos = read.csv("parasitos.csv", header = TRUE, sep = ";") # cargamos los datos
fit.larvas = lm(Longitud ~ Pescado, data = datos) #hacemos la regresión

summary(fit.larvas) ## pedimos los contrastes
## 
## Call:
## lm(formula = Longitud ~ Pescado, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.990 -0.645  0.050  0.615  1.600 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.5200     0.3041  18.153 1.59e-15 ***
## PescadoPejerrey     0.7900     0.4300   1.837  0.07861 .  
## PescadoPescadilla   1.4600     0.4300   3.395  0.00239 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9122 on 24 degrees of freedom
## Multiple R-squared:  0.325,  Adjusted R-squared:  0.2687 
## F-statistic: 5.776 on 2 and 24 DF,  p-value: 0.008954

La función hace I - 1 contrastes, es decir, tantas comparaciones como le permite los grados de libertad del \(\mathrm{SC_{Entre}}\), en este caso 2.

  • (Intercept) corresponde a la media del grupo base (el primero que aparece en el data frame), en nuestro caso corresponde a Lenguado. Ver boxplot.
  • PescadoPejerrey contabiliza la diferencia entre la media de Lenguado y la media de Pejerrey. Para obtener la media de Pejerrey solamente le sumamos a la media de lenguado (el Incercept) ésta diferencia (0.79). Ver boxplot.
  • PescadoPescadilla corresponde a la diferencia entre Pescadilla y Lenguado (el Intercept). Ver boxplot.

De esta tabla podemos ver que la media de Pescadilla difiere de la de Lenguado, (p < 0.001). Sin embargo, Pejerrey no difiere estadísticamente del Lenguado (p > 0.05). Pero, ¿Pejerrey difiere de Pescadilla? No lo sabemos, la función lm() solo hace comparaciones entre el Incercept y cada una de las poblaciones (recorda que la función solo puede hacer tantas comparaciones como grados de libertad haya).

Para poder realizar todas las comparaciones posibles debemos hacer contrastes a posteriori, uno de los más usado es la prueba de significancia honesta de Tukey (post-hoc). Para ello primero debemos realizar la prueba de ANOVA con la función aov() que funciona de la misma forma que anova(). Luego aplicamos la función TukeyHSD().

datos = read.csv("parasitos.csv", header = TRUE, sep = ";") # cargamos los datos
fit.larvas2 = aov(Longitud ~ Pescado, data = datos) ## similar a anova() -chequealo-

contrastes = TukeyHSD(fit.larvas2) ## hacemos TODOS los contrastes posibles
contrastes
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Longitud ~ Pescado, data = datos)
## 
## $Pescado
##                     diff       lwr      upr     p adj
## Pejerrey-Lenguado   0.79 -0.283898 1.863898 0.1791591
## Pescadilla-Lenguado 1.46  0.386102 2.533898 0.0065005
## Pescadilla-Pejerrey 0.67 -0.403898 1.743898 0.2828008

Ahora sí podemos ver las comparaciones de todos con todos. Como se ve, si bien las diferencias son las mismas en ambas funciones, los p-values son diferentes. La única diferencia que dió significativa es la que hay entre Pescadilla y Lenguado (p < 0.05), el resto dio no significativo. La ventaja que tiene esta prueba es que nos ofrece un intervalo de confianza. La diferencia entre Pescadilla y Pejerrey (la que nos faltaba) dio no significativa.

Ejercicios para ANOVA de un factor

A continuación planteo un par de ejercicios para ANOVA simple.

  1. En 1942 Geoffrey Beall realizó una serie de experimentos en Canadá (Chatham, Ontario) para evaluar el efecto de distintos insecticidas sobre distintas especies de polillas y mariposas (Lepidoptera). Realizó un total de 7 experimentos a campo donde delimitó determinadas unidades experimentales (imaginemonos que son quadrats), aplicó el insecticida y contabilizó el número de insectos vivos luego de un determinado intervalo de tiempo. Nosotros vamos a trabajar con el séptimo experimento, donde aplicó 6 tipos diferentes de insecticidad (A-F) a (claramente) 7 poblaciones estadísticas de la pollilla Manduca quinquemaculata. Los datos a utilizar ya vienen incluidos en el propio R bajo la biblioteca InsectSprays. A partir de estos averigue si existen diferencias estadísticas que nos permitan afirmar si existen diferencias significativas entre los distintos insecticidas. Una aclaración importante, al recopilar los datos Beall se dió cuenta que las varianzas entre las distintas poblaciones no eran iguales por lo que no se podría realizar un Análisis de la Varianza (violación al supuesto de Homocedasticidad), por lo que propuso que esto se podría solucionar aplicando la raiz cuadrada a los datos del conteo del número de polillas vivas.

  2. El departamento de Psicología de una Universidad de Castilla-La Mancha (España) ha realizado un estudio sobre hábitos, preferencias y satisfacción sexual en estudiantes universitarios. Utilizaremos los datos que recogieron en sus encuestas para conocer si existen diferencias entre la frecuencia mensual de relaciones sexuales de estudiantes universitarios pertenecientes a tres carreras universitarias diferentes. Los datos se pueden descargar del siguiente link (archivo sexo.csv).

  3. El 11 de marzo de 2020 la OMS declara a la COVID-19 (enfermedad producida por el virus SARS-CoV-2) como una pandemia. Desde ese entonces los casos de contagios y muertes producidos por este virus han aumentado drásticamente. Hoy en día existen varias vacunas en circulación que atacan directamente a dicho virus, por lo que los casos graves de COVID-19 han empezado a mermar. Una de las variables que condiciona en cierta medida los contagios por COVID-19 es la edad, y varios artículos afirman que las personas más propensas a infectarse son los adultos mayores ¿Es esto realmente así? El presente ejercicio tiene el objetivo de evaluar si existen diferencias entre los contagiados y fallecidos de COVID-19 de acuerdo a su grupo etario. Para evaluar esto se utilizarán datos provenentes de Colombia recopilados por el Instituto Nacional de la Salud en el periodo comprendido entre el 6 y 10 de julio de 2020. Trataremos de reproducir los resultados presentes en el artíiculo original. Realiza un ANOVA para evaluar si existen diferencias significativa entre los contagiados y fallecidos de acuerdo a los grupos etarios. Posteriormente llevá a cabo comparaciones múltiples post-hoc mediante la prueba de significancia honesta de Tukey (el artículo original utiliza la prueba de Duncan). Podes descargar los datos en este link (archivo COVID.csv).

  4. Spartina alterniflora es una especie de graminea invasora reconocida como una gran ingeniera de ecosistemas, debido a que es capaz de crear, alterar o mantener el hábitat e incidir en la composición de las comunidades, principalmente por modificación de las propiedades abióticas del ecosistema. Si Spartina tiene efectos positivos o negativos en las comunidades es discutido. Algunos autores apoyan la idea de que el desarrollo de marismas de especies exóticas provoca una disminución de la abundancia y diversidad de invertebrados bentónicos en relación a la planicie de marea adyacente, en contraposición a otro grupo que predice un efecto positivo. En el año 2017 una bióloga de la Universidad Nacional del Sur (Argentina) realizó un estudio en el estuario de Bahía Blanca con el fin de evaluar el efecto que tienen las marismas de S. alterniflora sobre la estructura comunitaria de la zona. En particular se compararon zonas con y sin marismas a lo largo de un gradiente altitudinal. A continuación se presenta parte de ese estudio donde se evalua la variación en la abundancia de Heleobia australis en 4 zonas diferentes del estuario. Las zonas rotuladas como A y B corresponden a áreas vegetadas y C y D a zonas con suelo desnudo. Los datos pueden descargarse del siguiente link (archivo Heleobia.csv). ¿Que se puede decir del efecto de la marisma sobre la abundancia de esta especie de molusco?