TEMA 6. RELACIÓN ENTRE UNA VARIABLE CUALITATIVA Y OTRA CUANTITATIVA

INTRODUCCION

El análisis de este tipo de asociación implica comparar la distribución de una variable para los distintos valores que toma otra u otras. Normalmente, se considera a la variable cualitativa como independiente o condicionante, mientras que la variable cuantitativa es la dependiente o condicionada.

Supongamos que, tras medir el peso de 100 individuos (50 mujeres -grupo A- y 50 hombres -grupo B-) queremos saber si existe relación entre peso y género. Esta comparación se puede realizar tanto para los valores medios como para su dispersión. Si el peso de varones y mujeres es el mismo, podríamos afirmar que el peso no dependería del género. En este caso podríamos plantear las siguientes hipótesis:

⚠️ ** RECUERDA**:

• Si el valor-p α>0,05 se acepta la Hipótesis nula (H0): el peso es independiente del género (no hay relación entre ellas).

• Si el valor-p α<0,05 se rechaza la hipótesis nula y se acepta Hipótesis alternativa (H1): las dos variables están relacionadas (no son independientes), pues el peso depende del género.

Para elegir la prueba que permita realizar este contraste, debemos tener en cuenta los siguientes aspectos:

• Si las variables siguen una distribución normal o no.

• Si ambas variables provienen de:

  • Muestras independientes a aquellas en las que los individuos en una y otra son distintos. Este es el caso anteriormente mencionado del peso de los 100 individuos, divididos en dos grupos que no relacionados (es decir, son independientes o no pareados).

  • Se considera muestras dependientes a aquellas que se han tomado sobre los mismos individuos, pero en diferentes circunstancias(pe. en diferentes momentos del tiempo, tras algún fenómeno, por ejemplo, antes y después del inicio de un tratamiento, o de la promulgación de una norma que regula el tráfico en una barrio). En este caso, se considera que son dependentes o pareadas.

• El número de grupos.

COMPARACIÓN DEL PROMEDIO DE DOS GRUPOS

Está basado en la prueba t, un procedimiento estadístico que comprueba si existe una diferencia significativa entre las medias de dos grupos. Hay tres tipos diferentes de pruebas t:

• La prueba t de una muestra.

• La prueba t para muestras independientes.

• La prueba t para muestras relacionadas.

