Práctica 5: El modelo ANOVA (ANalysis Of VAriance) de uno y dos factores

1. Ejemplo introductorio

Se quiere comparar la velocidad de tres algoritmos A, B y C para procesar imágenes. Se les pasa una batería de 10 imágenes, las procesan y se anota el tiempo que han tardado con cada una. Los tiempos medios por algoritmo son:

Algoritmo	Tiempo medio
A	5.1
B	5.0
C	4.9

El algoritmo C ha sido “globalmente” más rápido con la muestra de imágenes, ¿pero lo seguirá siendo a largo plazo? Veamos dos situaciones extremas con esos valores medios:

CONCLUSIÓN: SE BUSCA TÉCNICA “OBJETIVA” (JUSTIFICADA ESTADÍSTICAMENTE) QUE PERMITA COMPARAR Y DECIDIR SOBRE SI LOS PROMEDIOS DE MÁS DE 2 POBLACIONES SON SIGNIFICATIVAMENTE DISTINTOS O NO (A PARTIR DE OBSERVAR MUESTRAS DE ELLAS). ES DECIR, QUE PERMITA DECIDIR ENTRE:

\(H_0\): \(\mu_A = \mu_B = \mu_C\) (todas las medias iguales)
\(H_1\): no \(H_0\) (alguna media es distinta)

2. Objetivos del ANOVA

OBJETIVO 1:
- Decidir si:
  - el valor medio de \(Y\) depende del nivel del factor \(X\) o no
  - el factor \(X\) influye sobre el valor medio de \(Y\)
- Es decir, resolver el contraste:
  - \(H_0\): \(\mu_1 = \mu_2 = \ldots = \mu_a\) (todas las medias iguales)
  - \(H_1\): no \(H_0\) (alguna media distinta)
- Interpretación:
  - Aceptar \(H_0\) es “todas las medias son iguales”, es decir, “el factor \(X\) no influye sobre el valor medio de \(Y\)”, es decir “\(Y\) no depende de \(X\)”.
  - Rechazar \(H_0\) es “alguna media es distinta”, es decir “\(X\) sí influye sobre \(Y\)”, es decir “\(Y\) sí depende de \(X\)”.
OBJETIVO 2:
- Si se rechaza \(H_0\) interesa hacer un ránking de los niveles del factor \(X\) según el valor medio que aportan a la respuesta \(Y\).
  - Comparaciones a posteriori. Hay muchos métodos: LSD de Fisher, HSD de Tuckey,…
  - Grupos homogéneos. Se codifican los niveles de \(X\) con letras (a, b,…). Letras compartidas por niveles del factor \(X\) significa que las medias de \(Y\) no son lo bastante diferentes.

3. El modelo ANOVA de un factor \(Y = \mu_X + N(0,\sigma^2)\)

\(Y\): variable numérica de gran interés.
\(X\): factor que podría influir sobre el valor final de \(Y\).
- \(x_1\), \(x_2\), \(\ldots\), \(x_a\): niveles del factor \(X\) (\(a\) = número de niveles).
Parámetros: \(\mu_1\), \(\mu_2\), \(\ldots\), \(\mu_a\) y \(\sigma^2\)
Muestreo:
- Se toma \(X=x_i\)
- Se observa, \(y_{ij} = \mu_i + e_{ij}\) para \(j=1,2,\ldots\), donde
  - \(\mu_i\) es el valor medio de \(Y\) provocado por el nivel \(X=x_i\) del factor,
  - \(e_{ij}\) es un valor aleatorio, sigue el modelo \(N(0,\sigma^2)\), es independiente de todos los valores muestreados.
