## Warning: package 'knitr' was built under R version 4.1.3

PASOS PARA RESOLVER ESTOS EJERCICIOS. En general, se trata de analizar si la medida de centralización (media o mediana) de cierta variable cuantitativa se comporta de la misma manera en ciertos grupos.

OJO: en las soluciones he ordenado las medias/medianas de distintas formas (diagrama, con desigualdades, explicándolo de palabra…). NO es necesario hacerlo de todas las formas posibles, aunque tienes que entenderlas todas (porque se usan indistintamente unas u otras).

Ejercicio 1

Considera el ya famoso fichero robles. Contrasta si la concentración media de Magnesio es la misma en las zonas muestreadas.

  1. Descarga y lee el fichero robles. Puedes usar el botón Import Datase, o bien la orden

    datos = 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           1

    y 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    1

  2. Decide si hay al menos una zona en la que la concentración media de Magnesio es diferente del resto al nivel de siginificación \(\alpha = 0.05\). 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.367

    Se 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, se suele usar este valor como referencia independientemente del valor que se use para el contraste ANOVA), 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 ' ' 1

    que 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
resumen = summary(modelo) 
resumen$adj.r.squared
    [1] 0.08303731

    que resulta ser muy bajo, indicando que el modelo tiene poco poder explicativo.

  3. En caso afirmativo, ordena las zonas de acuerdo con la concentración media de Magnesio (puedes usar Bonferroni o Tukey). Ahora nos interesa ordenar las medias poblacionales. 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 de ellos. 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 

  1       2       3      
2 0.02569 -       -      
3 0.61144 1.00000 -      
4 0.05561 0.00042 0.01943

P value adjustment method: bonferroni

    Se puede usar los boxplots de la exploración inicial, o bien

    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

    para calcular las medias muestrales y tener una idea de cómo podrían estar ordenadas las medias poblacionales. Al nivel de significación \(\alpha = 0.05\) se tiene (aqui tienes el diagrama de medias ordenadas)

    • No hay evidencias muestrales para rechazar que las concentraciones medias de Magnesio en las zonas 2 y 3 cumplen \(\mu_2=\mu_3\), que son las medias mayores.
    • No hay evidencias muestrales para rechazar \(\mu_1=\mu_3\) pero \(\mu_1<\mu 2\).
    • No hay evidencias muestrales para rechazar \(\mu_4=\mu_1\) pero \(\mu_4<\mu_2\) y \(\mu_4<\mu_3\).\[ \] Al nivel de significación \(\alpha = 0.01\) se tiene
    • No hay evidencias muestrales para rechazar \(\mu_1=\mu_2=\mu_3\).
    • No hay evidencias muestrales para rechazar \(\mu_1=\mu_3 = \mu_4\).
    • Pero \(\mu_2>\mu_4\). \[ \] Si trabajamos con el método de Tuckey hay que cuidar un detalle técnico; convierte la segunda columna, que es numérica (números enteros), en un factor
        datos[ , 1] = as.factor(datos[, 1])

    y ahora, al nivel de significación \(\alpha = 0.05\)

    # 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,10,1))
plot(cld(glh, decreasing = FALSE, level = 0.05), 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 estadísticamente significativa -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
    • No hay diferencias significativas entre las concetraciones medias de Magnesio en las zonas 2 y 3.
    • Tampoco hay diferencias entre las concentraciones medias en las zonas 1 y 3, pero sí entre las zonas 1 y 2.
    • Tampoco entre las concentraciones medias en las zonas 1 y 4, pero sí entre las zonas 3 y 4.
    • Esto se puede simbolizar como \[ X1,\,X4 \leq X1,\, X3 \leq X2,\, X3\] usando ahora el nivel de significación \(\alpha = 0.01\)
    # 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,10,1))
plot(cld(glh, decreasing = FALSE, level = 0.01), col = "tan")

    • Ahora no hay diferencias significativas en la concentración media de Magnesio en las zonas 1, 2 y 3.
    • No hay diferencias significativas entre las zonas 1 y 4.
    • La concetración media en las zonas 2 y 3 es mayor que en la zona 4.
    • Esto se puede resumir con \[X1,\, X4, \leq X1,\, X2,\, X3\]

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:

  1. Descarga y lee el fichero de datos. Puedes usar el botón Import Datase, o bien la orden

    datos = read.table(file = "Practica10-NOx-30-dic-2020_limpio.csv", sep = ";", header = T)

    y asegurarnos de que han sido bien leídos

    head(datos)
               Fecha Getafe Leganes Alcala Alcobendas Fuenlabrada Mostoles Torrejon
1 29/12/20 17:00     33      42     32         15          24       17       25
2 29/12/20 18:00     55      61     58         25          53       33       34
3 29/12/20 19:00     73      42     43         26          41       22       47
4 29/12/20 20:00     45      33     80         23          30       17       63
5 29/12/20 21:00     28      31    138         15          21       18       47
6 29/12/20 22:00     23      29    107         12          16       13       28
  Alcorcon Coslada
1       22      34
2       42      49
3       28      89
4       24      42
5       20      47
6       18      27

    Nos quedamos con las columnas 2, 3, 4, 5 y 6 (la primera contiene la fecha y hora de la recogida de datos) ya que se pide las primeras 5 ciudades de la tabla; y hacemos una pequeña exploración inicial para asegurar la integridad de los datos:

    datos = datos[ , 2:6]
summary(datos)
    ##      Getafe         Leganes          Alcala         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.0

  2. Analiza 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 grá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}\]

    ###--- 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%)

    # 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 |     Alcala   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
         |
 Leganes |   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

    aggregate(datos$Respuesta ~ datos$Tratamiento, FUN = median)
      datos$Tratamiento datos$Respuesta
1            Getafe              27
2           Leganes              29
3            Alcala              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 (ojo: con Dunnet se compara el p-valor de la comparación de cada par de medianas con \(\alpha/2\)). Observa que (aqui está el diagrama con las medianas ordenadas)

    • Alcalá tiene la mediana muestral más grande y todas las comparaciones 2 a 2 son significativas (mira los p-valores de la primera columna de la tabla de Dunnet), por lo que esa diferencia se traslada a nivel poblacional: podemos suponer que la concentración media de SOx es mayor en la estación de Alcalá de Henares que en el resto de estaciones.
    • 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 \[ \] Ahora, al nivel de significación es del 5%
    # 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 = .05)
      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 |     Alcala   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
         |
 Leganes |   2.692975  -1.707022  -1.489230  -0.026488
         |    0.0035*     0.0439     0.0682     0.4894

alpha = 0.05
Reject Ho if p <= alpha/2

    Observa que ahora

    • Alcalá sigue teniendo la concentración mediana (poblacional) mayor que el resto.
    • No hay motivos para considerar que las concetraciones medianas del resto de estaciones sean diferentes entre sí.
    • Se podría resumir esta información con el código de letras que ya conoces:
      • Alcalá B
      • Leganés A, Getafe A, Fuenlabrada A, Alcobendas A