Fuente: DataLab (https://datatab.es/tutorial/t-test)

Comparación de la media de una única muestra (Prueba t de una muestra o “one sample” test)

Utilizamos la prueba t de una muestra cuando queremos comparar la media de una muestra con una media de referencia conocida.

Fuente: DataLab (https://datatab.es/tutorial/t-test)

Supongamos que queremos saber si el peso medio de los alumnos de la UC encuestados es diferente del peso medio recomendado por la OMM que sería 75 kg.

La hipótesis nula y alternativa de este ejemplo.

• H0: µ = 75 el peso medio de nuestros alumnos es igual al valor de referencia dado (no hay diferencia) –> valor-p α>0,05.

• H1: µ ≠ 75 el peso medio no es igual al valor de referencia dado (hay una diferencia) –> valor-p α<0,05.

Primero cargaríamos el dataframe alumnos_UC.RData.

load("D:/G2040/TEMA_3_ESDA_Graficos/datos/alumnos_uc.RData")
str(alumnos_uc)

## 'data.frame':    100 obs. of  9 variables:
##  $ peso                 : num  88.1 81.2 89.1 64.3 87.3 ...
##  $ edad                 : int  18 22 19 19 23 20 20 25 25 21 ...
##  $ procedencia          : Factor w/ 6 levels "Bahía","Besaya",..: 3 2 1 2 5 3 5 5 1 1 ...
##  $ genero               : chr  "Femenino" "Masculino" "Masculino" "Femenino" ...
##  $ rama_conocimiento    : int  4 1 3 3 1 5 4 3 1 1 ...
##  $ rama_conocimiento_let: Factor w/ 5 levels "Arte y Humanidades",..: 4 1 3 3 1 5 4 3 1 1 ...
##  $ uso_movil            : int  1 3 1 3 3 3 2 2 2 3 ...
##  $ uso_movil_let        : Ord.factor w/ 3 levels "Escaso"<"Moderado"<..: 1 3 1 3 3 3 2 2 2 3 ...
##  $ grupo_sanguineo      : Factor w/ 8 levels "A-","A+","AB-",..: 2 8 8 1 2 1 2 8 8 8 ...

A continuación aplicamos la prueba t:

library(stats)
t.test(alumnos_uc$peso, 
       mu = 75)                                               # Variable de referencia

## 
##  One Sample t-test
## 
## data:  alumnos_uc$peso
## t = 4.9377, df = 99, p-value = 3.208e-06
## alternative hypothesis: true mean is not equal to 75
## 95 percent confidence interval:
##  78.16056 82.40721
## sample estimates:
## mean of x 
##  80.28389

El significado de la salida es el siguiente:

• t: el estadístico t, calculado como (x – μ) / (s√n) = 4.9377.

• df: los grados de libertad, calculados como n-1 = 99.

• p-value: El valor-p- que corresponde a un estadístico t de 4.9377 con 11 grados de libertad (en este caso, p-value = 3.208e-06).

• 95 percent confidence interval: El intervalo de confianza con una probabilidad del 95%.

Dado que el valor-p obtenido es menor que 0.05 (α =.05), tenemos que rechazar la hipótesis nula, es decir, el peso medio de nuestros estudiantes es diferente del peso medio recomendado por la OMN.

📝 Actividad de evaluación continua:

• Importa el dataframe zonas_verdes y comprueba si la superficie de los parques es significativamente equivalente a 7.5 \(km^2\).

• Explica el significado de esos valores

Comparación de la media de dos muestras independientes

La prueba t para muestras independientes se utiliza para comparar las medias de dos grupos o muestras independientes y determinar si hay una diferencia significativa entre esas medias. Es una de las pruebas más populares y utilizadas.

Fuente: DataLab (https://datatab.es/tutorial/t-test)

Supongamos que queremos saber si el género de los alumnos influye en el peso, es decir, si hombres y mujeres tienen el mismo peso o no. Antes de proceder a la prueba, es conveniente realizar algunos análisis previos para verificar ciertos supuestos:

• Conviene realizar una visualización inicial de los datos.

• Suposición 1: ¿Son las dos muestras independientes?

• Suposición 2: ¿Los datos de cada uno de los 2 grupos siguen una distribución normal?

• Suposición 3: tienen las dos muestras la misma varianza?

• Elección de pruebas unilaterales o bilaterales.

Visualización de los datos

Para detectar diferencias, se puede utilizar un gráfico de caja y bigotes incorporando una variable cuantitativa que aparecerá en el eje horizontal. Este gráfico requiere dos argumentos:

boxplot(formula, data)

• formula: utiliza el formato y ~ x, donde y es el nombre de la variable continua y x la variable cualitativa usada como agrupación.

• data: de dónde se extraen las variables.

boxplot(peso ~ genero,
        alumnos_uc,
        notch = TRUE)
abline(h = mean(alumnos_uc$peso), col = 2, lwd = 2)                               # Añade el valor medio

Con el argumento notch = TRUE se representan los intervalos de confianza al 95% para la mediana (hendiduras en torno a la mediana); si no se superponen las medianas son diferentes estadísticamente y podemos continuar con la prueba. Por defecto, los diagramas de caja se dibujan con el orden con el que aparecen los factores en los datos. Los datos se puede reordenar aplicando la función reorder() a cualquier medida, como la mediana o la media de los datos.

# De más bajo a más alto
mediana <- with(alumnos_uc, reorder(genero, peso, median)) 
boxplot(alumnos_uc$peso ~ mediana, las = 1)

# De más alto a más bajo
mediana <- with(alumnos_uc, reorder(genero, -peso, median)) # Equivalente
boxplot(alumnos_uc$peso ~ mediana, las = 1)

Si se asigna el gráfico a un objeto, se puede obtener una lista con diferentes componentes.

res <- boxplot(peso ~ genero,
               alumnos_uc)

Estos componentes son:

res

• stats: cada columna representa el bigote inferior, el primer cuartil, la mediana, el tercer cuartil y el bigote superior de cada grupo.

• n: número de observaciones de cada grupo.

• conf: cada columna representa los extremos inferior y superior del intervalo de confianza de la mediana.

• out: número total de valores atípicos.

• group: número total de grupos.

• names: nombres de cada grupo.

Finalmente, se puede recrear el mismo boxplot (res) con la función bxp().

library(graphics)
bxp(res)

Podemos agregar puntos representando a los datos

boxplot(peso ~ genero,                                                       
        data = alumnos_uc,                                         # Dataframe de origen
        col = "white")                                     # Color del fondo de las cajas

stripchart(peso ~ genero,                                  # Equivalente
           data = alumnos_uc,
           vertical = TRUE,                                # Diagrama vertical
           method = "jitter",                              # Evita dibujar los datos atípicos. Distribuye aleatoriamente los puntos
           pch = 19,                                       # Tipo de símbolo
           add = TRUE,                                     # Superpone el stripchart al diagrama de caja y bigotes
           col = "red")                                    # Color del símbolo

Gráfico de violín con datos de dos grupos en lados diferentes. Primero, extraemos las dos variables:

borrame <- alumnos_uc[  , c(1,4)]

A continuación se dividen los datos en dos grupos

subset_masculino <- borrame[borrame$genero == "Masculino",  ]
subset_femenino <- borrame[borrame$genero == "Femenino",  ]

library("vioplot")

vioplot(subset_masculino,
        plotCentre = "line",                               # Mediana con una línea
        side = "right",                                    # Lado derecho
        col = "#5773CC",
        xlab = "Géneros",  
        ylim = c(60, 120))                                  # Color del lado derecho
vioplot(subset_femenino,
        plotCentre = "line",                               # Mediana con una línea
        side = "left",                                     # Lado izquierdo
        col = "#FFB900",                                   # Color del lado izquierdo
        add = TRUE)                                        # Sobre el gráfico anterior

legend("topleft",
       legend = c("Masculino", "Femenino"),
       fill = c("#5773CC", "#FFB900"))

Suposición 1: ¿Son las dos muestras independientes?

Sí, ya que se supone que hombres y mujeres poseen diferente configuración corporal.

Suposición 2: ¿Los datos de cada uno de los 2 grupos siguen una distribución normal?

Para ello utilizaremos la prueba de normalidad de Shapiro-Wilk. En esta prueba:

• Hipótesis nula H0: los datos se distribuyen normalmente

• Hipótesis alternativa H1: los datos no se distribuyen normalmente

Por lo tanto:

• Si el valor-p es superior a α =.05 (α>0,05), aceptamos la hipótesis nula H0: la variable sigue una distribución normal.

• Si el valor-p es inferior a α =.05 (α<0,05), rechazamos la hipótesis nula H0 y aceptamos la hipótesis alternativa H1: la variable no sigue una distribución normal.

✅ EJEMPLOS:

Ejemplo 1: resultados sobre datos con una distribución normal

set.seed(0)
datos <- rnorm(100)                              # 100 valores aleatorios con una distribución normal
hist(datos, col='steelblue')

shapiro.test(datos)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos
## W = 0.98957, p-value = 0.6303

El `p-value’ es 0.6303. Dado que es superior a .05, se asume que la muestra proviene de una población que tiene una distribución normal.

Ejemplo 2: resultados sobre datos con una distribución no normal

datos <- rpois(n=100, lambda=3)                                 # 100 valores aleatorios con una distribución de Poisson
hist(datos, col='steelblue')

shapiro.test(datos)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos
## W = 0.90768, p-value = 3.308e-06

El `p-value’ es 0.8.123e-07. Dado que es inferior a .05, se asume que la muestra proviene de una población que no tiene una distribución normal.

En el caso de trabajar con las variables de un dataframe, es posible replicar ese procedimiento dibujando sendos histogramas aplicando la función with() y posteriormente la función hist().

par(mfrow = c(1, 2))

with(alumnos_uc, hist(peso[genero == "Masculino"], main= "Masculino", col='steelblue'))
with(alumnos_uc, hist(peso[genero == "Femenino"], main= "Masculino", col='steelblue'))  

par(mfrow = c(1, 1))

La prueba de Shapiro se obtendría de la siguiente manera.

with(alumnos_uc, shapiro.test(peso[genero == "Masculino"]))                     

## 
##  Shapiro-Wilk normality test
## 
## data:  peso[genero == "Masculino"]
## W = 0.93866, p-value = 0.02566

with(alumnos_uc, shapiro.test(peso[genero == "Femenino"]))               

## 
##  Shapiro-Wilk normality test
## 
## data:  peso[genero == "Femenino"]
## W = 0.93184, p-value = 0.002899

Dado que en ambos casos el p-value es inferior a 0.05 (α<0,05), rechazamos la hipótesis nula H0 y aceptamos la hipótesis alternativa H1: la distribución de ambas variables no es normal. Si los datos no se distribuyen normalmente, es posible optar por

• La prueba de rango Wilcoxon (no paramétrica).

• Proceder a una de las siguientes transformaciones para que la variable de respuesta se acerque a una distribución normal.

  • Transformación logarítmica de la variable respuesta mediante la función log(y).

  • Transformación calculando la raíz cuadrada de la variable respuesta √y.

  • Transformación en la raíz cúbica.

Suposición 3: tienen las dos muestras la misma varianza?

Para ello se usa la prueba F mediante la función var.test(). Las hipótesis nula y alternativa son las siguientes:

• Hipótesis nula H0: las varianzas de ambas variables son iguales.

• Hipótesis alternativa H1: las varianzas de ambas variables son diferentes.

Por lo tanto:

• Si el valor de p es superior a α =.05 (α >.05), aceptamos la hipótesis nula: las varianzas de ambas variables son iguales.

• Si el valor de p es inferior a α =.05 (α <.05), se rechaza la hipótesis nula/aceptamos la alternativa: las varianzas de ambas variables no son iguales (son diferentes).

res.ftest <- var.test(peso ~ genero, 
                      data = alumnos_uc)
res.ftest

## 
##  F test to compare two variances
## 
## data:  peso by genero
## F = 1.0953, num df = 57, denom df = 41, p-value = 0.7672
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6082483 1.9193511
## sample estimates:
## ratio of variances 
##           1.095283

La interpretación de esos resultados es la siguiente

• F: el valor del estadístico F, en este caso 1,0953.

• num df, denom df: grados de libertad del numerador (57) y del denominador (41) para el estadístico F, calculados como \(n_1-1\) \(n_2-1\), respectivamente.

• p-value: valor-p que corresponde a un valor de F de 1,0953 con 57 (numerador) y 41 (denominador) grados de libertad, que es 0,7672.

• 95 percent confidence interval: el intervalo de confianza al 95% correspondiente a la ratio entre las varianzas de ambos grupos[.6082483, 1.9193511]. Dado que la ratio (1.095283) se encuetnra dentro de este intervalo, es probable que el ratio sea verdadero, o sea, que las varianzas sean iguales.

• sample estimates: representa la relación de varianzas entre cada grupo. Si calculáramos las varianzas de ambos grupos y las dividiéramos entre sí, la relación de varianzas debería ser 1.095283.

El valor p de la prueba F es 1,0953, con un p-value de 0.7672; este valor superior al valor crítico correspondiente a un nivel de significación α =.05. En consecuencia, aceptamos la hipótesis nula: las varianzas de ambos conjuntos de datos son iguales.

Elección de pruebas unilaterales o bilaterales.

Si la Hipótesis nula H0 supone que las dos variables bajo análisis tienen igual media, la hipótesis alternativa puede tener 3 posibilidades:

• Hipótesis alternativa H1 (dos colas): las medias de ambas variables no son iguales.

• Hipótesis alternativa H1 (una cola-izquierda): la media de la variable 1 es menor que la media de la variable 2.

• Hipótesis alternativa H1 (una cola-derecha): la media de la variable 1 es mayor que la media de la variable 2.

¿Cuándo elegir “two.sided” (default), “greater” or “less”?

• Pruebas bilaterales o de dos colas. Una prueba de dos colas se asocia a una hipótesis alternativa para la cual se desconoce el signo de la potencial diferencia. Por ejemplo, deseamos comparar las medias de dos muestras A y B, pero no sabemos si A debería ser superior a B o a la inversa. Esto nos lleva a elegir una prueba de dos colas, asociada a la siguiente hipótesis alternativa: H1: media(A) ≠ media(B). Las pruebas de dos colas son con diferencia las más utilizadas.

• Pruebas unilaterales o de una cola. Una prueba de una cola normalmente está asociada a una hipótesis alternativa para la cual se conoce el signo de la potencial diferencia antes de ejecutar el experimento y la prueba. En el ejemplo descrito más arriba, la hipótesis alternativa referida a una prueba de una cola podría redactarse así: media(A) < media(B) o media(A) > media(B).

Cálculo de la prueba t para dos muestras independientes.

Para realizar la prueba t comparando las medias de dos muestras independientes, la función R t.test() se usa de la siguiente manera:

t.test(x, y, alternative = “two.sided”, mu=0, paired = FALSE, var.equal = FALSE, conf.level=0.95)

En este caso

• x,y: vectores numéricos

• alternative: La hipótesis alternativa. El valor permitido es uno de “two.sided” (default), “greater” or “less”.

• var.equal: variable lógica que indica si se deben tratar las dos varianzas como iguales. Si es VERDADERO, entonces la varianza agrupada se utiliza para estimar la varianza, de lo contrario se utiliza la prueba de Welch.

• mu: valor que se asume es la verdadera diferencia en la media.

• paired: si se uso o no la prueba t-test por pares.

• conf.level: el nivel de confianza para usar la prueba.

La aplicación de la prueba responderá a la cuestión ¿hay diferencias significativas en el peso de hombres y mujeres?

prueba_t <- t.test(peso ~ genero, 
              data = alumnos_uc, 
              var.equal = TRUE)
prueba_t

## 
##  Two Sample t-test
## 
## data:  peso by genero
## t = 0.66771, df = 98, p-value = 0.5059
## alternative hypothesis: true difference in means between group Femenino and group Masculino is not equal to 0
## 95 percent confidence interval:
##  -2.862956  5.766466
## sample estimates:
##  mean in group Femenino mean in group Masculino 
##                80.89362                79.44187

La salida consta de los siguientes elementos:

• t: valor del estadístico t (t = 0.66771).

• df: grados de libertad (df = 98).

• p-value: nivel de significación de la prueba t (p-value = 0.5059).

• conf.int: intervalo de confianza de la media con una probabilidad del 95% (conf.int = [-2.862956 5.766466]).

• sample estimate: promedios de los dos grupos sometidos a comparación (caso de un test independiente) o la diferencia en las medias (caso de un test pareado).

En consecuencia, si las hipótesis nula y alternativa son:

• H0: µ1 = µ2 (las medias son iguales).

• HA: µ1 ≠µ2 (las medias no son iguales).

Entonces

• Si el valor de p es inferior a 0.05 (α <.05)se rechaza la hipótesis nula/acepta la hipótesis alternativa: las medias de ambas variables no son iguales.

• Si el valor de p es mayor que 0.05 (α >.05)se acepta la hipótesis nula: las medias de ambas variables son iguales.

Al crear un objeto denominado res podemos acceder a cada elemento de la salida por separado:

prueba_t$p.value                                    # valor p

## [1] 0.5058904

prueba_t$estimate                                   # promedios

##  mean in group Femenino mean in group Masculino 
##                80.89362                79.44187

prueba_t$conf.int                                   # intervalo de confianza

## [1] -2.862956  5.766466
## attr(,"conf.level")
## [1] 0.95

En caso de que los resultados ofrecieran una diferencia significativa, también se podría analizar si el peso de los hombres es superior al de las mujeres o viceversa.

• El peso de los hombres es inferior al de las mujeres:

t.test(peso ~ genero, data = alumnos_uc,
        var.equal = TRUE, alternative = "less")

## 
##  Two Sample t-test
## 
## data:  peso by genero
## t = 0.66771, df = 98, p-value = 0.7471
## alternative hypothesis: true difference in means between group Femenino and group Masculino is less than 0
## 95 percent confidence interval:
##      -Inf 5.062194
## sample estimates:
##  mean in group Femenino mean in group Masculino 
##                80.89362                79.44187

• El peso de los hombres es superior al de las mujeres:

t.test(peso ~ genero, data = alumnos_uc,
        var.equal = TRUE, alternative = "greater")

## 
##  Two Sample t-test
## 
## data:  peso by genero
## t = 0.66771, df = 98, p-value = 0.2529
## alternative hypothesis: true difference in means between group Femenino and group Masculino is greater than 0
## 95 percent confidence interval:
##  -2.158684       Inf
## sample estimates:
##  mean in group Femenino mean in group Masculino 
##                80.89362                79.44187

📝 Actividad de evaluación continua:

ACTIVIDAD 1. Aplica la prueba t al análisis de la relación entre la superficie según el parque esté o no en obras. Para ello:

• Comprueba si la variable obras es un factor o no. Si no lo es, conviértela en factor.

zonas_verdes$obras <- factor(zonas_verdes$obras, 
                             levels = c(0,1),
                             labels = c("Sin obras", "Con obras"))

• Elabora un gráfico de caja y bigotes

boxplot(superficie ~ obras,
        zonas_verdes,
        notch = TRUE)
abline(h = mean(zonas_verdes$superficie), col = 2, lwd = 2)                               # Añade el valor medio

• ¿Siguen los dos grupos (con y sin obras) una distribución normal?.

with(zonas_verdes, shapiro.test(superficie[obras == "Sin obras"]))                     
with(zonas_verdes, shapiro.test(superficie[obras == "Con obras"]))               

• Comprueba si tienen las mismas varianzas.

var.test(superficie ~ obras, 
        data = zonas_verdes)

• ¿Hay diferencias significativas en el promedio de la superficie de los jardines con y sin obras?

t.test(superficie ~ obras, 
      data = zonas_verdes, 
      var.equal = TRUE)

• ¿Es la superficie de los jardines con obras mayor que la de los jardines sin obras?

t.test(superficie ~ obras, 
      data = zonas_verdes,
      var.equal = TRUE, 
      alternative = "greater")

ACTIVIDAD 2. Aunque en España no son habituales los conches de transmisión automática, se considera que son más cómodos y tienen una menor propensión a sufrir daños mecánicos. En este apartado analizaremos además si tienen ventaja sobre los coches de transmisión manual en los apartados consumo (mpg) y aceleración (qsec)

• En primer lugar, carga el fichero mtcars y crea con él un objeto denominado coches.

data(mtcars)
coches <- mtcars

• A continuación, transforma en factor la variables am.

coches$am <- factor(coches$am, 
                    levels = c(0,1),
                    labels = c("automático","manual"))

• Elabora un gráfico de caja y bigotes que relacione las variable amy mpg. Incluye una línea que corresponda al valor mediano de mpg

boxplot(mpg ~ am,
        coches)
abline(h = median(coches$mpg), col = 2, lwd = 2)                               

• Construye un gráfico de violín que represente el consumo en función del tipo de transmisión. Para ello, deberás crear dos subconjuntos en función del tipo de transmisión. Posteriormente, elabora el gráfico de violín. Incluye una leyenda a tu gusto.

coches.automáticos <- coches[coches$am == "automático",  ]
coches.manuales <- coches[coches$am == "manual",  ]

Gráfico de violín.

• Aplica la prueba de la t de Student para determinar si hay diferencias significativas en el consumo de los vehículos en función del tipo de transmisión. ¿Es estadísticamente significativo el consumo medio en función del tipo de transmisión? En consecuencia, ¿qué vehículo recomendarías?

t.test(mpg ~ am, 
       data = coches, 
       var.equal = TRUE)

• Repite la aplicación de la prueba t de Student para verificar si el tipo de transmisión (variable am) influye también en la aceleración del vehículo (variable qsec). A partir de los resultados de la prueba, ¿influye o no?

t.test(qsec ~ am, 
       data = coches, 
       var.equal = TRUE)

• Comprueba que la variable mpg sigue una distribución normal aplicando la prueba de Shapiro

shapiro.test(coches$mpg) 

Prueba de Wilcoxon de dos muestras independientes

Es conocida también como prueba de la suma de rangos de Wilcoxon o prueba de Mann-Whitney. Es la alternativa no paramétrica a la prueba t de dos muestras no pareada, y se usa para comparar dos grupos independientes de datos cuando éstos no se distribuyen normalmente. La idea que subyace es que si las dos muestras comparadas proceden de la misma población, al juntar todas las observaciones y ordenarlas de menor a mayor, cabría esperar que las observaciones de una y otra muestra estuviesen intercaladas aleatoriamente. Por lo contrario, si una de las muestras pertenece a una población con valores mayores o menores que la otra población, al ordenar las observaciones, estas tenderán a agruparse de modo que las de una muestra queden por encima de las de la otra.

wilcox.test(x, y, alternative = “two.sided”)

donde:

• x,y: vectores numéricos.

• alternative: hipótesis alternativa, siendo los valores permitidos “two.sided” (default), “greater” or “less”.

Las hipótesis nula y alternativa se enuncian así:

• Hipótesis nula H0: los rangos de ambas poblaciones son iguales.

• Hipótesis alternativa Ha: los rangos de ambas poblaciones no son iguales.

Por lo tanto

• Si el valor de p es inferior a α =.05, se rechaza la hipótesis nula/acepta la hipótesis alternativa: las medias de ambas variables no son iguales.

• Si el valor de p es mayor a α =.05, se acepta la hipótesis nula: las medias de ambas variables son iguales.

Para realizarlo

prueba_wilcoxon <- wilcox.test(peso ~ genero, 
                   data = alumnos_uc)
prueba_wilcoxon

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  peso by genero
## W = 1312, p-value = 0.5138
## alternative hypothesis: true location shift is not equal to 0

Si queremos imprimir sólo el valor p:

prueba_wilcoxon$p.value

## [1] 0.5137662

El valor p de la prueba es 0.5138, mayor que el nivel de significación alfa = 0,05. Se acepta la hipótesis nula: el peso medio de los hombres no es significativamente diferente del peso de las mujeres (el tamaño medio no es igual).

También es posible replicar si el peso medio de los hombres es menor que el de las mujeres:

wilcox.test(peso ~ genero, data = alumnos_uc, 
        alternative = "less")

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  peso by genero
## W = 1312, p-value = 0.7454
## alternative hypothesis: true location shift is less than 0

O la opción contraria:

wilcox.test(peso ~ genero, data = alumnos_uc,
        alternative = "greater")

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  peso by genero
## W = 1312, p-value = 0.2569
## alternative hypothesis: true location shift is greater than 0

📝 Actividad de evaluación continua:

ACTIVIDAD 3

• Repite, ahora aplicando el test de Wilcoxon, el análisis de la relación entre la superficie según el parque esté o no en obras.

wilcox.test(superficie ~ obras, 
            data = zonas_verdes)

• Además, para estar seguro de los resultados, comprueba que existe una relación significativa entre las variables mpg y am.

wilcox.test(mpg ~ am, 
            data = coches)

Comparación de las medias de dos grupos por parejas (prueba t para muestras relacionadas)

La prueba t para muestras relacionadas, también conocida como de muestras dependientes, se utiliza para comparar las medias de dos grupos dependientes.

• En una muestra pareada (muestra dependiente) los valores medidos están disponibles por pares. Los pares se crean, por ejemplo, mediante mediciones repetidas sobre los mismos sujetos.

• Las muestras independientes (muestra no pareada) son el resultado de personas y mediciones independientes entre sí.

Fuente: DataLab (https://datatab.es/tutorial/t-test)

Como ejemplo de un conjunto de datos que se puede someter a un análisis por grupos pareados, supongamos que en una localidad se instala una nueva industria potencialmente contaminante. Para determinar si su actividad afecta a la calidad del agua, se midió la concentración de un determinado contaminante en el agua durante las 10 semanas anteriores al comienzo de la actividad y durante las 10 semanas posteriores. Los datos son los siguientes:

antes <- c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)
despues <- c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)

A partir de ellos se crea un nuevo dataframe, que comprende una variable cualitativa (grupo) y una variable cuantitativa (concentración).

datos_calidad_agua <- data.frame(
  grupo = rep(c("antes", "despues"), each = 10),
                concentracion = c(antes,  despues))

rm(antes, despues)

La función de R para llevar a cabo una prueba t comparando las medias de dos parejas de grupos es la siguiente:

t.test(x, y, paired = TRUE, alternative = “two.sided”)

en la que

• x,y: son los dos grupos.

• paired: una valor lógico que indica (TRUE) que se desea una prueba pareada.

• alternative: la hipótesis alternativa con tres posibilidades: “two.sided” (por defecto), “greater” or “less”.

Con estos datos podemos calcular unos estadísticos básicos:

promedio <- aggregate(concentracion ~ grupo, 
                      data = datos_calidad_agua, 
                      FUN = mean)
colnames(promedio) <- c("Grupo", "Media")
promedio

##     Grupo  Media
## 1   antes 199.16
## 2 despues 393.65

Al igual que ocurre en el caso de las muestras independientes, es necesario comprobar las mismas asumpciones, pero en este caso serán omitidas. Para realizar la prueba t para muestras pareadas:

prueba_muestras_pareadas <- t.test(concentracion ~ grupo, 
                                   data = datos_calidad_agua, 
                                   paired = TRUE)
prueba_muestras_pareadas

## 
##  Paired t-test
## 
## data:  concentracion by grupo
## t = -20.883, df = 9, p-value = 6.2e-09
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -215.5581 -173.4219
## sample estimates:
## mean difference 
##         -194.49

Los componentes del cuadro superior son idénticos al obtenido al analizar muestras independientes:

• t: el valor del estadístico t (t = 20.88).

• df: los grados de libertad (df= 9).

• p-value: el nivel de significación de la prueba t obtenido (p-value = 6.210^{-9}).

• conf.int: el intervalo de confianza de la media de las diferencias con un nivel de significación del 95% (173.42, 215.56).

• sample estimates: la media de las diferencias entre las parejas (mean = 194.49).

El valor de p de la prueba \(6.2*10^{-9}\) es menor que el nivel de significación (alpha = 0.05). Por lo tanto, podemos rechazar la hipótesis nula y concluir que el valor medio de la concentración de ese contaminante después de la instalación de la empresa fue diferente de su valor medio antes de su instalación.

También es posible acceder a los elementos de esta lista por separado

prueba_muestras_pareadas$p.value                         

## [1] 6.200298e-09

prueba_muestras_pareadas$estimate

## mean difference 
##         -194.49

prueba_muestras_pareadas$conf.int

## [1] -215.5581 -173.4219
## attr(,"conf.level")
## [1] 0.95

Si quisiéramos averiguar si la concentración antes es inferior a la concentración después utilizaríamos la siguiente sintaxis:

t.test(concentracion ~ grupo, 
       data = datos_calidad_agua, 
       paired = TRUE,
       alternative = "less")

## 
##  Paired t-test
## 
## data:  concentracion by grupo
## t = -20.883, df = 9, p-value = 3.1e-09
## alternative hypothesis: true mean difference is less than 0
## 95 percent confidence interval:
##       -Inf -177.4177
## sample estimates:
## mean difference 
##         -194.49

Alternativamente, si lo que se quiere es si fue mayor antes que después:

t.test(concentracion ~ grupo, 
       data = datos_calidad_agua, 
       paired = TRUE,
       alternative = "greater")

## 
##  Paired t-test
## 
## data:  concentracion by grupo
## t = -20.883, df = 9, p-value = 1
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
##  -211.5623       Inf
## sample estimates:
## mean difference 
##         -194.49

✅ EJEMPLO DE ANÁLISIS DE LA CALIDAD DEL AIRE:

El fichero calidad_aire_covid.csv contiene los valores medios mensuales de varios contaminantes en varias localidades de Cantabria durante los años 2019 y 2020. Dado que en este último año, a causa de la epidemia de COVID, se produjeron periodos de confinamiento, es posible que la ausencia de actividad haya afectado a los niveles de los contaminantes. A título de ejemplo, se analizará una estación concreta y un contaminante.

Primero, los alumnos deberán cargar el fichero *.csv. Se utiliza la función read.csv2() dado que el fichero original está escrito en un formato “europeo”

calidad_aire_covid <- read.csv2("https://personales.unican.es/rasillad/docencia/G2040/TEMA_6/calidad_aire_covid.csv")

A continuación, como ejemplo, se selecciona la estación de contaminación denominada Tetuan (situada en Santander) y el contaminante no2 (dióxido de nitrógeno).

tetuan_no2 <- subset(calidad_aire_covid,
                site == "TETUAN" & contaminante == "no2")

Sobre ese subconjunto de datos, se calcula el valor medio según la variable grupo (en realidad, separa los periodos anterior y posterior al confinamiento por COVID-19) y se representa mediante un gráfico de caja y bigotes.

media <- aggregate(valor ~ grupo, 
                   data = tetuan_no2, 
                   FUN = mean)
colnames(media) <- c("Grupo", "Media")
media

##   Grupo    Media
## 1  2019 12.97732
## 2  2020 11.74122

boxplot(valor ~ grupo,
        tetuan_no2)
abline(h = mean(tetuan_no2$valor), col = 2, lwd = 2) 

Por último, se someten esos datos a una prueba t.

t.test(valor ~ grupo, 
      data = tetuan_no2, 
      paired = TRUE)

## 
##  Paired t-test
## 
## data:  valor by grupo
## t = 1.4198, df = 11, p-value = 0.1834
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.6801795  3.1523847
## sample estimates:
## mean difference 
##        1.236103

Test de Wilcoxon

Respecto a la aplicación del test de Wilcoxon, en este caso debe incluirse el argumento paired = TRUE en la sintaxis ya conocida de esta prueba:

wilcox.test(x, y, paired = TRUE, alternative = “two.sided”)

• x, y: los valores a comparar.

• paired: valor lógico especificando que se quiere realizar un test pareado.

• alternative: la hipótesis alternativa con tres posibles valores: “two.sided” (default), “greater” or “less”.

Como ejemplo

prueba_wilcoxon_pares <- wilcox.test(valor ~ grupo, 
                                     data = tetuan_no2, 
                                     paired = TRUE)
prueba_wilcoxon_pares

## 
##  Wilcoxon signed rank exact test
## 
## data:  valor by grupo
## V = 61, p-value = 0.09229
## alternative hypothesis: true location shift is not equal to 0

El valor p de la prueba es 0.09229, mayor que el nivel de significación (alpha = 0.05). En consecuencia, la mediana de la concentración antes de la instalación de la empresa no fue significativamente diferente de la concentración después del confinamiento.

📝 ACTIVIDAD DE EVALUACIÓN CONTINUA:

ACTIVIDAD 4

Utilizando los datos del fichero [calidad_aire_covid.csv], cada alumno deberá elegirá una localidad de Cantabria, comparando los niveles medios de cada contaminante en 2019 y 2020, para determinar si el confinamiento tuvo algún efecto sobre la calidad del aire, aplicando las pruebas correspondientes. La actividad deberá ser enviada al profesor en forma de fichero word.

COMPARACIÓN DE LOS PROMEDIOS DE MÁS DE DOS GRUPOS: Análisis de varianza (ANOVA):

Para comparar las medias de más de dos grupos se utiliza el Análisis de varianza, conocido como ANOVA, que tiene dos variantes:

Prueba análisis de varianza unidireccional (one-way ANOVA)

El análisis de varianza unidireccional (one-way ANOVA), también conocido como ANOVA de un factor, es una extensión de la prueba t de dos muestras independientes. Se utiliza para determinar si existe o no una diferencia estadísticamente significativa entre las medias de tres o más grupos independientes. Se llama ANOVA unidireccional porque estamos analizando cómo una variable predictora impacta una variable de respuesta.

###¿Cuáles son las hipótesis de prueba de ANOVA?

• Hipótesis nula H0: los promedios de todas las variables es igual.

• Hipótesis alternativa H1: al menos el promedio de una variable es diferentes al de las restantes

Por lo tanto

• Si el valor de p es inferior a 0.05 se rechaza la hipótesis nula/acepta la hipótesis alternativa: la media de una variable es diferente de las restantes (no es igual).

• Si el valor de p es mayor a 0.05, se acepta la hipótesis nula: las medias de todas las variables son iguales.

¿Cómo funciona la prueba ANOVA unidireccional?

• Calcular la varianza común, que se denomina varianza (\(S2_{dentro}\)) de las muestras o varianza residual.

• Calcular la varianza entre las medias muestrales de la siguiente manera:

  • Calcular la media de cada grupo

  • Calcular la varianza entre medias muestrales (\(S2_{entre}\)).

• Calcular el estadístico F como la relación de \(S2{entre}\)/\(S2{dentro}\).

Pasos iniciales

En el ejemplo que trabajaremos a continuación, se analizará si el nivel de uso de móvil tiene relación con el peso de los alumnos.

Es conveniente iniciar el análisis calculando estadísticos de resumen por grupos: número de casos, media, desviación típica etc…

aggregate(peso ~ uso_movil_let, data = alumnos_uc, length)

##   uso_movil_let peso
## 1        Escaso   37
## 2      Moderado   30
## 3       Intenso   33

aggregate(peso ~ uso_movil_let, data = alumnos_uc, mean)

##   uso_movil_let     peso
## 1        Escaso 83.35045
## 2      Moderado 79.43244
## 3       Intenso 77.61965

Igualmente, es conveniente visualizar los datos.

boxplot(peso ~ uso_movil_let, 
        data = alumnos_uc) 

Implementación de la prueba ANOVA unidireccional

Para comprobar si existe alguna diferencia significativa entre el valor promedio en cada categoría se usa la función aov(). El resultado no es la tabla del ANOVA, ya que para obtenerla hay que aplicar la función summary al resultado de aov.

prueba_aov <- aov(peso ~ uso_movil_let, 
               data = alumnos_uc)                 

summary(prueba_aov)                                               

##               Df Sum Sq Mean Sq F value Pr(>F)  
## uso_movil_let  2    604   302.0   2.729 0.0703 .
## Residuals     97  10733   110.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En una tabla ANOVA:

• En la primera columna, dos etiquetas: el nombre del factor, en este caso uso_movil_let, y Residuals, que representa los errores o residuos del ANOVA.

• La segunda columna, etiquetada Df, proporciona los grados de libertad correspondientes al factor (su número de niveles menos 1) y a los residuos (el número de individuos en la tabla, menos el número de niveles del factor).

• La tercera columna, Sum Sq, muestra las sumas de los cuadrados del factor, \(SS_{Tr}\), y de los residuos, \(SS_E\).

• La cuarta columna, Mean Sq, contiene las medias de los cuadrados del factor, \(MS_{Tr}\), y de los residuos, \(MS_E\).

• La quinta columna, F value, nos da el valor del estadístico de contraste.

• En la sexta columna, **Pr(>F)*, aparece el p-valor del contraste.

• La séptima columna, sin etiqueta, indica el nivel de significación del p-valor según el código usual, explicado en la última línea del resultado.

A mayor número de asteriscos, más significativo es el p-valor y por lo tanto es más fuerte la evidencia de que las medias comparadas no son todas iguales.

¿Cómo interpretamos el resultado de la prueba one-way ANOVA?

Como el valor- p es mayor que el nivel de significación 0,05, se acepta la hipótesis nula: los promedios de la variable peso son iguales para cada categoría de la variable uso de móvil. En consecuencia, no existen diferencias significativas entre los grupos resaltados con “*” en el resumen del modelo.

Pero debe tenerse en cuenta que este resultado no significa que se haya obtenido evidencia de que cada categoría es significativamente diferente de las restantes. Para ello es necesario realizar comparaciones por pares, determinando si la diferencia media entre pares específicos de grupos es estadísticamente significativa.

Una posibilidad es la prueba HSD de Tukey. La función TukeyHD() toma como argumento el ANOVA ajustado.

TukeyHSD(prueba_aov)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = peso ~ uso_movil_let, data = alumnos_uc)
## 
## $uso_movil_let
##                       diff        lwr       upr     p adj
## Moderado-Escaso  -3.918010 -10.069284 2.2332644 0.2879009
## Intenso-Escaso   -5.730798 -11.725676 0.2640806 0.0641228
## Intenso-Moderado -1.812788  -8.128784 4.5032083 0.7738521

Este procedimiento proporciona una tabla con las siguientes columnas:

• diff: diferencias entre las medias de cada pareja de grupos.

• lwr, upr: valor más bajo/más alto del intervalo de confianza al 95 % (por defecto)

• p adj: valor-p después del ajuste para comparaciones múltiples.

De acuerdo con esta tabla, las diferencias entre las medias de los pares de categorías no son estadísticamente significativas, ya que el valor-p ajustado es superior a 0.05.

Otra posibilidad es usar la función glht(), incluida en el paquete multcomp. La sintaxis es la siguiente:

glht(model, lincft)

• model: modelo ANOVA, por ejemplo un objeto obtenido tras aplicar la función aov().

• lincft: especificación de las hipótesis lineales a ser probadas. Las comparaciones múltiples en los modelos ANOVA se especifican mediante objetos devueltos por la función mcp().

library(multcomp)

summary(glht(prueba_aov, 
             linfct = mcp(uso_movil_let = "Tukey")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = peso ~ uso_movil_let, data = alumnos_uc)
## 
## Linear Hypotheses:
##                         Estimate Std. Error t value Pr(>|t|)  
## Moderado - Escaso == 0    -3.918      2.584  -1.516   0.2878  
## Intenso - Escaso == 0     -5.731      2.519  -2.275   0.0641 .
## Intenso - Moderado == 0   -1.813      2.654  -0.683   0.7738  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Obsérvese que los resultados son iguales, aunque tabla es algo diferente, ya que, además de la probabilidad (Pr(>|t|)), incluye el valor de la prueba t en cada pareja, el error típico etc…

Por último, la función pairewise.t.test() también puede usarse para calcular comparaciones pareadas entre los niveles de los grupos con diferentes correcciones.

pairwise.t.test(alumnos_uc$edad, 
                alumnos_uc$uso_movil_let,
                 p.adjust.method = "BH")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  alumnos_uc$edad and alumnos_uc$uso_movil_let 
## 
##          Escaso Moderado
## Moderado 0.94   -       
## Intenso  0.94   0.94    
## 
## P value adjustment method: BH

El resultado es una tabla de valores p ajustados por el método “Benjamini-Hochberg”.

Al igual que ocurría con la prueba t, la prueba ANOVA se puede aplicar solo cuando:

• Las observaciones son obtenidas de manera independiente y aleatoriamente de una población definida por una serie de grupos.

• Los datos de cada categoría o grupo poseen una distribución normal.

• Estas poblaciones tienen una varianza similar (para comprobarlo se utiliza el test de Levene).

Los datos de cada nivel de factor se distribuyen normalmente.

Para verificar si los datos se distribuyen normalmente, se puede realizar un gráfico de normalidad de los residuos, en el que se representan los cuantiles de los residuos frente a los cuantiles de la distribución normal, acompañada de una línea de referencia. Los residuos deben seguir aproximadamente una línea recta, superpuesta a la línea de referencia.

plot(prueba_aov, 2)

Esa conclusión puede reforzarse aplicando la prueba de Shapiro-Wilk test a esos mismos resíduos:

aov_residuals <- residuals(prueba_aov)

shapiro.test(x = aov_residuals )

## 
##  Shapiro-Wilk normality test
## 
## data:  aov_residuals
## W = 0.9578, p-value = 0.00281

El valor-p es claramente inferior a 0.05, por lo que la variable analizada no tiene una distribución normal (supuesto de normalidad).

Las poblaciones normales tienen una varianza común.

La gráfica de residuos versus valores ajustados se puede usar para verificar la homogeneidad de las varianzas.

plot(prueba_aov, 1)

En ese gráfico no se observan relaciones evidentes entre los residuos y los valores ajustados (la media de cada grupo), por lo que es posible asumir la homogeneidad de las varianzas. No obstante, algunos casos (viviendas) son considerados “outliers”, lo cual podría explicar la ausencia de normalidad y afectar a la homogeneidad de varianza de los datos. En estos casos, es recomendable eliminar esos outliers.

También es posible utilizar la prueba de Bartlett o la prueba de Levene para verificar la homogeneidad de las varianzas. El último es más recomendable, ya que es menos sensible a las desviaciones de la distribución normal.

El test plantea las siguientes hipótesis

• Hipótesis nula (H0): las varianzas entre los grupos son iguales.

• Hipótesis alternativa (HA): la varianzas entre los grupos no son iguales.

Por consiguiente,

• Si el valor-p es inferior a 0.05 se rechaza la hipótesis nula/acepta la hipótesis alternativa: las varianzas de los grupos no son iguales.

• Si el valor-p es mayor a 0.05, se acepta la hipótesis nula: las varianzas de los grupos son iguales.

En R existe la función leveneTest() del paquete car.

library(car)

leveneTest(peso ~ uso_movil_let, 
           data = alumnos_uc)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  2.3513 0.1007
##       97

Del resultado anterior podemos ver que el valor p es mayor que el nivel de significanción de 0.05. Esto significa que no hay evidencia que sugiera que la varianza entre los grupos sea significativamente diferente desde el punto de vista estadístico. Por lo tanto, podemos asumir la homogeneidad de las varianzas en los diferentes grupos o categorías.

¿Podemos realizar la prueba ANOVA cuando se viola el supuesto de homogeneidad de varianza?

Prueba de Welch.

Sí. Una posibilidad es aplicar un procedimiento alternativo (es decir: prueba unidireccional de Welch), que no requiere que se haya implementado esa suposición en la función oneway.test().

oneway.test(peso ~ uso_movil_let, 
            data = alumnos_uc)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  peso and uso_movil_let
## F = 3.0436, num df = 2.000, denom df = 61.542, p-value = 0.0549

De nuevo se confirma la existencia de diferencias signficativas, incluso no asumiendo varianzas iguales.

Otra alternativa es un Pairwise t-tests

pairwise.t.test(alumnos_uc$peso, 
                alumnos_uc$uso_movil_let,
                p.adjust.method = "BH", 
                pool.sd = FALSE)

## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  alumnos_uc$peso and alumnos_uc$uso_movil_let 
## 
##          Escaso Moderado
## Moderado 0.218  -       
## Intenso  0.058  0.524   
## 
## P value adjustment method: BH

Prueba de Kruskall-Wallis

Otra posibilidad es utilizar una alternativa no-paramétrica, la prueba de Kruskal-Wallis (KW rank sum test). Esta prueba usa las siguientes hipótesis nula y alternativa:

• Hipótesis nula (H0): el valor de la mediana es igual entre todos los grupos.

• Hipótesis alternativa (H1): el valor de la mediana no es igual entre todos los grupos.

Por lo tanto

• Si el valor de p es inferior a 0.05 se rechaza la hipótesis nula/acepta la hipótesis alternativa: al menos un grupo tiene una mediana cuyo valor no es igual al de los otros grupos.

• Si el valor de p es mayor que 0.05 se acepta la hipótesis nual: el valor de la mediana es igual entre los grupos.

kruskal.test(peso ~ uso_movil_let, 
             data = alumnos_uc)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  peso by uso_movil_let
## Kruskal-Wallis chi-squared = 5.0048, df = 2, p-value = 0.08189

📝 ACTIVIDAD DE EVALUACIÓN CONTINUA:

ACTIVIDAD 5. El conjunto de datos airquality muestra los valores diarios de ozono y de diferentes variables meteorológicas en Nueva York, entre mayo y septiembre de 1973. El objetivo de esta práctica es comprobar si existen diferencias significativas entre esas variables en función del mes en el que se obtuvieron.

• Para ello, comenzamos cargando el dataset y creando un objeto denominado aq.

• Transforma en factor la variables Month.

aq$Month <- as.factor(aq$Month)

• Dibuja un gráfico de caja y bigotes que muestre los valores mensuales de la variable Ozone y su valor medio durante el periodo de análisis.

boxplot(Ozone ~ Month,
        data = aq)
abline(h = mean(aq$Ozone), col = 2, lwd = 2)                               

• Para comprobar si existe alguna diferencia significativa entre el valor promedio de ozono en cada mes se debe aplicar un análisis ANOVA. Realiza este análisis e interpreta su resultado.

resultados.anova.ozono <- aov(Ozone ~ Month, data = aq)                 
summary(resultados.anova.ozono)                                               

• ¿Qué meses muestran una diferencia significativa con los restantes?

TukeyHSD(resultados.anova.ozono)

• Comprueba si los residuos del model ANOVA se distribuyen según una curva normal.

plot(resultados.anova.ozono, 2)

• Esa conclusión puede reforzarse aplicando la prueba de Shapiro-Wilk test a esos mismos resíduos. ¿Qué criterio sirve para establecer la significación estadística de esta prueba?

residuals.anova.ozono <- residuals(resultados.anova.ozono)
shapiro.test(residuals.anova.ozono)

kruskal.test(Ozone ~ Month,
        data = aq)

ACTIVIDAD 6 El dataset “Ozone”, incluido en el paquete mlbench, consiste en un dataframe con 366 observaciones de 13 variables realizadas en la ciudad de Los Ángeles, en el año 1976, correspondiendo cada observación a un día. Se puede acceder al mismo desde R con los siguientes comandos:

• De ese dataframe extrae un subconjunto con las variables • [4] Valor máximo horario de ozono.

• [5] Altura del geopotencial (en metros) de 500 hPa.

• [6] Velocidad del viento (millas por hora9 en superficie.

• [8] Temperatura del aire en superficie (ºF).

• [10] Altura de la inversión térmica (feet).

• [11] Gradiente de presión en superficie (mm Hg).

• [12] Temperatura en la base de la inversión (ºF).

• [13] Visiblidad (millas).

• Calcula el valor de la mediana de la variable Daily maximum one-hour-average ozone reading y crea una nueva variable en la que el valor 1 corresponda a registros inferiores a la mediana, y el valor 2 consista en los registros iguales o superiores a la mediana.

mediana_V4 <- median(Ozone$V4,
       na.rm = TRUE)

Ozone <- within(Ozone,{
  categorias_ozono <- NA  
  categorias_ozono[V4 < mediana_V4] <- "Menor"         
  categorias_ozono[V4 >= mediana_V4]<- "Mayor o igual"})

• Calcula si existen diferencias significativas, en función de las categorías de esa nueva variable, en los valores medios de las siguientes variables:

• [5] 500 millibar pressure height (m) measured at Vandenberg AFB.

aov_z500 <- aov(V5 ~ categorias_ozono, 
               data = Ozone)
summary(aov_z500) 
TukeyHSD(aov_z500)

• [10] Inversion base height (feet) at LAX.

aov_inversion <- aov(V10 ~ categorias_ozono, 
               data = Ozone)
summary(aov_inversion) 
TukeyHSD(aov_inversion)

• [13] Visibility (miles) measured at LAX.

aov_visib <- aov(V13 ~ categorias_ozono, 
               data = Ozone)
summary(aov_visib) 
TukeyHSD(aov_visib)

• Representa gráficamente esas diferencias mediante gráficos de caja y bigotes.

par(mfrow=c(2, 2))
boxplot(V5 ~ categorias_ozono,
        data = Ozone,
        main = "z500 vs ozono",
        xlab = "Categorías",
        ylab = "hPa")

boxplot(V10 ~ categorias_ozono,
        data = Ozone,
        main = "Inversion vs ozono",
        xlab = "Categorías",
        ylab = "Metros")

boxplot(V13 ~ categorias_ozono,
        data = Ozone,
        main = "Visibilidad vs ozono",
        xlab = "Categorías",
        ylab = "Km")

dev.off()

• Repite el procedimiento anterior con la variable [V5], pero en vez de usar la mediana, utiliza los cuartiles, de tal manera que la nueva variable tendrá 4 niveles (1, 2, 3, 4). Calcula existen diferencias significativas, en función de las categorías de esa nueva variable, en los valores medios y en la varianza de:

cuartiles_V5 <- quantile(Ozone$V5,                    
                      probs = c(0,0.25,0.50,0.75,1),  
                      na.rm=TRUE)                      

Ozone$categorias_z500 <- cut(Ozone$V5,
                             breaks = cuartiles_V5, 
                             include.lowest = TRUE,
                             labels = c("q1","q2","q3","q4"))

• [6] Wind speed (mph) at Los Angeles International Airport (LAX).

• Valores medios

aov_viento <- aov(V6 ~ categorias_z500, 
               data = Ozone)
summary(aov_viento) 
TukeyHSD(aov_viento)

• Varianza (Prueba de Levene)

library(car)
levene_V6 <- leveneTest(V6 ~ categorias_z500, 
                  data = Ozone)

levene_V6

• [8] Temperature (degrees F) measured at Sandburg, CA.

• Valores medios

aov_tem <- aov(V8 ~ categorias_z500, 
               data = Ozone)
summary(aov_tem) 
TukeyHSD(aov_tem)

• Varianza (Prueba de Levene)

levene_V8 <- leveneTest(V8 ~ categorias_z500, 
                  data = Ozone)

levene_V8

• [12] Inversion base temperature (degrees F) at LAX.

• Valores medios

aov_inversion <- aov(V8 ~ categorias_z500, 
               data = Ozone)
summary(aov_inversion) 
TukeyHSD(aov_inversion)

• Varianza (Prueba de Levene)

levene_V8 <- leveneTest(V8 ~ categorias_z500, 
                  data = Ozone)

levene_V8

ACTIVIDAD 7. Utilizando los datos del dataframe zonas_verdes, comprueba si es verdad que existen diferencias sustanciales en cuanto al equipamiento de zonas verdes en los diferentes barrios de la ciudad (superficie vs barrio), e identifica, en su caso, si existen diferencias significativas en algún barrio concreto.

• Carga los datos, si no estuvieran en el Global Environment.

• Calcula el valor medio de la superficie de los parques en cada barrio y representa gráficamente su distribución mediante un gráfico de caja y bigotes.

aggregate(superficie ~ barrio, data = zonas_verdes, mean)

boxplot(superficie ~ barrio, data = zonas_verdes) 

• Aplica la prueba ANOVA e interpreta su resultado.

prueba_aov_barrios <- aov(superficie ~ barrio, data = zonas_verdes)                 

summary(prueba_aov_barrios)     

• Verifica, aplicando la prueba HSD de Tukey, si existen diferencias entre alguno de los barrios con respecto a los restantes.

TukeyHSD(prueba_aov_barrios)

• Verifica si los datos de superficie se distribuyen normalmente creando un gráfico de normalidad de los residuos y aplicando la prueba de Shapiro-Wilk a los residuos.

plot(prueba_aov_barrios, 2)

prueba_aov_barrios_residuals <- residuals(prueba_aov_barrios)
shapiro.test(x = prueba_aov_barrios_residuals )

• Comprueba igualmente el requisito de la homogeneidad de varianzas representando la gráfica de residuos versus valores ajustados y aplicando la prueba de Levene

plot(prueba_aov_barrios, 1)

leveneTest(superficie ~ barrio, 
           data = zonas_verdes)

Finalmente, comprueba, aplicando la prueba de Kruskall-Wallis, que los resultados son similares a ANOVA.

kruskal.test(superficie ~ barrio, 
             data = zonas_verdes)

Two-Way ANOVA Test in R

La prueba ANOVA bidireccional se usa para evaluar simultáneamente el efecto de dos variables de agrupación (dos variables cualitativas A y B) sobre una variable respuesta (variable cuantitativa). Las variables de agrupación también se conocen como factores, mientras que sus diferentes categorías (grupos) se denominan niveles. El número de niveles puede variar entre factores. Las combinaciones de niveles de factores se denominan celdas.

• Cuando los tamaños de muestra dentro de las celdas son iguales, tenemos el llamado diseño “equilibrado”. En este caso, se puede aplicar la prueba ANOVA estándar bidireccional.

• Cuando los tamaños de muestra dentro de cada nivel de las variables independientes no son los mismos (caso de diseños no equilibrados), la prueba ANOVA debe manejarse de manera diferente.

Hipótesis prueba ANOVA bidireccional (dos vías).

Las hipótesis nula son las siguientes:

• No hay diferencia en las medias del factor A

• No hay diferencia en las medias del factor B

• No hay interacción entre los factores A y B

La hipótesis alternativa para los casos 1 y 2 es: las medias no son iguales.

La hipótesis alternativa para el caso 3 es: hay una interacción entre A y B.

Supuestos de la prueba ANOVA bidireccional

La prueba ANOVA bidireccional, como todas las pruebas de ANOVA, asume que las observaciones dentro de cada celda se distribuyen normalmente y tienen varianzas iguales.

Para el siguiente ejercicio seleccionaremos del conjunto de datos original uno nuevo consistente en 3 variables, 2 cualitativas (uso_movil_let y genero) y una cuantitativa (edad).

nuevos_datos <- subset(alumnos_uc, 
                       select=c("uso_movil_let","genero", "edad")) 

El objetivo de este ejercicio consiste en saber si el precio de una vivienda está condicionado por la inclusión en ella de una terraza y un garaje.

Prueba ANOVA bidireccional para diseños equilibrados

Los diseños equilibrados corresponden a la situación en la que tenemos tamaños de muestra iguales dentro de los niveles de nuestros niveles de agrupación independientes.

Análisis inicial de los datos

Una primera posibilidad es la creación de una tabla de frecuencias para conocer la distribución de los datos.

table(nuevos_datos$uso_movil_let, nuevos_datos$edad)

##           
##            18 19 20 21 22 23 24 25
##   Escaso    5  6  5  4  1  6  6  4
##   Moderado  4  3  9  3  3  1  2  5
##   Intenso   4  5  4  4  5  4  3  4

En segundo lugar, se puede representar la dispersión de los datos en un gráfico de caja y bigotes en los que aparezcan los datos agrupados según los niveles de las dos variables cualitativas.

boxplot(edad ~ uso_movil_let * genero, 
        data=nuevos_datos, 
        frame = FALSE, 
        col = c("#00AFBB", "#E7B800"), 
        horizontal = TRUE, 
        las = 2)

En tercer lugar, también se puede aplicar un gráfico denominado de “interacción bidireccional”, que representa la media (u otro estadístico similar) de la variable respuesta para combinaciones bidireccionales de factores, ilustrando las posibles interacciones. En este gráfico se representa la variable cuantitativa según grupos de una de las dos variables cualitativas, al tiempo que se colorea de acuerdo a la segunda variable cualitativa.

interaction.plot(x.factor = nuevos_datos$uso_movil_let,               # Factor representado en el eje x.
                 trace.factor = nuevos_datos$genero,            # Factor representado como líneas.
                 response = nuevos_datos$edad,                # Variable cuantitativa.
                 fun = mean, 
                 type = "b",                                    # p representa sólo puntos, l líneas y b ambos.
                 legend = TRUE, 
                 xlab = "Uso móvil", 
                 ylab="Edad",
                 pch=c(1,19), 
                 col = c("#00AFBB", "#E7B800"))

ANOVA bidireccional

Como ha sido señalado previamente, se quiere saber si la variable cuantitativa depende de los valores que tomen las variables cualitativas. La función aov() sirve para este propósito, mientras que la función summary.aov()resumen el modelo.

prueba_aov2 <- aov(edad ~ uso_movil_let + genero, 
                data = nuevos_datos)
summary(prueba_aov2)

##               Df Sum Sq Mean Sq F value  Pr(>F)   
## uso_movil_let  2    1.4    0.68   0.131 0.87723   
## genero         1   39.4   39.37   7.549 0.00717 **
## Residuals     96  500.7    5.22                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La salida incluye las columnas F value y Pr(>F correspondientes al valor p de la prueba.

De la tabla ANOVA podemos concluir que el uso de móvil no tiene efecto sobre el peso, pero sí el género.

ANOVA bidireccional con efecto de interacción. En el análisis anterior, se considera que el efecto de ambas variables cualitativas es independiente (modelo aditivo), que debe utilizarse si se estima que no existe interacción entre esas variables. Si esas dos variables podrían interactuar creando un efecto sinérgico, debe reemplazarse el símbolo más (+) por un asterisco (*)

prueba_aov3 <- aov(edad ~ uso_movil_let * genero, 
                data = nuevos_datos)
summary(prueba_aov3)

##                      Df Sum Sq Mean Sq F value  Pr(>F)   
## uso_movil_let         2    1.4    0.68   0.131 0.87709   
## genero                1   39.4   39.37   7.559 0.00716 **
## uso_movil_let:genero  2   11.1    5.54   1.063 0.34962   
## Residuals            94  489.6    5.21                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A partir de los resultados del cuadro superior, se puede concluir que, a partir de los valores p y para un nivel de significación de 0,05:

• el valor p de uso_movil_let es 0.87709 (no significativo), lo que indica que los niveles de esta variable no están asociados al peso. significativamente diferentes.

• el valor p de genero es 0.00716 (significativo), lo que indica que los niveles de esta variable sí están asociados con el peso.

• el valor p para la interacción entre uso_movil_let y genero es 0.34962 (no significativo), lo que indica que las relaciones entre el peso y el uso de móvil y el género no interactúan.

Cálculo de algunos estadísticos de resumen

También es posible usar la función model.tables() de la siguiente manera:

model.tables(prueba_aov3, type="means", se = TRUE)

## Tables of means
## Grand mean
##       
## 21.31 
## 
##  uso_movil_let 
##     Escaso Moderado Intenso
##      21.41    21.13   21.36
## rep  37.00    30.00   33.00
## 
##  genero 
##     Femenino Masculino
##        21.84     20.57
## rep    58.00     42.00
## 
##  uso_movil_let:genero 
##              genero
## uso_movil_let Femenino Masculino
##      Escaso   22.048   20.563   
##      rep      21.000   16.000   
##      Moderado 22.000   20.000   
##      rep      17.000   13.000   
##      Intenso  21.500   21.154   
##      rep      20.000   13.000

Obsérvese que el programa envía un aviso señalando que Design is unbalanced - use se.contrast() for se's.

Comparación múltiple por pares entre las medias de los grupos.

En la prueba ANOVA, un valor p significativo indica que algunas de las medias de los grupos son diferentes, pero no sabemos qué pares de grupos son diferentes. Es posible realizar múltiples comparaciones por pares para determinar si la diferencia media entre pares específicos de grupos es estadísticamente significativa.

** Método de Tukey para realizar múltiples comparaciones por pares**

Como la prueba ANOVA es significativa, podemos calcular el estadístico HSD de Tukey (Diferencias significativas honestas) con la función TukeyHSD() para realizar múltiples comparaciones por pares entre las medias de los grupos. Esta función toma el ANOVA ajustado como argumento. No necesitamos realizar la prueba para la variable “soporte” porque tiene solo dos niveles, que ya se ha demostrado que son significativamente diferentes mediante la prueba ANOVA. Por lo tanto, la prueba Tukey HSD se realizará solo para la variable de factor “genero”.

TukeyHSD(prueba_aov3, 
         which = "genero")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = edad ~ uso_movil_let * genero, data = nuevos_datos)
## 
## $genero
##                         diff       lwr        upr     p adj
## Masculino-Femenino -1.270406 -2.188496 -0.3523161 0.0071995

• diff: diferencia entre las medias de los dos grupos

• lwr,upr: los valores inferior y superior del intervalo de confianza al 95 % (predeterminado).

• p adj: valor de p después del ajuste para las comparaciones múltiples.

Puede verse en la salida que todas las comparaciones por pares son significativas con un valor p ajustado < 0,05.

Comparaciones múltiples utilizando el paquete multcomp

La función glht() (paquete multcomp) realiza múltiples procedimientos de comparación para un ANOVA. El formato simplificado es el siguiente:

• glht(modelo, lincft)

• modelo: un modelo ajustado, por ejemplo, un objeto devuelto por aov().

• lincft(): una especificación de las hipótesis lineales a probar. Las comparaciones múltiples en los modelos ANOVA se especifican mediante objetos devueltos por la función mcp().

library(multcomp)

summary(glht(prueba_aov3, 
             linfct = mcp(uso_movil_let = "Tukey")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = edad ~ uso_movil_let * genero, data = nuevos_datos)
## 
## Linear Hypotheses:
##                         Estimate Std. Error t value Pr(>|t|)
## Moderado - Escaso == 0  -0.04762    0.74457  -0.064    0.998
## Intenso - Escaso == 0   -0.54762    0.71304  -0.768    0.723
## Intenso - Moderado == 0 -0.50000    0.75285  -0.664    0.785
## (Adjusted p values reported -- single-step method)

Prueba t por pares

La función pairwise.t.test() también se puede utilizar para calcular comparaciones por pares entre niveles de grupo con correcciones para pruebas múltiples.

pairwise.t.test(nuevos_datos$edad, 
                nuevos_datos$genero,
                p.adjust.method = "BH")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  nuevos_datos$edad and nuevos_datos$genero 
## 
##           Femenino
## Masculino 0.0066  
## 
## P value adjustment method: BH

Comprobando los requisitos de ANOVA

ANOVA asume que los datos se distribuyen normalmente y que la varianza entre los grupos es homogénea. Podemos verificar eso con algunas herramientas de diagnóstico.

Comprobando la homogeneidad de varianzas

La gráfica de residuos versus valores ajustados se utiliza para comprobar la homogeneidad de las varianzas. En el gráfico a continuación, no hay relaciones evidentes entre los residuos y los valores ajustados (la media de cada grupo), lo cual permite suponer la homogeneidad de las varianzas.

plot(prueba_aov3, 1)

Los puntos 344, 346 y 347 se detectan como valores atípicos, lo que puede afectar gravemente la normalidad y la homogeneidad de la varianza. Puede ser útil eliminar los valores atípicos para cumplir con los supuestos de la prueba.

Use la prueba de Levene para verificar la homogeneidad de las varianzas. Se utilizará la función leveneTest()(paquete cars).

library(car)

leveneTest(edad ~ uso_movil_let * genero, 
           data = nuevos_datos)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  0.8809 0.4972
##       94

Del resultado anterior podemos ver que el valor p no es menor que el nivel de significancia de 0.05. Esto significa que no hay evidencia que sugiera que la varianza entre los grupos sea significativamente diferente desde el punto de vista estadístico. Por lo tanto, podemos asumir la homogeneidad de las varianzas en los diferentes grupos de tratamiento.

Comprobación de la condición de normalidad

En el siguiente gráfico, los cuantiles de los residuos se representan frente a los cuantiles de la distribución normal. También se traza una línea de referencia de 45 grados. La gráfica de probabilidad normal de residuos se utiliza para verificar la suposición de que los residuos se distribuyen normalmente, que deben seguir aproximadamente una línea recta.

plot(prueba_aov3, 2)

Como no todos los puntos se sitúan aproximadamente a lo largo de esta línea de referencia, no se cumple la condición de normalidad. La conclusión anterior está respaldada por la prueba de Shapiro-Wilk sobre los residuos de ANOVA (W = 0,98, p = 0,5) que ofrece indicios de que se haya violado la normalidad.

prueba_aov3_residuals <- residuals(object = prueba_aov3)                  # Extrae los resíduos
shapiro.test(x = prueba_aov3_residuals)                               # Ejecuta la prueba de Shapiro-Wilk sobre los residuos

## 
##  Shapiro-Wilk normality test
## 
## data:  prueba_aov3_residuals
## W = 0.97015, p-value = 0.02268

Cálculo de la prueba two-way ANOVA para diseños no equilibrados

Un diseño desequilibrado tiene un número desigual de sujetos en cada grupo. Hay tres formas fundamentalmente diferentes de ejecutar un ANOVA bajo estas condiciones y se conocen como sumas de cuadrados tipo I, tipo II y tipo III. Para simplificar las cosas, el método recomendado son las sumas de cuadrados de tipo III. Los tres métodos dan el mismo resultado cuando el diseño está equilibrado, pero en condiciones de desequilibrio no proporcionan los mismos resultados. La función Anova() del paquete car se puede usar para calcular la prueba ANOVA de dos vías para diseños no equilibrados.

library(car)
prueba_aov3_bis <- aov(edad ~ genero * uso_movil_let, 
                 data = nuevos_datos)
Anova(prueba_aov3_bis, type = "III")

## Anova Table (Type III tests)
## 
## Response: edad
##                       Sum Sq Df   F value    Pr(>F)    
## (Intercept)          27463.7  1 5273.0407 < 2.2e-16 ***
## genero                  39.4  1    7.5601  0.007157 ** 
## uso_movil_let            3.7  2    0.3505  0.705250    
## genero:uso_movil_let    11.1  2    1.0628  0.349618    
## Residuals              489.6 94                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

MANOVA Test: Multivariate Analysis of Variance

Por último, cuando hay múltiples variables (cuantitativas) de respuesta, puede probarlas simultáneamente usando un análisis de varianza multivariante (MANOVA).

Por ejemplo, en el siguiente ejercicio seleccionaremos 3 variables, 1 cualitativa (genero) y dos cuantitativas (peso y edad).

nuevos_datos2 <- subset(alumnos_uc, 
                       select=c("genero","peso", "edad")) 

MANOVA se utilizar bajo ciertos supuestos:

• Las variables dependientes deben distribuirse normalmente dentro de los grupos.

• Homogeneidad de las varianzas en toda la gama de predictores.

• Linealidad entre todos los pares de variables dependientes, todos los pares de covariables y todos los pares de variables dependientes-covariables en cada celda

Interpretación de MANOVA

Si la prueba multivariante global es significativa, concluimos que el efecto correspondiente (tratamiento) es significativo. En ese caso, la siguiente pregunta es determinar si el tratamiento afecta solo el peso, solo la altura o ambos. En otras palabras, queremos identificar las variables dependientes específicas que contribuyeron al efecto global significativo.

prueba_manova <- manova(cbind(peso, edad) ~ genero, 
                     data = nuevos_datos2)
summary(prueba_manova)

##           Df   Pillai approx F num Df den Df  Pr(>F)  
## genero     1 0.078001   4.1031      2     97 0.01947 *
## Residuals 98                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos verificar en qué medida difieren

summary.aov(prueba_manova)

##  Response peso :
##             Df  Sum Sq Mean Sq F value Pr(>F)
## genero       1    51.3  51.341  0.4458 0.5059
## Residuals   98 11285.4 115.158               
## 
##  Response edad :
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## genero       1  39.50  39.501   7.713 0.00657 **
## Residuals   98 501.89   5.121                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Del resultado anterior, se puede ver que las dos variables tienen un peso diferentes en la variable respuesta.

COMPARACIÓN DE VARIANZAS

Comparación de la varianza de dos poblaciones.

Como hemos señalado en el apartado 2.1.4, la prueba F se utiliza para:

• Verificar la igualdad de las varianzas de las dos muestras.

• Comparar la variabilidad de un nuevo método de medición con uno antiguo. ¿El nuevo método reduce la variabilidad de la medida?

Las preguntas típicas de investigación son:

• si la varianza del grupo A (σ2A) es igual a la varianza del grupo B (σ2B)?

• si la varianza del grupo A (σ2A) es menor que la varianza del grupo B (σ2B)?

• si la varianza del grupo A (σ2A) es mayor que la varianza del grupo B (σ2B)?

De acuerdo con lo señalado anteriormente, las hipótesis a verificar son

• Hipótesis nula H0: las varianzas de ambas variables son iguales. Podemos definir la hipótesis nula correspondiente (H0) de la siguiente manera:

  • H0:σ2A=σ2B

  • H0:σ2A≤σ2B

  • H0:σ2A≥σ2B

• Hipótesis alternativa H1: las varianzas de ambas variables no son iguales. Las hipótesis alternativas correspondientes son las siguientes:

  • Ha:σ2A≠σ2B (diferente)

  • Ha:σ2A>σ2B (mayor)

  • Ha:σ2A<σ2B (menos)

De acuerdo con ello

• Si el valor de p es inferior a .05, se rechaza la hipótesis nula/aceptamos la alternativa: las varianzas de ambas variables no son igulaes.

• Si el valor de p es superior a .05, aceptamos la hipótesis nula: las varianzas de ambas variables son iguales.

Tenga en cuenta que:

• Las hipótesis 1) se denominan pruebas de dos colas

• Las hipótesis 2) y 3) se denominan pruebas de una cola

