Pruebas de hipótesis
En este ejemplo veremos cómo aplicar pruebas de hipótesis usando R. Para ello usaremos el conjunto de datos cars de R y nos centraremos en la variable speed que representa la velocidad en mph.
Supongamos que se sabe que la variable sigue una distribución Normal y queremos probar la hipótesis
\[ H_0: \mu=0 \text{ vs. } H_A:\mu \neq 0 \] Primero cargamos los datos y usamos la función t.test
t.test(cars$speed)
##
## One Sample t-test
##
## data: cars$speed
## t = 20.594, df = 49, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 13.89727 16.90273
## sample estimates:
## mean of x
## 15.4
Observamos el p-valor y conclimos que se rehaza \(H_0\). ¿Qué otro tipo de hipótesis podemos probar?
t.test(cars$speed, mu=17)
##
## One Sample t-test
##
## data: cars$speed
## t = -2.1397, df = 49, p-value = 0.03739
## alternative hypothesis: true mean is not equal to 17
## 95 percent confidence interval:
## 13.89727 16.90273
## sample estimates:
## mean of x
## 15.4
t.test(cars$speed, mu=17, alternative = "less")
##
## One Sample t-test
##
## data: cars$speed
## t = -2.1397, df = 49, p-value = 0.01869
## alternative hypothesis: true mean is less than 17
## 95 percent confidence interval:
## -Inf 16.6537
## sample estimates:
## mean of x
## 15.4
¿Qué concluiría en cada uno de los casos anteriores?
Una forma de visualizar el p-valor, o equivalente, las regiones de aceptación y rechazo es la siguiente:
library(gginference)
ggttest(t.test(cars$speed, mu = 17, alternative = "less"))
## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
ggttest(t.test(cars$speed, mu = 17))
## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
Dos personas miden el tiempo que hacen de su casa a la escuela y quieren saber si hay diferencia. En este caso no se sabe nada acerca de la distribución de los datos, así que podemos usar una prueba:
datos<-read.csv("tiempos.csv")
#Prueba de normalidad
shapiro.test(datos$tiempo[1:10])
##
## Shapiro-Wilk normality test
##
## data: datos$tiempo[1:10]
## W = 0.9209, p-value = 0.3645
shapiro.test(datos$tiempo[11:20])
##
## Shapiro-Wilk normality test
##
## data: datos$tiempo[11:20]
## W = 0.91621, p-value = 0.3264
Con el p-valor podemos afirmar que hay evidencias suficientes para afirmar que los datos siguen una distribución normal. Ahora veremos qué pasa con la varianza:
var.test(tiempo~persona, datos)
##
## F test to compare two variances
##
## data: tiempo by persona
## F = 1.053, num df = 9, denom df = 9, p-value = 0.9399
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.2615592 4.2395166
## sample estimates:
## ratio of variances
## 1.053036
Y finalmente hacemos la prueba sobre la media
t.test(datos$tiempo[1:10]-datos$tiempo[11:20], altervative="g", mu=12)
##
## One Sample t-test
##
## data: datos$tiempo[1:10] - datos$tiempo[11:20]
## t = -11.78, df = 9, p-value = 9.011e-07
## alternative hypothesis: true mean is not equal to 12
## 95 percent confidence interval:
## 0.1202329 3.9477671
## sample estimates:
## mean of x
## 2.034
Ahora usaremos la base diamonds que contiene datos y atributos de 54,000 diamantes. Queremos saber si existen evidencias suficientes para afirmar que el precio de los diamantes de corte premium, claridad I1 entre 0.7 y 1.4 g y el precio de los diamantes de corte ideal con las mismas características cuestan aproximadamente lo mismo. Primero extraemos los datos que nos interesan
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.2 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
ideal<-diamonds%>%filter(cut=="Ideal"& clarity=="I1" & carat<1.4 & carat>0.7)
premium<-diamonds%>%filter(cut=="Premium"& clarity=="I1" & carat<1.4 & carat>0.7)
Veamos graficamente si es posible inferir algo
hist(ideal$price, col=rgb(0,0,1,1/4))
hist(premium$price, col=rgb(1,0,0,1/4),add=T)
Ahora una prueba de hipótesis formal
t.test(ideal$price,premium$price)
##
## Welch Two Sample t-test
##
## data: ideal$price and premium$price
## t = 6.2695, df = 206.6, p-value = 2.081e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 420.491 806.259
## sample estimates:
## mean of x mean of y
## 3526.542 2913.167
ggttest(t.test(ideal$price,premium$price))
## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
¿Qué concluiría en este ejemplo?
Ahora queremos averiguar si existen evidencias suficientes para afirmar que la proporción de diamantes calidad “very good” D es la misma que la proporción de diamantes “ideal” D. Primero extraemos los datos
diamonds%>%filter(cut=="Ideal"& color=="D")
## # A tibble: 2,834 x 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.3 Ideal D SI1 62.5 57 552 4.29 4.32 2.69
## 2 0.3 Ideal D SI1 62.1 56 552 4.3 4.33 2.68
## 3 0.71 Ideal D SI2 62.3 56 2762 5.73 5.69 3.56
## 4 0.71 Ideal D SI1 61.9 59 2764 5.69 5.72 3.53
## 5 0.71 Ideal D SI2 61.6 55 2767 5.74 5.76 3.54
## 6 0.76 Ideal D SI2 62.4 57 2770 5.78 5.83 3.62
## 7 0.73 Ideal D SI2 59.9 57 2770 5.92 5.89 3.54
## 8 0.75 Ideal D SI2 61.3 56 2773 5.85 5.89 3.6
## 9 0.72 Ideal D SI1 60.8 57 2782 5.76 5.75 3.5
## 10 0.64 Ideal D VS1 61.5 56 2787 5.54 5.55 3.41
## # … with 2,824 more rows
# 21551/ 2834
diamonds%>%filter(cut=="Very Good"& color=="D")
## # A tibble: 1,513 x 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Very Good D VS2 60.5 61 357 3.96 3.97 2.4
## 2 0.23 Very Good D VS1 61.9 58 402 3.92 3.96 2.44
## 3 0.26 Very Good D VS2 60.8 59 403 4.13 4.16 2.52
## 4 0.24 Very Good D VVS1 61.5 60 553 3.97 4 2.45
## 5 0.26 Very Good D VVS2 62.4 54 554 4.08 4.13 2.56
## 6 0.26 Very Good D VVS2 62.8 60 554 4.01 4.05 2.53
## 7 0.26 Very Good D VVS1 62.1 60 554 4.03 4.12 2.53
## 8 0.75 Very Good D SI1 63.2 56 2760 5.8 5.75 3.65
## 9 0.61 Very Good D VVS2 59.6 57 2763 5.56 5.58 3.32
## 10 0.71 Very Good D SI1 63.6 58 2764 5.64 5.68 3.6
## # … with 1,503 more rows
#12,082/1,513
#?prop.test
prop.test(c(1513,2834),c(12082,21551),alternative="two.sided")
##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(1513, 2834) out of c(12082, 21551)
## X-squared = 2.6527, df = 1, p-value = 0.1034
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.013767833 0.001219019
## sample estimates:
## prop 1 prop 2
## 0.1252276 0.1315020
Si quisieramos una prueba de una sola muestra
prop.test(x=1513,n=12082,p=0.12,alternative="two.sided")
##
## 1-sample proportions test with continuity correction
##
## data: 1513 out of 12082, null probability 0.12
## X-squared = 3.0774, df = 1, p-value = 0.07939
## alternative hypothesis: true p is not equal to 0.12
## 95 percent confidence interval:
## 0.1194042 0.1312909
## sample estimates:
## p
## 0.1252276