Práctica 9. ANOVA. Soluciones
Universidad de Alcalá.
Estadística (650008). Grado en Biología sanitaria.
Actualizado: 2021-01-16
Ejercicio 1
Resuelve el ejemplo ANOVA, cortesía de la unidad de Bioestadística Clínica del hospital Ramón y Cajal.
En concreto:
Crea un fichero csv con los datos de los tratamientos: aquí tienes la tabla de datos.
Decide si hay algún tratamiento que produce una respuesta distinta de los demás. Lo primero es leer los datos
datos = read.table(file = "Practica09-RamonCajal.csv", sep = ";", header = T)y asegurarnos de que han sido bien leídos
head(datos)## X1 X2 X3 X4 X5 ## 1 180 172 163 158 147 ## 2 173 158 170 146 152 ## 3 175 167 158 160 143 ## 4 182 160 162 171 155 ## 5 181 175 170 155 160a partir de aquí, vamos a ir usando distintos fragmentos de la plantilla, y comentando los resultados obtenidos:
####--- poner datos en formato correcto datos = stack(datos) datos = datos[, c(2,1)] colnames(datos) <- c("Tratamiento", "Respuesta") head(datos)## Tratamiento Respuesta ## 1 X1 180 ## 2 X1 173 ## 3 X1 175 ## 4 X1 182 ## 5 X1 181 ## 6 X2 172###--- EXPLORACION # boxplots bp = boxplot(Respuesta ~ Tratamiento, data=datos, col="tan") stripchart(Respuesta ~ Tratamiento, data=datos, col="red", vertical = TRUE, method = "jitter", cex=0.8, add=TRUE, pch=19)
Observa que hay pocos datos en cada grupo, lo que permite que los boxplots puedan tener formas más o menos “raras”. Nosotros asumimos la hipótesis independencia en la toma de las muestras, y tenemos que comprobar la normalidad y la homocedasticidad de la muestra. Podemos usar métodos gráficos
###--- Graficos diagnostico condiciones ANOVA modelo = lm(datos$Respuesta ~ datos$Tratamiento) par(mfrow = c(1, 2), mar = c(5,5,2,1)) for(i in 1:2){ plot(modelo, which = i) }
par(mfrow = c(1, 1))El qqplot (panel derecho) no muestra desviaciones respecto de la normalidad. En cuanto a la homocedasticidad, fíjate (panel izquierdo) en que el tratamiento que tiene una media mayor (algo más de 175) parece tener una dispersión menor que el resto.
Podemos usar los contrastes de hipótesis para salir de dudas. A continuación, se incluye todo el bloque de contrastes que hay en el script: se hace uno de normalidad (aunque el método gráfico no hace sospechar nada) y otro de homocedasticidad (ahí, el método grafico puede generar alguna duda)
###--- contrastes normalidad # test Shapiro-Wilks # H0: las muestras provienen de poblaciones normales shapiro.test(modelo$residuals)## ## Shapiro-Wilk normality test ## ## data: modelo$residuals ## W = 0.98911, p-value = 0.9927# test Kolmogorov-Smirnov-Lilliefors ## H0: las muestras provienen de poblaciones normales # library(nortest) # lillie.test(modelo$residuals) ###--- contrastes homocedasticidad ## H0: varianzas iguales (caso de poblaciones normales) bartlett.test(Respuesta ~ Tratamiento, data = datos)## ## Bartlett test of homogeneity of variances ## ## data: Respuesta by Tratamiento ## Bartlett's K-squared = 2.6811, df = 4, p-value = 0.6125## H0: varianzas iguales (caso de poblaciones cualesquiera) # library(car) # leveneTest(Respuesta ~ Tratamiento, data = datos)En todo caso, ambos p-valores son grandes, por lo que no hay porqué dudar de que se cumplan las hipótesis.
Procedemos a hacer el contraste ANOVA:
###--- TABLA ANOVA # H0: las medias son iguales # H1: al menos una media es diferente modelo = lm(datos$Respuesta ~ datos$Tratamiento) # p-valor de contraste anova(modelo)## Analysis of Variance Table ## ## Response: datos$Respuesta ## Df Sum Sq Mean Sq F value Pr(>F) ## datos$Tratamiento 4 2010.6 502.66 11.24 6.062e-05 *** ## Residuals 20 894.4 44.72 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1que resulta ser significativo: al menos un tratamiento produce una respuesta diferente de la de los demás. Podemos calcular el coeficiente de determinación \(r^2\)
# coeficiente de determinacion r2 summary(modelo)## ## Call: ## lm(formula = datos$Respuesta ~ datos$Tratamiento) ## ## Residuals: ## Min 1Q Median 3Q Max ## -12.0 -4.4 0.6 3.8 13.0 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 178.200 2.991 59.586 < 2e-16 *** ## datos$TratamientoX2 -11.800 4.229 -2.790 0.011305 * ## datos$TratamientoX3 -13.600 4.229 -3.216 0.004339 ** ## datos$TratamientoX4 -20.200 4.229 -4.776 0.000115 *** ## datos$TratamientoX5 -26.800 4.229 -6.337 3.48e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.687 on 20 degrees of freedom ## Multiple R-squared: 0.6921, Adjusted R-squared: 0.6305 ## F-statistic: 11.24 on 4 and 20 DF, p-value: 6.062e-05que es Multiple R-squared: 0.6921; bastante bueno; el modelo tiene bastante poder explicativo.
En caso afirmativo, ordena los tratamientos de acuerdo con su efectividad para reducir la hipertensión (puedes usar Bonferroni o Tukey). Ahora nos interesa ordenar las medias poblacionales, es decir, la respuesta que producen en la presión arterial los distintos tratamientos. Esencialmente, se hacen comparaciones de medias 2 a 2, pero teniendo la precaución de repartir el error de tipo I entre todas las comparaciones para que el error de tipo I acumulado no se dispare. Hay distintas formas de hacer esas comparaciones, nosotros analizamos dos. Podemos usar el ajuste de Bonferroni
# Ajuste Bonferroni (ptt=pairwise.t.test(datos$Respuesta, datos$Tratamiento, p.adj="bonferroni", pool.sd=FALSE))## ## Pairwise comparisons using t tests with non-pooled SD ## ## data: datos$Respuesta and datos$Tratamiento ## ## X1 X2 X3 X4 ## X2 0.1914 - - - ## X3 0.0210 1.0000 - - ## X4 0.0472 1.0000 1.0000 - ## X5 0.0016 0.0982 0.0906 1.0000 ## ## P value adjustment method: bonferronique produce una salida un poco complicada de resumir (puedes ver un ejemplo de cómo hacerlo en el ejercicio 2 para el caso de comparaciones múltiples de medianas con el test de Dunnet), o el de Tuckey
# Ajuste Tukey, representaciones graficas # install.packages("multcomp") # por si no lo tienes instalado library(multcomp) datos.aov = aov(Respuesta ~ Tratamiento, data = datos) glh = summary(glht(model=datos.aov, linfct = mcp(Tratamiento = "Tukey"))) par(mar=c(4,4,8,1)) plot(cld(glh, decreasing = FALSE), col = "tan")
- Los grupos que comparten letra incluyen tratamientos entre los que no hay una respuesta significativamente diferente a nivel poblacional (ojo, que las medias muestrales sí son distintas, pero la diferencia no es relevante).
- X1 y X2 forman un grupo, X2, X3 y X4 otro grupo, y X4 y X5 el último grupo.
- Efectivamente, hay tratamientos que están en 2 grupos, y eso puede parecer raro. Por ejemplo, X2:
- Se puede considerar que la respuesta de X2 no es significativamente diferente de la de X1 y, a la vez, no es significativamente distinta de la X3 y X4.
- Sin embargo, X1 no comparte letra con X3 y X4, luego entre X1 y X3, X4 sí hay diferencias significativas.
- Podemos pensar X2 como un tratamiento con respuesta intermedia; no muy diferente de X1, no muy diferente de X3 y X4, pero esa “poca” diferencia acumulada hace que entre X1 y X3, X4 sí haya diferencias significativas.
- Esto se puede simbolizar como \[ X4,\,X5 \leq X2,\, X3,\, X4 \leq X1,\, X2\]
- Recuerda que eran tratamientos para reducir la hipertensión, por tanto:
- Hay 2 mejores tratamientos: X4 y X5.
- Si un paciente no pudiera tomar X5 (por cualquier motivo), entonces tomaría X4. Pero ese tratamiento no produce respuesta significativamente distinta de X2 y X3, por lo que podría administrasele también cualquiera de estos dos.
Ejercicio 2
Este fichero de datos fue extraído de aquí el 16 de diciembre de 2020 a las 9:15, y contiene los datos de concentración de NOx medidos en 5 estaciones de la Comunidad de Madrid. En esta página hay datos de contaminantes atmosféricos en tiempo real D.G. del Medio Ambiente Area de Calidad Atmosférica - Red de Calidad del Aire de la comunidad de Madrid.
Se pide:
Analiza si la concentración del contaminante es análoga en todos ellos.
En caso contrario, ordena los municipios de más contaminado a menos contaminado (puedes usar Bonferroni o Tukey si se cumplen las condiciones para ANOVA, o el test de Dunnet en caso contrario). Usa un nivel de significación del 10%.
datos = read.table(file = "Practica09-NOx-16-dic.csv", sep = ";", header = T)y asegurarnos de que han sido bien leídos
head(datos)## Fecha Getafe Leganés Alcalá.de.Henares Alcobendas Fuenlabrada ## 1 15/12/20 09:00 44 38 29 50 21 ## 2 15/12/20 10:00 51 *** N 33 32 17 ## 3 15/12/20 11:00 32 19 24 31 13 ## 4 15/12/20 12:00 18 19 21 25 15 ## 5 15/12/20 13:00 18 21 13 23 17 ## 6 15/12/20 14:00 21 21 15 15 15 ## Móstoles Torrejón.de.Ardoz Alcorcón Coslada ## 1 32 61 32 60 ## 2 17 71 21 45 ## 3 12 72 12 55 ## 4 13 46 14 24 ## 5 14 34 16 24 ## 6 10 33 12 31Nos quedamos con las columnas 2, 3, 4, 5 y 6; y hacemos una pequeña exploración inicial para asegurar la integridad de los datos:
datos = datos[ , 2:6] summary(datos)## Getafe Leganés Alcalá.de.Henares Alcobendas ## Min. :18.00 Length:24 Min. : 13.00 Min. : 12.00 ## 1st Qu.:28.75 Class :character 1st Qu.: 20.50 1st Qu.: 28.75 ## Median :47.50 Mode :character Median : 32.00 Median : 43.50 ## Mean :49.62 Mean : 44.21 Mean : 50.92 ## 3rd Qu.:66.75 3rd Qu.: 41.50 3rd Qu.: 64.50 ## Max. :97.00 Max. :295.00 Max. :130.00 ## Fuenlabrada ## Min. : 11.00 ## 1st Qu.: 15.00 ## Median : 24.00 ## Mean : 33.08 ## 3rd Qu.: 34.50 ## Max. :103.00Debería llamarte la atención que la segunda columna contienen datos de tipo
character, pero cabría esperar que se tratara de números. ¿Qué ha podido pasar? Cabe la posibilidad de que haya alguna celda vacía, o algún caracter extraño. Como la tabla no es grande, podemos visualizar la completadatos## Getafe Leganés Alcalá.de.Henares Alcobendas Fuenlabrada ## 1 44 38 29 50 21 ## 2 51 *** N 33 32 17 ## 3 32 19 24 31 13 ## 4 18 19 21 25 15 ## 5 18 21 13 23 17 ## 6 21 21 15 15 15 ## 7 28 21 16 12 15 ## 8 25 24 15 20 19 ## 9 40 37 19 23 33 ## 10 65 54 31 35 49 ## 11 72 77 37 43 59 ## 12 81 70 40 53 36 ## 13 97 68 45 58 28 ## 14 76 67 45 64 27 ## 15 57 49 41 66 32 ## 16 51 56 21 102 30 ## 17 54 42 19 130 13 ## 18 43 27 26 106 13 ## 19 28 25 43 69 11 ## 20 29 27 37 44 21 ## 21 34 48 41 30 34 ## 22 52 66 48 41 76 ## 23 92 79 107 54 103 ## 24 83 89 295 96 97Seguro que has dertectado que aparece un
***N; parece que a esa hora/estación no se registró correctamente la medida, y se usa ese código para indicarlo. Como es una palabra, y el contenido de las columnas de un data.frame debe ser del mismo tipo, el programa ha convertido los valores de esas columna en palabras. Tendremos que eliminar esa fila (la 2)datos = datos[-c(2), ]poner datos en formato correcto
datos = stack(datos) datos = datos[, c(2,1)] colnames(datos) <- c("Tratamiento", "Respuesta")y convertir las palabras en números
datos$Respuesta = as.numeric(datos$Respuesta)A partir de aquí, seguiremos el script
##--- EXPLORACION # boxplots bp = boxplot(Respuesta ~ Tratamiento, data=datos, col="tan") stripchart(Respuesta ~ Tratamiento, data=datos, col="red", vertical = TRUE, method = "jitter", cex=0.8, add=TRUE, pch=19)
Los boxplots sugieren que no se va a cumplir ni la normalidad (bigote inferior mucho más corto que el superior) ni la homocedasticidad (boxplots de anchuras muy diferentes). Se puede corroborar mediante las herramientas gáficas y numéricas (contrastes) que conoces.
###--- Graficos diagnostico condiciones ANOVA modelo = lm(datos$Respuesta ~ datos$Tratamiento) par(mfrow = c(1, 2), mar = c(5,5,2,1)) for(i in 1:2){ plot(modelo, which = i) }
par(mfrow = c(1, 1))Llama la atención un valor muy atípico en Alcalá de Henares que, tal vez, cabría eliminar. De momento, lo mantenemos.
En todo caso, ni las columnas el gráfico de la izquierda tienen alturas similares (no hay homocedasticidad) ni la mayor parte de los puntos se disponen sobre la recta en el qqplot (gráfico de la derecha, no hay normalidad)
###--- contrastes normalidad # test Shapiro-Wilks # H0: las muestras provienen de poblaciones normales shapiro.test(modelo$residuals)## ## Shapiro-Wilk normality test ## ## data: modelo$residuals ## W = 0.70991, p-value = 9.104e-14# test Kolmogorov-Smirnov-Lilliefors ## H0: las muestras provienen de poblaciones normales # library(nortest) # lillie.test(modelo$residuals) ###--- contrastes homocedasticidad ## H0: varianzas iguales (caso de poblaciones normales) bartlett.test(Respuesta ~ Tratamiento, data = datos)## ## Bartlett test of homogeneity of variances ## ## data: Respuesta by Tratamiento ## Bartlett's K-squared = 30.102, df = 4, p-value = 4.666e-06## H0: varianzas iguales (caso de poblaciones cualesquiera) # library(car) # leveneTest(Respuesta ~ Tratamiento, data = datos)Los p-valores son muy pequeños, por lo que hacer un ANOVA no tiene sentido. Es importante resaltar que, aunque la teoría no respalde el uso de ANOVA, matemáticamente es posible hacer los cálculos. Es decir, si se lo pides, el ordenador te hará el contraste ANOVA (tú mandas) aunque del resultado que arroje no se puedan extraer conclusiones.
Optamos, pues, por hacer un contraste de Kruskal-Wallis. Recuerda que ahora la hipótesis nula es
Todas las medianas son igualesfrente a la hopótesis alternativaAlguna de las medianas es diferente de las demás###--- NO PARAMETRICOS: si falla la normalidad # KRUSKAL-WALLIS equivale al ANOVA # H0: medianas iguales # H1: al menos una mediana diferente kruskal.test(datos$Respuesta, datos$Tratamiento)## ## Kruskal-Wallis rank sum test ## ## data: datos$Respuesta and datos$Tratamiento ## Kruskal-Wallis chi-squared = 11.871, df = 4, p-value = 0.01834el p-valor no es excesivamente pequeño: al 1% de significación no se rechazaría H0, pero al 5% y al 10% (que es el que indica el enunciado) sí, por lo que se decide que rechazar H0. Para ordenar las medianas y determinar qué estaciones registran concentraciones de NOx significativamente diferentes (superiores/inferiores) usaremos el test de Dunnet (usamos un nivel de significación es del 10%)
# Dunnett es equivalente a ordenar las medias # install.packages("dunn.test") # por si no lo tienes instalado library(dunn.test) dunn.test(datos$Respuesta, datos$Tratamiento, alpha = .1)## Kruskal-Wallis rank sum test ## ## data: x and group ## Kruskal-Wallis chi-squared = 11.8712, df = 4, p-value = 0.02 ## ## ## Comparison of x by group ## (No adjustment) ## Col Mean-| ## Row Mean | Alcalá.d Alcobend Fuenlabr Getafe ## ---------+-------------------------------------------- ## Alcobend | -1.864647 ## | 0.0311* ## | ## Fuenlabr | 0.721085 2.585733 ## | 0.2354 0.0049* ## | ## Getafe | -1.999574 -0.134927 -2.720660 ## | 0.0228* 0.4463 0.0033* ## | ## Leganés | -1.501892 0.362754 -2.222978 0.497681 ## | 0.0666 0.3584 0.0131* 0.3094 ## ## alpha = 0.1 ## Reject Ho if p <= alpha/2junto con el valor muestral de la mediana de cada problación
aggregate(datos$Respuesta ~ datos$Tratamiento, FUN = median)## datos$Tratamiento datos$Respuesta ## 1 Getafe 44 ## 2 Leganés 42 ## 3 Alcalá.de.Henares 31 ## 4 Alcobendas 44 ## 5 Fuenlabrada 27Observa que desde un punto de vista muestral, de mayor a menor, las medianas se ordenan como sigue: Alcalá, Leganés, Getafe, Alcobendas y Fuenlabrada. Sin embargo, a nivel poblacional, no todas las diferencias son significativas, porque en la tabla de Dunnet no todos los p-valores tienen asignado un asterisco. Observa que, referido a la mediana,
- Getafe y Alcobendas tienen la mediana muestral más grande (coinciden). Entre las medianas de Getafe y Alcobendas no hay diferencias significativas (a nivel poblacional), como muestra el p-valor de la tercera fila y seunda columna.
- Getafe presenta diferencias significativas con Alcalá y Fuenlabrada, pero no con Alcobendas (fila 3 de la tabla de Dunnet).
- Lo mismo le sucede a Alcobendas (entrada 1 y 1 de la tabla y columna 2).
- Leganés prese nta diferencias significativas con Fuenlabrada, pero no con Alcalá (fila 4 de la tabla).
- Y Alcalá no presenta diferencias significativas con Fuenlabrada.
- Se podría resumir esta información con el código de letras que ya conoces:
- Getafe C, Alcobandas C, Leganés C, B
- Alcalá B, A
- Fuenlabrada A