La estadística de prueba se puede obtener calculando la razón de las dos varianzas S2A y S2B.

F=S2AS2B

Los grados de libertad son nA−1 (para el numerador) y nB−1 (para el denominador).

Tenga en cuenta que, cuanto más se desvíe esta relación de 1, más fuerte será la evidencia de varianzas poblacionales desiguales. Tenga en cuenta que la prueba F requiere que las dos muestras se distribuyan normalmente. La prueba F es muy sensible a la desviación del supuesto de normalidad, por lo que debe verificarse si los datos se distribuyen normalmente antes de usarla. La prueba de Shapiro-Wilk se puede utilizar para comprobar si se cumple esta suposición, al igual que el gráfico Q-Q (gráfico cuantil-cuantil).

Si hay dudas sobre la normalidad, la mejor opción es utilizar la prueba de Levene o la prueba de Fligner-Killeen, menos sensibles a la desviación del supuesto normal.

La función de R var.test() se puede usar para comparar

var.test(peso ~ genero, 
        data = alumnos_uc,
        alternative = "two.sided")        # Hipótesis alternativa: “two.sided” (por defecto), “greater” o “less”)

## 
##  F test to compare two variances
## 
## data:  peso by genero
## F = 1.0953, num df = 57, denom df = 41, p-value = 0.7672
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6082483 1.9193511
## sample estimates:
## ratio of variances 
##           1.095283