Un ejemplo de muestra: \(\begin{array}{c|cccc} X & Y & & & \\ \hline x_1 & y_{11} & y_{12} & \cdots & y_{1n_1} \\ x_2 & y_{21} & y_{22} & \cdots & y_{2n_2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_a & y_{a1} & y_{a2} & \cdots & y_{an_a} \end{array}\),
En el ejemplo de los algoritmos:
- La variable de interés es \(Y\) = “tiempo de procesado”.
- El factor que puede influir sobre el tiempo es \(X\) = “algoritmo”.
  - Tiene 3 niveles, “A”, “B” y “C”.
- Desconocemos si el tiempo de cada algoritmo sigue el modelo normal, se puede chequear.
- Desconocemos si los 10 tiempos de cada algoritmo son independientes entre sí, e independientes de los tiempos de los otros algoritmos, pero podemos suponerlo, y luego ya veremos si es compatible con los resultados.
El \(p\)-valor del estadístico \(F\) de la tabla ANOVA ayuda a decidir si el efecto del factor \(X\) es significativo o no.

4. El modelo ANOVA de dos factores

4.1. Sin interacción \(Y = \mu_A + \mu_B + N(0,\sigma^2)\)

\(Y\): variable numérica de gran interés.
\(A\): un factor que podría influir sobre el valor final de \(Y\).
- \(A_1\), \(A_2\), \(\ldots\), \(A_a\): niveles del factor \(A\) (\(a\) = número de niveles).
\(B\): otro factor que podría influir sobre el valor final de \(Y\).
- \(B_1\), \(B_2\), \(\ldots\), \(B_b\): niveles del factor \(B\) (\(b\) = número de niveles).
Parámetros:
- \(\mu_{A_1}\), \(\mu_{A_2}\), \(\ldots\), \(\mu_{A_a}\) (correspondientes al factor \(A\))
- \(\mu_{B_1}\), \(\mu_{B_2}\), \(\ldots\), \(\mu_{B_b}\) (correspondientes al factor \(B\))
- \(\sigma^2\) la dispersión de los datos de \(Y\) en cada combinación de niveles de los factores
Muestreo:
- Se toma \(A=A_i\), \(B=B_j\)
- Se observa, \(y_{ijk} = \mu_{A_i} + \mu_{B_j} + e_{ijk}\) para \(k=1,2,\ldots\), donde
  - \(\mu_{A_i}\) es el valor medio de \(Y\) provocado por el nivel \(A=A_i\) del factor \(A\),
  - \(\mu_{B_j}\) es el valor medio de \(Y\) provocado por el nivel \(B=B_j\) del factor \(B\),
  - \(e_{ijk}\) es un valor aleatorio, sigue el modelo \(N(0,\sigma^2)\), es independiente de todos los valores muestreados.
Un ejemplo de muestra: \(\begin{array}{c|cccc} & B & & & \\ A & B_1 & B_2& \cdots & B_b\\ \hline A_1 & [y_{111}, \ldots, y_{11n_{11}}] & [y_{121}, \ldots, y_{12n_{12}}] & \cdots & [y_{1b1}, \ldots, y_{1bn_{1b}}] \\ A_2 & [y_{211}, \ldots, y_{21n_{21}}] & [y_{221}, \ldots, y_{22n_{22}}] & \cdots & [y_{2b1}, \ldots, y_{2bn_{2b}}] \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ A_a & [y_{a11}, \ldots, y_{a1n_{a1}}] & [y_{a21}, \ldots, y_{a2n_{a2}}] & \cdots & [y_{ab1}, \ldots, y_{abn_{ab}}] \end{array}\).
El ANOVA plantea dos contrastes, donde cada uno chequea el efecto (principal) de cada factor sobre la respuesta \(Y\). El \(p\)-valor de cada estadístico \(F\) ayuda a decidir sobre la significatividad del efecto de cada factor sobre la respuesta \(Y\).

4.2. Con interacción \(Y = \mu_A + \mu_B + \gamma_{(A,B)}+ N(0,\sigma^2)\)

El planteamiento es el mismo, con los dos factores \(A\) y \(B\).
El efecto “conjunto” de una combinación de niveles \(A\) y \(B\) NO ES “la suma de los efectos marginales”, es decir:
- \(y_{ijk} = \mu_{A_i} + \mu_{B_j} + \gamma_{A_i, B_j} + e_{ijk}\)
- con \(\gamma_{A_i, B_j} \neq 0\) al menos para alguna combinación de niveles \(A_i\), \(B_j\).
El ANOVA resuelve tres contrastes, uno para cada factor, y otro para la posible interacción, sobre la respuesta \(Y\). El \(p\)-valor de cada estadístico \(F\) ayuda a decidir sobre la significatividad del efecto de cada factor, y el de la interacción, sobre la respuesta \(Y\).
La interacción se puede chequear también gráficamente.

5. El modelo ANOVA con R

5.1. Los datos

Los datos deben ser en dos variables (una para \(Y\) y otra para \(X\)) o tres (una para \(Y\), otra para el factor \(A\) y otra para el factor \(B\)). Pueden ser:
- Columnas de una misma hoja de datos (mejor)
- Vectores sueltos (ok, pero se podría formar una hoja de datos con ellos)
A veces hay varios vectores con datos de \(Y\), un vector “por cada nivel” del factor \(X\). En ese caso conviene “fabricar” la hoja de datos que tenga:
- Una columna con todos los datos de \(Y\): concatenar todos los vectores de datos de \(Y\) con la función c().
- Otra columna con los niveles del factor \(X\) que acompañan a cada dato de \(Y\): repetir “cada nivel” de \(X\) (número o palabra) tantas veces como datos de \(Y\) asociados a ese nivel, y concatenar todos al final. Por ejemplo c(rep("x1", 10), rep("x2", 30),...).
- Atención! la columna (o columnas) de la variable factor (\(X\)) debe ser de tipo factor en R. Comprobarlo con la función str() aplicada a la hoja de datos, y si no lo es, forzarlo con la función factor(). Por ejemplo: data.frame(X=factor(vectornivelesX), Y=vectordatosY)

5.2. Chequeo gráfico de los datos

Diagramas de puntos o de caja son muy informativos de cómo puede resultar la decisión final.
Función boxplot(formula, data) realiza cajas de la variable dependiente separada por niveles de la variable (o variables) independientes.
Argumentos:
- formula: expresión de la forma y ~ x o bien y ~ x1 + x2 donde y es la etiqueta de la columna numérica (dependiente) y x (o x1 y x2) es la etiqueta de la columna factor (independiente) que hace la separación en grupos.
- data: hoja de datos con las columnas referidas en la formula.
Ejemplos:
- boxplot(formula=Sepal.Length ~ Species, data=iris)
- boxplot(formula=vocab~nativeBorn, data=carData::GSSvocab)
- boxplot(formula=vocab~nativeBorn+gender, data=carData::GSSvocab)

5.3. Cálculo de estadísticos por grupos

Tenemos una hoja de datos con una columna cuantitativa de la que queremos obtener estadísticos por grupos que vienen determinados por otras columnas.
Función aggregate(formula, data, FUN)
Argumentos:
- formula: expresión de la forma y ~ x o bien y ~ x1 + x2 donde y es la etiqueta de la columna numérica (dependiente) y x (o x1 y x2) es la etiqueta de la columna factor (independiente) que hace la separación en grupos.
- data: hoja de datos con las columnas referidas en la formula.
- FUN: función a calcular sobre los datos de la columna y para cada grupo indicado por la columna x (o columnas x1 y x2)
Ejemplos:
- aggregate(formula=Sepal.Lengty ~ Species, data=iris, FUN="mean")
- aggregate(formula=Sepal.Lengty ~ Species, data=iris, FUN="sd")
- aggregate(formula=vocab~nativeBorn+gender, data=carData::GSSvocab, FUN="mean")

5.4. La tabla ANOVA y la decisión del contraste

Es una tabla que ayuda a calcular el estadístico \(F\) y su \(p\)-valor para decidir el contraste de igualdad de medias.
Se consigue en dos pasos: primero aov() y luego summary()
Función aov( formula, data,... ):
- Argumentos:
  - formula: expresión que indica, primero, la variable depemndiente, después el símbolo ~ (que indica la dependencia), y después la (o las) variables factor, de la forma:
    - y ~ x, en el modelo de un factor
    - y ~ x1 + x2, en el modelo de dos factores sin interacción,
    - y ~ x1 * x2, en el modelo de dos factores con interacción,
  - data: si los datos están en dos vectores suelto, no utilizar. Si los datos están en dos columnas de una hoja de datos, poner el nombre de la hoja de datos.
- Devuelve: un objeto tipo lista, y por pantalla, las sumas de cuadrados (inter e intra) y sus grados de libertad, además de la estimación (insesgada) de \(\sigma\)
Función anova(object):
- Argumento: object, el objeto devuelto por aov()
- Devuelve: la tabla ANOVA (como hoja de datos), con las columnas de: grados de libertad (Df), las sumas de cuadrados (Sum Sq) inter (con el nombre de cada factor) e intra (con el nombre Residuals), las medias de cuadrados (Mean Sq), el valor (o valores) del estadístico de contraste (F) para cada factor, y su \(p\)-valor (Pr(>F)). Si se quiere acceder al \(p\)-valor exacto, poner anova(object)[1,5] (etc.)
Decisión del contraste:
- “Rechazar \(H_0\) si \(p\)-valor \(< \alpha\)”
- Recuerda:
  - Aceptar \(H_0\) significa…
    - Aceptar \(\mu_1 = \mu_2 = \cdots = \mu_a\), que es lo mismo que
    - Aceptar que \(Y\) se comporta igual, sea cual sea el valor de \(X\), que es lo mismo que
    - Aceptar que \(Y\) NO depende del factor \(X\), que es lo mismo que
    - Aceptar que \(X\) no ejerce influencia sobre \(Y\).
  - Rechazar \(H_0\) significa…
    - Lo contrario de todo lo anterior
    - En este caso se motiva continuar el análisis estableciendo un ranking de niveles del factor \(X\) (ver sección 4.4)

5.5. Chequeo gráfico de la interacción (dos factores)

Se puede inspeccionar gráficamente la ausencia o presencia de interacción.
Función interaction.plot(x.factor, trace.factor, response,...).
Argumentos:
- x.factor: columna de datos del factor (cuyos niveles formarán el eje x).
- trace.factor: columna de datos del otro factor.
- response: columna de datos de la variable respuesta \(Y\).
Si no hay interacción, las “curvas” de medias para cada nivel de un factor, fijando el otro factor, forman un perfil, que varía “paralelamente” cuando se cambia el nivel del otro factor.

5.6. Adecuación del modelo

Para poder aplicar ANOVA y confiar en sus resultados, las \(Y\)’s deben cumplir las condiciones del modelo ANOVA (errores “normales”, independientes entre sí, independientes del nivel del factor, y con varianza común en todos los niveles).
Cuando hay pocos datos por nivel, es imposible detectar nada, y se puede aceptar todo. Cuando hay muchos datos, se pueden comprobar las condiciones.
Normalidad de los datos de todos los niveles:
- Se pueden representar histogramas y ver si encajan en “campanas de Gauss”.
- Se pueden lanzar contrastes de normalidad shapiro.test() y ver si se acepta la hipótesis nula, que es la normalidad
- Se puede lanzar el ANOVA con R (aov()), hacerle un plot(), cuya figura Normal Q-Q (la segunda), relaciona los cuantiles teóricos de la normal con los de los errores estimados y:
  - Si los puntos quedan “cerca” de la diagonal, entonces, aparentemente, se está cumpliendo la normalidad de los errores.
  - Si “bastantes” puntos se alejan de la diagonal, entonces, aparentemente, está fallando la normalidad de los errores.
Igualdad de varianzas en todos los niveles:
- Se puede contrastar de manera objetiva: contraste de Bartlett de igualdad de varianzas
  - Función bartlett.test(formula, data) o bien bartlett.test(x,...):
    - formula: la misma que en aov()
    - data: el mismo que aov()
    - x,...: los vectores con los datos de \(Y\) para los diversos niveles del factor \(X\)
  - La hipótesis nula es la igualdad de varianzas en todos los niveles
- Se puede lanzar el ANOVA con R (aov()), hacerle un plot(), cuya figura Residuals vs Factor Levels (la cuarta) muestra las desviaciones de los datos respecto de la media del nivel, permitiendo comparar las varianzas (amplitud vertical de los puntos).

5.7. Comparaciones a posteriori

Si la tabla ANOVA conduce a aceptar que \(Y\) no depende de \(X\), ya está todo terminado.
En caso contrario, interesa conocer qué niveles de \(X\) dan mayor media de \(Y\). Hacemos un ranking de los niveles de \(X\).
Método LSD de Fisher:
- Función LSD.test(y, trt, alpha, console), perteneciente al paquete agricolae, que se debe cargar (e instalar si no lo está).
- Argumentos:
  - y: objeto devuelto por la función aov().
  - trt: nombre (etiqueta entrecomillada) del factor (la \(X\)) o vector de nombres de los factores si hay varios.
  - alpha: nivel de significación deseado (\(0.05\) por defecto).
  - console: poner a TRUE si se quiere ver en la R Console el resultado.
- Devuelve: un objeto de tipo lista con estadísticos y tablas,
  - Una tabla con los niveles del factor y, para cada uno,
    - la media de la variable \(Y\),
    - error estándar,
    - número de observaciones,
    - intervalo de confianza, etc.;
  - Estadísticos varios;
  - LSD: valor de la mínima distancia significativa, y
  - Tabla con
    - el ranking: código (a, b, etc.) asignado por valor de media,
    - nivel del factor,
    - media de \(Y\) para dicho nivel.
    - Puede aparece un nivel con dos códigos distintos (cuando no se puede distinguir del nivel superior ni del inferior, pero ellos sí se distinguen entre sí).

6. Ejercicios evaluables

Cambia el siguiente bloque a eval=TRUE para que cargue los datos necesarios.

# despeja esta linea cuando tengas el archivo
load("mt1021-2021-la-5-anova.RData")
ls()

Problema 1 (25%)

La variable tiempos recoge los tiempos de realización de una misma tarea informática de varios operadores (de características muy similares) bajo 3 sistemas operativos, para comparar sobre cuál es más rápida.

Realiza una gráfica “sencilla” que compare los tiempos por cada sistema operativo sobre el tiempo.

# Escribe aquí el codigo R de lo que se pide