Práctica 8. ANOVA. Soluciones
Universidad de Alcalá.
Estadística. Grado en Ciencias Ambientales.
Actualizado: 2021-12-20
Ejercicio 1
Considera el ya famoso fichero robles. Contrasta si la concentración media de Magnesio es la misma en las zonas muestreadas.
Descarga y lee el fichero robles. Puedes usar el botón
Import Datase, o bien la ordendatos = read.table(file = "Practica06-robles.csv", sep = ";", header = T, dec = ".")y asegurarnos de que han sido bien leídos
head(datos)## Num Hierro Manganeso Zinc Calcio Magnesio Potasio Fosforo Nitrogeno Zona ## 1 1 0.058 0.0303 0.0089 2.365 0.400 2.632 0.145 2.776 1 ## 2 2 0.060 0.0294 0.0109 2.745 0.432 2.495 0.161 2.918 1 ## 3 3 0.058 0.0289 0.0090 2.513 0.349 2.396 0.169 4.826 1 ## 4 4 0.059 0.0275 0.0090 2.361 0.349 1.979 0.155 4.893 1 ## 5 5 0.061 0.0283 0.0087 2.330 0.350 2.197 0.139 3.386 1 ## 6 6 0.027 0.0110 0.0097 2.450 0.367 2.233 0.148 3.776 1 ## Variedad Tratamiento ## 1 A 0 ## 2 A 0 ## 3 A 0 ## 4 A 0 ## 5 A 0 ## 6 B 1y nos quedamos con las columnas que indica el enunciado
datos = datos[, c("Magnesio", "Zona")] head(datos)## Magnesio Zona ## 1 0.400 1 ## 2 0.432 1 ## 3 0.349 1 ## 4 0.349 1 ## 5 0.350 1 ## 6 0.367 1Decide si hay al menos una zona en la que la concentración media de Magnesio es diferente del resto. A partir de aquí, seguiremos la plantilla ANOVA, e iremos comentando los resultados obtenidos. Los datos ya están en el formato correcto, es decir, una columna con la variable respuesta y otra con la explicativa; pondremos la explicativa primero y la respuesta después
datos = datos[, c(2,1)]y asignamos a las columnas los nombres que se usan en el script
colnames(datos) <- c("Tratamiento", "Respuesta") head(datos)Tratamiento Respuesta 1 1 0.400 2 1 0.432 3 1 0.349 4 1 0.349 5 1 0.350 6 1 0.367Se exploran los datas
#--- 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)
No hay muchos 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) muestra desviaciones respecto de la normalidad, y el izquierdo la dispersión de lso residuos (hace referencia a la homocedasticidad). Los gráficos no son del todo inequívocos, de modo que haremos los contrastes de hipótesis correspondientes para salir de dudas. A continuación, se incluye todo el bloque de contrastes que hay en el script: se hace uno de normalidad y otro de homocedasticidad
###--- 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.93733, p-value = 0.03435# 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 = 4.7222, df = 3, p-value = 0.1933## H0: varianzas iguales (caso de poblaciones cualesquiera) # library(car) # leveneTest(Respuesta ~ Tratamiento, data = datos)En todo caso, ambos p-valores son suficientemente grandes (mayores que 0.01), por lo que no hay porqué dudar de que se cumplan las hipótesis. Procedemos a hacer el contraste ANOVA: \[H_0:\,\text{ Las medias son iguales en las cuatro zonas}\] \[H_1:\,\text{ Al menos una media es distinta en alguna de las cuatro zonas}\]
###--- 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 1 0.022917 0.0229168 4.3506 0.04414 * Residuals 36 0.189630 0.0052675 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1que resulta ser significativo: al menos en una zona hay una media (poblacional) 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 -0.097361 -0.063752 -0.009666 0.047259 0.159595 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.45045 0.02811 16.025 <2e-16 *** datos$Tratamiento -0.02152 0.01032 -2.086 0.0441 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.07258 on 36 degrees of freedom Multiple R-squared: 0.1078, Adjusted R-squared: 0.08304 F-statistic: 4.351 on 1 and 36 DF, p-value: 0.04414que 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))que produce una salida que puede resultar 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. Un detalle técnico; convierte la segunda columna, que es numérica (números enteros), en un factor
datos[ , 1] = as.factor(datos[, 1])Ahora
# Ajuste Tukey, representaciones graficas # install.packages("multcomp") # por si no lo tienes instalado library(multcomp)Warning: package 'multcomp' was built under R version 4.1.2Warning: package 'TH.data' was built under R version 4.1.2datos.aov = aov(Respuesta ~ Tratamiento, data = datos) glh = summary(glht(model=datos.aov, linfct = mcp(Tratamiento = "Tukey"))) par(mar=c(4,4,10,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). Podemos calcular las medias de cada muestra
aggregate(datos$Respuesta ~ datos$Tratamiento, FUN = mean)datos$Tratamiento datos$Respuesta 1 1 0.381400 2 2 0.461400 3 3 0.429125 4 4 0.323300- Entre las medias poblacionales de los niveles 1 y 4 no hay diferencia significativa.
- Entre las medias poblacionales de los niveles 1 y 3 no hay diferencia significativa.
- Pero entre las medias poblacionales de los niveles 1 y 3 sí hay diferencia significativa (Por eso no comparten las dos letras).
- De la misma manera, entre los niveles 2 y 3 no hay diferencias significativas, pero sí las hay entre los niveles 2 y 3.
- Esto se puede simbolizar como \[ X1,\,X4 \leq X1,\, X3, \leq X2,\, X3\]
- 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 30 de diciembre de 2020 a las 16:02, 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:
Descarga y lee el fichero de datos. Puedes usar el botón
Import Datase, o bien la ordendatos = read.table(file = "Practica08-NOx-30-dic-2020_limpio.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 29/12/20 17:00 33 42 32 15 24 2 29/12/20 18:00 55 61 58 25 53 3 29/12/20 19:00 73 42 43 26 41 4 29/12/20 20:00 45 33 80 23 30 5 29/12/20 21:00 28 31 138 15 21 6 29/12/20 22:00 23 29 107 12 16 Móstoles Torrejón.de.Ardoz Alcorcón Coslada 1 17 25 22 34 2 33 34 42 49 3 22 47 28 89 4 17 63 24 42 5 18 47 20 47 6 13 28 18 27Nos quedamos con las columnas 2, 3, 4, 5 y 6 (la primera contiene la fecha y hora de la recogida de datos); 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. : 6.00 Min. : 6.00 Min. : 15.00 Min. : 5.00 ## 1st Qu.:10.50 1st Qu.:10.00 1st Qu.: 33.00 1st Qu.: 9.00 ## Median :27.00 Median :29.00 Median : 44.00 Median :12.00 ## Mean :29.05 Mean :31.58 Mean : 93.16 Mean :19.84 ## 3rd Qu.:46.50 3rd Qu.:42.00 3rd Qu.:122.50 3rd Qu.:22.00 ## Max. :73.00 Max. :83.00 Max. :321.00 Max. :82.00 ## Fuenlabrada ## Min. : 3.0 ## 1st Qu.: 5.5 ## Median :18.0 ## Mean :20.0 ## 3rd Qu.:30.0 ## Max. :53.0Analiza si la concentración del contaminante es análoga en todos ellos. Empezamos por poner datos en formato correcto (una columna con la variable respuesta y otra con la explicativa)
datos = stack(datos) datos = datos[, c(2,1)] colnames(datos) <- c("Tratamiento", "Respuesta")A partir de aquí, seguiremos el script para explorar los datos
##--- 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))Efectivamente, 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.8008, p-value = 5.312e-10# 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 = 88.714, df = 4, p-value < 2.2e-16## 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 son \[H_0:\,\text{ Las medianas son iguales en las cuatro zonas}\] \[H_1:\,\text{ Al menos una mediana es distinta en alguna de las cuatro zonas}\]
```r
###--- 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 = 24.657, df = 4, p-value = 5.897e-05
```
el p-valor es muy pequeño, por lo que hay 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%)
```r
# 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 = 24.6568, df = 4, p-value = 0
Comparison of x by group
(No adjustment)
Col Mean-|
Row Mean | Alcalá.d Alcobend Fuenlabr Getafe
---------+--------------------------------------------
Alcobend | 4.399997
| 0.0000*
|
Fuenlabr | 4.182205 -0.217792
| 0.0000* 0.4138
|
Getafe | 2.719463 -1.680534 -1.462741
| 0.0033* 0.0464* 0.0718
|
Leganés | 2.692975 -1.707022 -1.489230 -0.026488
| 0.0035* 0.0439* 0.0682 0.4894
alpha = 0.1
Reject Ho if p <= alpha/2
```
junto con el valor muestral de la mediana de cada problación
```r
aggregate(datos$Respuesta ~ datos$Tratamiento, FUN = median)
```
```
datos$Tratamiento datos$Respuesta
1 Getafe 27
2 Leganés 29
3 Alcalá.de.Henares 44
4 Alcobendas 12
5 Fuenlabrada 18
```
Observa 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
* Alcalá tiene la mediana muestral más grande y todas las comparacioens 2 a 2 son significativas (mira los p-valores de la primera columna de la tabla de Dunnet), por qlo que esa diferencia se traslada a nivel poblacional: podemos suponer que la concentració n de SOx es mayor en Alcalá de Henares.
* La siguiente mediana más grande (muestralmente) es la de Leganés. Si atendemos a la última fila de la tabla de Dunnet, se observa que hay diferencias significativas con todas las demás ciudades excepto con Getafe. Por tanto, la concentración de NOx en Leganés y Getafe es similar (no es diferente), pero es significativamente menor en ambas que en Alcalá.
* La siguiente mediana más grande (muestralmente) es la de Getafe (penúltima fila de la tabla): su mediana no es diferente (poblacionalmente hablando) que la de Fuenlabrada, pero sí mayor que la de Alcobendas. Recuerda que la de Alcobendas si es menor que la de Leganés.
* Finalmente, la segunda fila de la tabla indica que no hay diferencias significativas entre Fuenlabrada y Alcobendas.
* Se podría resumir esta información con el código de letras que ya conoces:
* Alcalá D
* Leganés C, Getafe C, B
* Fuenlabrada B, A, Alcobendas A