La función var.test() devuelve una lista que contiene los siguientes componentes:

• statistic: el valor del estadístico de la prueba F.

• parameter: los grados de libertad de la distribución F del estadístico de prueba.

• p.value: el valor p de la prueba.

• conf.int: un intervalo de confianza para el cociente de las varianzas de la población.

• estimar: el cociente de las varianzas de la muestra

Para conseguir estos valores podemos escribir:

res.ftest$estimate

## ratio of variances 
##           1.095283

res.ftest$p.value

## [1] 0.7672127

El valor p de la prueba F es p = 1.078e-05 es menor que el nivel de significancia de 0,05. En conclusión, existe una diferencia significativa entre las dos varianzas.

Comparación de las varianzas de múltiples muestras

Existen diferentes opciones para probar la igualdad (homogeneidad) de la varianza entre grupos, que incluyen:

• La prueba F: compara las varianzas de dos muestras.

• La prueba de Bartlett: compara las varianzas de \(k\) muestras, donde \(k\) puede ser más de dos muestras.

• La prueba de Levene: compara las varianzas de \(k\) muestras, donde \(k\) puede ser más de dos muestras. Es una alternativa a la prueba de Bartlett que es menos sensible a las desviaciones de la normalidad.

• La prueba de Fligner-Killeen: prueba no paramétrica muy robusta frente a desviaciones de la normalidad.

