Ejercicios - Pruebas de Contraste de Hipótesis para dos Muestras

logo

1. Enunciado

Para esta actividad se usará el conjunto de datos airquality. Este contiene 153 mediciones diarias de la calidad del aire en Nueva York, de mayo a septiembre de 1973, incluyendo las siguientes 6 variables:

[,1] Ozone numeric Ozone (ppb). Ozono medio en partes por billón de 1300 a 1500 horas en Roosevelt Island.
[,2] Solar.R numeric Solar R (lang). Radiación solar en Langleys en la banda de frecuencia 4000–7700 Angstroms de 08:00 a 12:00 horas en Central Park.
[,3] Wind numeric Wind (mph). Velocidad promedio del viento en millas por hora a las 07:00 y 10:00 horas en el Aeropuerto La Guardia.
[,4] Temp numeric Temperature (degrees F). Temperatura máxima diaria en grados Fahrenheit en el Aeropuerto La Guardia.
[,5] Month numeric Month (1–12)
[,6] Day numeric Day of month (1–31)

Se nos pide estudiar si la radiación solar es diferente en los meses de invierno respecto a los de verano.

2. Desarrollo

2.1 Estadística descriptiva

En el primer paso estudiaremos el conjuntos de datos desde el punto de vista de la Estadística Descriptiva, para ello construiremos un gráfico de caja, obteniendo algunas medidas de centralidad y dispersión.

#Preparación de datos
datos=airquality
period=rep("Verano",length(datos$Month))
period[which(datos$Month==5)]="Invierno"
period[which(datos$Month==6 & datos$Day<21)]="Invierno"
period[which(datos$Month==9 & datos$Day>23)]="Invierno"
datos$period = period

head(datos)

##   Ozone Solar.R Wind Temp Month Day   period
## 1    41     190  7.4   67     5   1 Invierno
## 2    36     118  8.0   72     5   2 Invierno
## 3    12     149 12.6   74     5   3 Invierno
## 4    18     313 11.5   62     5   4 Invierno
## 5    NA      NA 14.3   56     5   5 Invierno
## 6    28      NA 14.9   66     5   6 Invierno

#Construcción de Gráfico de Caja
library("ggplot2")
g=ggplot(datos, aes(period,Solar.R)) +  geom_boxplot(fill = "paleturquoise", color = "cadetblue4") + 
  labs(x="Estación del año", y="Radiación solar en Langleys") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(text = element_text(size = 12)) 
plot(g)

#Medidas estadísticas
library("psych")
estadisticos=describeBy(datos$Solar.R, datos$period, mat = F)
print(estadisticos)

## 
##  Descriptive statistics by group 
## group: Invierno
##    vars  n   mean     sd median trimmed    mad min max range  skew kurtosis
## X1    1 54 193.07 102.36    207  197.91 118.61   8 334   326 -0.36    -1.25
##       se
## X1 13.93
## ------------------------------------------------------------ 
## group: Verano
##    vars  n   mean    sd median trimmed   mad min max range  skew kurtosis   se
## X1    1 92 181.74 82.29    205  187.64 80.06   7 314   307 -0.55    -0.95 8.58

2.2 Evaluación de supuestos paramétricos

2.2.1 Evaluación de supuesto de normalidad

#Normalidad
library("nortest")
# Conjunto completo
t1a=lillie.test(datos$Solar.R)
t1b=shapiro.test(datos$Solar.R)
print(t1a)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  datos$Solar.R
## D = 0.12, p-value = 2.35e-05

print(t1b)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos$Solar.R
## W = 0.94183, p-value = 9.492e-06

# Verano
t1a=lillie.test(datos$Solar.R[which(datos$period=="Verano")])
t1b=shapiro.test(datos$Solar.R[which(datos$period=="Verano")])
print(t1a)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  datos$Solar.R[which(datos$period == "Verano")]
## D = 0.13553, p-value = 0.0002593

print(t1b)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos$Solar.R[which(datos$period == "Verano")]
## W = 0.92505, p-value = 5.407e-05

#Invierno
t1a=lillie.test(datos$Solar.R[which(datos$period=="Invierno")])
t1b=shapiro.test(datos$Solar.R[which(datos$period=="Invierno")])
print(t1a)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  datos$Solar.R[which(datos$period == "Invierno")]
## D = 0.13686, p-value = 0.01321

print(t1b)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos$Solar.R[which(datos$period == "Invierno")]
## W = 0.92258, p-value = 0.001874

Las pruebas de contraste indican que la distribución indicada no es normal. Podemos complementar el análisis con un gráfico QQ.

#QQplot
qqnorm(datos$Solar.R, pch = 19, col = "gray50")
qqline(datos$Solar.R)

2.2.2 Evaluación de homocedasticidad

Podemos aplicar la prueba F para evaluar el principio de homocedasticidad. En este caso ambos conjuntos tienen varianza homogénea.

#Prueba F
t1 = var.test(datos$Solar.R[which(datos$period=="Invierno")],datos$Solar.R[which(datos$period=="Verano")],conf.level=0.95)
print(t1)

## 
##  F test to compare two variances
## 
## data:  datos$Solar.R[which(datos$period == "Invierno")] and datos$Solar.R[which(datos$period == "Verano")]
## F = 1.5472, num df = 53, denom df = 91, p-value = 0.06739
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.9692546 2.5479263
## sample estimates:
## ratio of variances 
##            1.54722

2.2.3 Prueba de contraste de hipótesis para dos muestras.

Como los datos NO siguen una distribución normal, se utilizará una prueba de contraste de hipótesis para dos muestras NO paramétrico. En este caso, la prueba Wilcoxon rank-sum, considerando que las muestras son no pareadas.

#Prueba Wilcoxon rank-sum

res = wilcox.test(Solar.R ~ period, data = datos,conf.level = 0.95,
                   exact = FALSE,paired = F)
print(res)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Solar.R by period
## W = 2765.5, p-value = 0.2547
## alternative hypothesis: true location shift is not equal to 0

Los resultados indican que no hay diferencia en la radiación entre invierno y verano considerando un nivel de confianza del 95% (p<0.05).

Actividad
Buscar y resolver un ejemplo similar al estudiado.