Práctica 4: El modelo ANOVA (ANalysis Of VAriance)
Pablo Gregori - Universitat Jaume I de Castellón
MT1021 Ampliación de Estadística e Investigación Operativa
1. Ejemplo motivador
EJEMPLO 1: Tres trabajadores de una cadena de montaje van a comparar su efectividad. Se recoge el tiempo que tarda cada trabajador en realizar su tarea sobre las siguentes 10 unidades que pasan por su punto de la cadena.
A = c(4.44, 4.77, 6.56, 5.07, 5.13, 6.72, 5.46, 3.73, 4.31, 4.55)
B = c(6.22, 5.36, 5.40, 5.11, 4.44, 6.79, 5.50, 3.03, 5.70, 4.53)
C = c(3.93, 4.78, 3.97, 4.27, 4.37, 3.31, 5.84, 5.15, 3.86, 6.25)Tenemos que decidir si hay algún trabajador “más rápido” o no. Para ello valoramos a cada uno con su tiempo medio:
| Trabajador | Tiempo medio |
|---|---|
| A | 5.074 |
| B | 5.208 |
| C | 4.573 |
Por lo tanto, el más rápido ha sido el C.
Comparemos gráficamente los tiempos de cada trabajador (dibujamos nube de puntos y diagrama de caja superpuestos):
boxplot(A, B, C, names=c('A', 'B', 'C'), xlab='Trabajador', ylab='Tiempo')
points(x=rep(1,10), y=A); points(x=1, y=mean(A), pch=16, col='red')
points(x=rep(2,10), y=B); points(x=2, y=mean(B), pch=16, col='red')
points(x=rep(3,10), y=C); points(x=3, y=mean(C), pch=16, col='red')
Se ve que los tiempos de cada trabajador fluctúan, cada uno en torno a su media.
- PREGUNTA: Si volvieramos a medir los tiempos de los trabajadores sobre las 10 próximas unidades:
- ¿Saldrían exactamente los mismos tiempos? ¿Por qué?
- ¿Por quién apostarías que sería más rápido en promedio? ¿Por qué?
- ¿Cuánto confiarías en tu apuesta? ¿Por qué?
FIN EJEMPLO 1
EJEMPLO 2: Misma situación que en Ejemplo 1, pero con otros tres trabajadores distintos.
D = c(5.01, 5.04, 5.22, 5.07, 5.08, 5.24, 5.11, 4.94, 5.00, 5.02)
E = c(5.18, 5.09, 5.09, 5.06, 5.00, 5.23, 5.10, 4.86, 5.12, 5.01)
F = c(4.51, 4.59, 4.51, 4.54, 4.55, 4.45, 4.70, 4.63, 4.50, 4.74)De nuevo comparamos su rapidez usando los tiempos medios:
| Trabajador | Tiempo medio |
|---|---|
| D | 5.073 |
| E | 5.074 |
| F | 4.572 |
Por lo tanto, el más rápido ha sido el F.
Vemos ahora las fluctuaciones de los tiempos de cada uno en el diagrama de caja:
boxplot(D, E, F, names=c('D', 'E', 'F'), xlab='Trabajador', ylab='Tiempo')
points(x=rep(1,10), y=D); points(x=1, y=mean(D), pch=16, col='red')
points(x=rep(2,10), y=E); points(x=2, y=mean(E), pch=16, col='red')
points(x=rep(3,10), y=F); points(x=3, y=mean(F), pch=16, col='red')
- PREGUNTA: Si volvieramos a medir los tiempos de los trabajadores sobre las 10 próximas unidades:
- ¿Saldrían exactamente los mismos tiempos? ¿Por qué?
- ¿Por quién apostarías que sería más rápido en promedio? ¿Por qué?
- ¿Cuánto confiarías en tu apuesta? ¿Por qué?
FIN EJEMPLO 2
CONCLUSIÓN: SE BUSCA TÉCNICA “OBJETIVA” (JUSTIFICADA ESTADÍSTICAMENTE) QUE PERMITA COMPARAR Y DECIDIR SOBRE LOS PROMEDIOS DE MÁS DE 2 POBLACIONES (A PARTIR DE OBSERVAR MUESTRAS DE ELLAS).
EJERCICIO EVALUABLE 1: Ejecuta el bloque de código siguiente y contesta a las preguntas a continuación:
set.seed(123)
y1 = rnorm(10, 5, 1)
y2 = rnorm(10, 5, 1)
y3 = rnorm(10, 5, 1)
boxplot(y1, y2, y3, names=c("Mañana", "Tarde", "Noche"))
En el bloque de código anterior, se han generado 30 datos numéricos aleatorios: 10 por la mañana (y1), 10 por la tarde (y2), y 10 por la noche (y3). Se comparan las 3 muestras visualizando los 3 diagramas de caja comparativos.
- ¿De qué distribución de probabilidad viene cada una de las 3 muestras?
- CONTESTA AQUÍ
- ¿Cuál es la media de cada muestra?
| Franja | Media muestra |
|---|---|
| Mañana | ? |
| Tarde | ? |
| Noche | ? |
- ¿Cuál es la media de cada población (de la que viene cada muestra)?
| Franja | Media poblacional |
|---|---|
| Mañana | ? |
| Tarde | ? |
| Noche | ? |
- ¿Ha coincidido la media de cada muestra con la media de la población de la que proviene? ¿Cuál crees que es la razón?
- CONTESTA AQUÍ
- Si te interesa conseguir valores “generalmente más altos”, y miras solo los diagramas de caja, ¿en qué franja del día interesa obtener el nuevo dato?
- CONTESTA AQUÍ
- Si te interesa conseguir valores “generalmente más altos”, y miras solo el código R que ha generado los datos, ¿en qué franja del día interesa obtener el nuevo dato?
- CONTESTA AQUÍ
FIN EJERCICIO EVALUABLE 1
2. El modelo ANOVA de un factor
El modelo ANOVA de un factor involucra:
- Una variable de gran interés, \(Y\), numérica.
- Un factor \(X\) con valores (niveles) \(x_1\), \(x_2\), \(\ldots\), \(x_a\) (\(a\) = número de niveles).
- Unos parámetros \(\mu_1, \ldots, \mu_a\) (uno por cada nivel del factor \(X\)) y otro parámetro \(\sigma^2 > 0\)
- Cuando \(X=x_i\), entonces \(Y \sim \text{N}(\mu_i, \sigma^2)\).
- 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}\), donde para cada nivel \(x_i\) de \(X\), se han muestreado \(n_i\) valores de \(Y\), es decir, \(y_{i1}\), \(y_{i2}, \ldots,\) \(y_{in_i}\))
- Entonces \(y_{ij} = \mu_i + e_{ij}\) donde:
- \(\mu_i\) es el valor medio provocado por el nivel \(x_i\) del factor \(X\)
- \(e_{ij}\) es el error aleatorio normal de media \(0\) y varianza \(\sigma^2\), que es independiente a cada muestreo, tanto de \(x_i\) como de otros \(e_{ij}\) anteriores.
En el ejemplo de los trabajadores:
- La variable de interés es \(Y\) = “tiempo de tarea”.
- El factor que puede influir sobre el tiempo es \(X\) = “trabajador que hace la tarea”.
- Tiene 3 niveles, “A”, “B” y “C”.
- Desconocemos si el tiempo de cada trabajador sigue el modelo normal, pero podemos suponerlo, y luego ya veremos si es compatible con los resultados.
- Desconocemos si los 10 tiempos de cada trabajador son independientes entre sí, e independientes de los tiempos de los otros trabajadores, pero podemos suponerlo, y luego ya veremos si es compatible con los resultados.
3. La función aov(), que lo hace casi todo, y summary(), que te lo enseña
- Los datos de \(Y\) deben estar “colocados” en un solo vector (o columna de hoja de datos), y debe haber otro vector (o columna de hoja de datos) que indique el nivel del factor \(X\) para cada dato de la \(Y\). Si no se tienen, hay que crear el vector de datos de \(Y\) concatenando las muestras separadas de \(Y\), y el vector de niveles del factor \(X\), con repetición de los niveles tantas veces como datos de \(Y\) haya asociados a cada nivel.
EJEMPLO 3: Partiendo de los vectores A, B y C del EJEMPLO 1, construye un único vector con todos los tiempos (al que puedes llamar tiempos) y otro vector “paralelo” que contenga la etiqueta del trabajador de cada tiempo (al que puedes llamar trabajador).
# aquí el código RFIN EJEMPLO 3
- Función:
aov( formula, data,... ) - Argumentos:
data: hoja de datos con columnas (o nada, si los datos van en vectores sueltos)formula: expresión de la dependencia usando los nombres de los vectores o las etiquetas de las columnas, si van en hoja de datos:Y~X(también se puede hacer ANOVA de dos y más factores, poniendoY~X1+X2+...los puntos suspensivos indican más factores si hay).
Devuelve 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:
summary( obj,... ) - Argumentos:
obj: objeto resultante de llamar a la funciónaov().
Devuelve la típica tabla ANOVA, con los grados de libertad (
Df), las sumas de cuadrados (Sum Sq) inter (con el nombre del factor) e intra (Residuals), las medias de cuadrados (Mean Sq), el valor del estadístico de contraste (F), y su \(p\)-valor (Pr(>F)), para decidir el contraste.
4. Tareas del ANOVA
4.1. Decidir si \(Y\) es igual en todos los niveles de \(X\) o no
summary()>Pr(>F)y comparar con la significación \(\alpha\) para decidir.- Recuerda: “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.3)
- Aceptar \(H_0\) significa…
4.2. Decidir si los datos son adecuados para aplicarles el ANOVA
- Para poder aplicar ANOVA y disfrutar de sus resultados, las \(Y\)’s deben cumplir unas condiciones:
- Los residuos (errores estimados \(y_{ij} - \overline{y}_{i*}\)) deben ajustarse al modelo normal (hacer contraste de Lilliefors, o Shapiro-Wilks).
- Arrupados por nivel del factor, deberían tener varianzas similares: hacer contraste de Bartlett
- Todos independientes unos de otros e independientes del valor de la \(X\).
- Esta labor se suele realizar intuitivamente con dos figuras.
- Función
plot(...)aplicada sobre el objeto devuelto poraov(), devuelve varias figuras, entre las que :Residuals vs Factor Levels: se puede apreciar tanto si los errores tienen una “tendencia” que los hace “no independientes”, como si la varianza de es común en todos los niveles o no.- Eje X: los niveles del factor \(X\)
- Eje Y: errores estimados (tipificados)
Normal Q-Q: se relacionan los cuantiles teóricos de la normal con los de los errores estimados. Si se alejan “mucho” los puntos de la diagonal.
COMO SON DÍFÍCILES DE INTERPRETAR, SE RECOMIENDA APLICAR:
- Contraste de Bartlett de igualdad de varianzas (
help(bartlett.test)): a todos los datos juntos, ya que compara las varianzas de los niveles - Contrastes de normalidad (
help(shapiro.test)): a cada grupo de datos por separado
Se deberían aceptar todos para poder aplicar ANOVA con ciertas garantías.
4.3. Comparaciones a posteriori (por el método LSD de Fisher), sólo si no todos los niveles son iguales:
- Si la tabla ANOVA conduce a aceptar que \(Y\) no depende de \(X\), ya está todo terminado.
- En caso contrario, interesa conocer que niveles de \(X\) dan mayor media de \(Y\). Hacemos un ranking de los niveles de \(X\).
- Función
LSD.test(y, trt, alpha, console), perteneciente al paqueteagricolae, que se debe cargar (e instalar si no lo está). - Argumentos:
y: objeto devuelto por la funciónaov().trt: nombre (etiqueta entrecomillada) de la variable factor.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 del nivel.
- Puede aparece un nivel con dos códigos (cuando no se puede distinguir de otros dos niveles, que sí se distinguen entre sí).
- el ranking: código (
- Una tabla con los niveles del factor y, para cada uno,
5. Ejercicio de entrenamiento
Usa los datos de flores de R (carga con data(iris)) para ver si cada una de las 4 dimensiones medidas se pueden diferenciar por la “especie” de flor o no. Realiza todos los pasos descritos, con comentarios.
6. Más ejercicios evaluables
Carga el espacio de trabajo mt1021-1415-la-4-anova.RData, y en él encontrarás las variables que necesitas para los ejercicios.
EJERCICIO EVALUABLE 2: 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 un gráfico “sencillo” del tiempo en función del sistema operativo, que visualice la posible influencia del sistema operativo sobre el tiempo.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Calcula los tiempos medios bajo cada sistema operativo, y establece un ranking de “rapidez”.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Aplica el ANOVA para demostrar o refutar que el tiempo medio depende del sistema operativo. ¿Qué se puede decir usando un nivel de confianza del 5%?
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Comprueba si es razonable aplicar ANOVA a estos datos, o si no se debería aplicar. Figuras y comentario.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- En caso afirmativo, y si ANOVA resultó en que el sistema operativo SÍ influía en los tiempos medios, haz un nuevo ranking donde se vea si hay diferencias estadísticamente significativas (al 5%) entre los sistemas operativos.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDENFIN EJERCICIO EVALUABLE 2
EJERCICIO EVALUABLE 3: Las variables prodA, prodB, prodC y prodD recogen la producción (en miles de unidades) de 4 líneas en las que se aplican métodos distintos (A, B, C y D) que se quieren comparar.
- Realiza un gráfico “sencillo” de la producción en función del método, que visualice la posible influencia del método sobre la producción.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Calcula las producciones medias bajo cada método, y establece un ranking de “productividad”.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Aplica el ANOVA para demostrar o refutar que la productividad media depende del método. ¿Qué se puede decir usando un nivel de confianza del 1%?
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- Comprueba si es razonable aplicar ANOVA a estos datos, o si no se debería aplicar. Figuras y comentario.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS SI TAMBIÉN SE PIDEN- En caso afirmativo, y si ANOVA resultó en que el métedo SÍ influía en la productividad media, haz un nuevo ranking donde se vea si hay diferencias estadísticamente significativas (al 5%) entre los métodos.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS, QUE TAMBIÉN SE PIDENFIN EJERCICIO EVALUABLE 3
EJERCICIO EVALUABLE 4: Una asignatura tiene 5 grupos de laboratorio, con un profesor distinto en cada grupo. Alumnos de algún grupo se han quejado por creer que con su profesor, las notas son más bajas que en otros grupos. Utiliza los datos labo de las notas para saber qué puede aportar la estadistica a este caso. Describe el mismo esquema de trabajo de los ejericicios anteriores.
# Escribe aquí el codigo R de lo que se pide Y TUS COMENTARIOS, QUE TAMBIÉN SE PIDENFIN EJERCICIO EVALUABLE 4