Para todas estas pruebas es posible formular las siguientes pruebas estadísticas:

• La hipótesis nula es que todas las varianzas de las poblaciones son iguales.

• La hipótesis alternativa es que al menos dos de ellos difieren.

Prueba de Bartlett

La prueba de Bartlett se usa para probar la homogeneidad de las varianzas en varias muestras, y requiere que los datos sigan una distribución normal. La función bartlett.test() tiene la siguiente sintaxis:

bartlett.test(formula, data)

• formula: una sintaxis en la forma variable cuantitativa ~ variable cualitativa.

• data: matriz o dataframe a analizar.

La función devuelve una lista que contiene la siguiente información:

• statistic: estadístico de Bartlett.

• parameter: grados de libertad de la distribución de chi-cuadrado.

• p.value: el valor de p de la prueba.

Prueba de Bartlett con una variable independiente

bartlett.test(peso ~ genero, 
        data = alumnos_uc)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  peso by genero
## Bartlett's K-squared = 0.097218, df = 1, p-value = 0.7552

A partir de la salida de esta función, el valor p de 2.2e-16 es inferior al nivel de significancia de 0,05. Esto significa que hay evidencias que sugiere que el tamaño de las viviendas está relacionado con su accesiblidad desde el punto de vista estadístico.

Prueba de Bartlett con múltiples variables independientes: la función interaction() debe usarse para colapsar múltiples factores en una sola variable que contiene todas las combinaciones de los factores.

bartlett.test(peso ~ interaction(genero, procedencia),
        data = alumnos_uc)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  peso by interaction(genero, procedencia)
## Bartlett's K-squared = 11.556, df = 11, p-value = 0.3979

Prueba de Levene

Es una alternativa a la prueba de Bartlett cuando los datos no se distribuyen normalmente. Se puede utilizar la función leveneTest() (en el paquete car).

library(car)

Prueba de Levene con una variable independiente

leveneTest(peso ~ genero, 
           data = alumnos_uc)                     

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.7408 0.3915
##       98

Prueba de Leven con múltiples variables independientes

leveneTest(peso ~ genero*procedencia, data = alumnos_uc)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 11  1.5204 0.1384
##       88

Prueba de Fligner-Killeen

Es una de las muchas pruebas de homogeneidad de varianzas que es más robusta frente a las desviaciones de la normalidad. La función fligner.test() :

fligner.test(peso ~ genero, 
           data = alumnos_uc)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  peso by genero
## Fligner-Killeen:med chi-squared = 0.48995, df = 1, p-value = 0.4839

rm(res, res.ftest)