Inferencia estadística

Pruebas de hipótesis

Author

Edimer David Jaramillo

Published

November 14, 2024

Escenarios de inferencia

Bibliotecas

Code

library(readxl)
library(tidyverse)

Datos de encuesta

Code

datos <- read_excel("datos/encuesta_depurada.xlsx")
datos

Analizando la normalidad

Gráficos:
- Histogramas
- Densidades
- Gráficos cuantil-cuantil (Q-Q Norm)
Pruebas de hipótesis:
- Prueba de Shapiro Wilk
- Prueba de Anderson Darling

Gráficos

Densidades

Code

datos |> 
  ggplot(mapping = aes(x = promedio_academico)) +
  geom_density()

Cuantil-Cuantil

Code

datos |> 
  ggplot(mapping = aes(sample = promedio_academico)) +
  geom_qq() +
  geom_qq_line()

El mismo gráfico anterior con la biblioteca ggpubr:

Code

library(ggpubr)
ggqqplot(data = datos$promedio_academico)

El mismo gráfico con la biblioteca car:

Code

library(car)
qqPlot(x = datos$promedio_academico)

[1]  4 22

Pruebas de hipótesis

Juego de hipótesis:

\[H_0: X \sim N(\mu, \sigma)\] \[H_1: X \nsim N(\mu, \sigma)\]

Nivel de significancia: en este caso vamos a usar un nivel de significancia del 1%, es decir, del 0.01 (\(\alpha = 0.01\))

Shapiro Wilk

Code

shapiro.test(x = datos$promedio_academico)


    Shapiro-Wilk normality test

data:  datos$promedio_academico
W = 0.94241, p-value = 0.1056

Conclusión: como el valor p (0.1056) es mayor que el nivel de significancia (0.01) no existe evidencia para rechazar la hipótesis nula, es decir, que la variable aleatoria promedio_académico se distribuye de forma normal.

Anderson Darling

Code

library(nortest)
ad.test(x = datos$promedio_academico)


    Anderson-Darling normality test

data:  datos$promedio_academico
A = 0.60659, p-value = 0.1046

Conclusión: como el valor p (0.1046) es mayor que el nivel de significancia (0.01) no existe evidencia para rechazar la hipótesis nula, es decir, que la variable aleatoria promedio_académico se distribuye de forma normal.

Inferencia sobre una población

Inferencia sobre la media

Normalidad de la variable

Se comprobó previamente que la variable promedio_académico se distribuye de forma normal.

Juego de hipótesis

\[H_0: \mu = 3.5\] \[H_1: \mu \neq 3.5\]

Nivel de significancia

Vamos a utilizar un nivel de significancia del 5% (\(\alpha = 0.05\))

Calcular el estadístico observado

En este caso calculamos el estadístico T:

\[T = \frac{\bar{X} - \mu}{S/\sqrt{n}}\]

Code

x_barra <- mean(datos$promedio_academico)
mu_referencia <- 3.5
desviacion_muestral <- sd(datos$promedio_academico)
raiz_n <- sqrt(nrow(datos))

\[T = \frac{3.699 - 3.5}{0.2488643/5.477226} = 4.379768\]

Code

(x_barra - mu_referencia) / (desviacion_muestral / raiz_n)

[1] 4.379768

Región de rechazo

Tabla t-student
Podemos obtener los límites critícos con R:

Code

qt(p = 0.025, df = 29, lower.tail = TRUE)

[1] -2.04523

Code

qt(p = 0.025, df = 29, lower.tail = FALSE)

[1] 2.04523

Conclusión: como el estadístico calculado (4.379768) está dentro de las regiones de rechazo, existe evidencia para rechazar la hipótesis nula, es decir, que se rechaza que el promedio del promedio_académico sea igual a 3.5. Esto se concluye con un nivel de significancia del 5%.

Intervalo de confianza para \(\mu\)

Límite inferior del intervalo de confianza:

\[\bar{X} - t_{\alpha/2, n-1} \times \frac{s}{\sqrt{n}}\]

Code

x_barra - (2.045 * (desviacion_muestral / raiz_n))

[1] 3.606083

Límite superior del intervalo de confianza:

\[\bar{X} + t_{\alpha/2, n-1} \times \frac{s}{\sqrt{n}}\]

Code

x_barra + (2.045 * (desviacion_muestral / raiz_n))

[1] 3.791917

Conclusión: como el valor de referencia (3.5) no está dentro del intervalo de confianza, existe evidencia para rechazar la hipótesis nula, es decir, que se rechaza que el promedio del promedio_académico sea igual a 3.5. Esto se concluye con un nivel de significancia del 5%.

Valor p

Calcular el área que deja un valor de 4.371917 a la izquierda:

Code

pt(q = -4.371917, df = 29, lower.tail = TRUE)

[1] 7.228808e-05

Calcular el área que deja un valor de 4.371917 a la derecha:

Code

pt(q = 4.371917, df = 29, lower.tail = FALSE)

[1] 7.228808e-05

El valor p es la suma de las dor áreas anteriores:

Code

7.228808e-05 + 7.228808e-05

[1] 0.0001445762

Conclusión: como el valor p (0.0001445762) es menor que el nivel de significancia (0.05) existe evidencia para rechazar la hipótesis.

Solución con R

Utilizamos la función t.test() con los siguientes argumentos:
- x: la variable sobre la cual estamos haciendo inferencia. En este caso el promedio_académico
- alternative: tipo de hipótesis alternativa. En este es una prueba bilateral usamos “two.sided”
- conf.level: nivel de confianza (1 - nivel de significancia = 1 - 0.05 = 0.95)
- mu: valor promedio de referencia. En este caso es 3.5

Code

t.test(x = datos$promedio_academico,
       alternative = "two.sided",
       conf.level = 0.95,
       mu = 3.5)


    One Sample t-test

data:  datos$promedio_academico
t = 4.3798, df = 29, p-value = 0.0001415
alternative hypothesis: true mean is not equal to 3.5
95 percent confidence interval:
 3.606073 3.791927
sample estimates:
mean of x 
    3.699

Presentando resultados de la prueba anterior, de forma ordenada (tidy):

Code

prueba_t1 <- t.test(
  x = datos$promedio_academico,
  alternative = "two.sided",
  conf.level = 0.95,
  mu = 3.5
)

library(broom)
prueba_t1 |> tidy()

Inferencia sobre una proporción

Ejemplo: En este caso estamos interesados en la proporción de estudiantes que trabajan.

Juepo de hipótesis

\[H_0: p = 0.5\] \[H_1: p \neq 0.5\]

Nivel de significancia

En este caso usamos un nivel de significancia del 5%

Calcular el estadístico

Calculamos la proporción de los estudiantes que trabajan

Code

prop.table(table(datos$trabajo))[[2]]

[1] 0.3666667

Solución con R

Usamos la función prop.test() con los siguientes argumentos:
- x: el número de éxitos para el evento de interés. En este caso son estudiantes que trabajan (n = 11)
- n: total de ensayos (observaciones). En este caso equivale a 30.
- p: proporción de referencia. En este caso es 0.5
- alternative: tipo de hipótesis alternativa. En este caso es “two.sided”
- conf.level: nivel de confianza. En este caso es 1 - 0.05 = 0.95

Code

prop.test(
  x = 11,
  n = 30,
  p = 0.5,
  alternative = "two.sided",
  conf.level = 0.95
)


    1-sample proportions test with continuity correction

data:  11 out of 30, null probability 0.5
X-squared = 1.6333, df = 1, p-value = 0.2012
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.2054281 0.5609198
sample estimates:
        p 
0.3666667

Conclusión:
- como el valor p (0.2012) es mayor que el nivel de significancia (0.05), no existe evidencia para rechazar la hipótesis nula, es decir, que la proporción de estudiantes que trabajan es igual a 50% (0.5)
- como el intervalo de confianza (LI: 0.2054281, LS: 0.5609198) contiene el valor de referencia, no existe evidencia para rechazar la hipótesisis nula.

Inferencia sobre 2 poblaciones

Inferencia sobre \(\sigma^2\)

Datos

Code

datos_cebada <- read_excel("datos/datos_cebada.xlsx") |>
  mutate(year = as.factor(year))
datos_cebada

Hipótesis

\[H_0: \sigma^2_{1931} / \sigma^2_{1932} = 1\] \[H_1: \sigma^2_{1931} / \sigma^2_{1932} \neq 1\]

Nivel de significancia

En este caso vamos a usar un nivel de significancia del 5%.

Analizando la normalidad

Gráfico:

Code

ggqqplot(data = datos_cebada$yield)

Prueba de shapiro wilk:

Code

shapiro.test(x = datos_cebada$yield)


    Shapiro-Wilk normality test

data:  datos_cebada$yield
W = 0.96055, p-value = 0.05005

Solución con R

Vamos a usar la función var.test() con los siguientes argumentos:
- formula: y ~ x. En este caso “y” es la variable produccion y el “x” es el año.
- ratio: es el resultado del cociente de las dos varianzas. En este caso asumimos en la hipótesis nula el valor de “1”.
- alternative: tipo de prueba. En este es bilateral (“two.sided”)
- conf.level: nivel de confianza. En este caso es 0.95 (1 - 0.5)

Code

var.test(datos_cebada$yield ~ datos_cebada$year,
         ratio = 1,
         alternative = "two.sided",
         conf.level = 0.95)


    F test to compare two variances

data:  datos_cebada$yield by datos_cebada$year
F = 1.3952, num df = 29, denom df = 29, p-value = 0.375
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.6640874 2.9314037
sample estimates:
ratio of variances 
          1.395245

Conclusión:
- Como el valor p (0.375) es mayor que el nivel de significancia (0.05), no existe evidencia para rechazar la hipótesis, es decir, que la varianza para la producción del año 1931 es igual a la varianza para la producción del año 1932.
- Como el intervalo de confianza (LI: 0.6640874, LS: 2.9314037) contiene el valor de referencia (“1”), no existe evidencia para rechazar la hipótesis.

Prueba de Bartlett

Esta prueba es recomendada para comparar varianzas cuando la normalidad se cumple.

Code

bartlett.test(datos_cebada$yield ~ datos_cebada$year)


    Bartlett test of homogeneity of variances

data:  datos_cebada$yield by datos_cebada$year
Bartlett's K-squared = 0.78702, df = 1, p-value = 0.375

Prueba de Levene

Esta prueba es recomendada para comparar varianzas cuando la normalidad no se cumple:

Code

leveneTest(datos_cebada$yield ~ datos_cebada$year)

Inferencia sobre \(\mu_1 - \mu_2\)

Datos

Code

datos_trigo <- read_excel("datos/datos_trigo.xlsx")
datos_trigo

Hipótesis

\[H_0: \mu_{bajo} = \mu_{alto}\]

\[H1: \mu_{bajo} \neq \mu_{alto}\]

El juego de hipótesis anterior es equivalente al siguiente:

\[H_0: \mu_{bajo} - \mu_{alto} = 0\]

\[H1: \mu_{bajo} - \mu_{alto} \neq 0\]

Nivel de significancia

El nivel de significancia será del 5%.

Normalidad

Para la prueba t-student sobre dos poblaciones (medias) se debe garantizar que la variable (en cada grupo) se distribuye de forma normal.

Code

nitro_bajo <- datos_trigo %>% filter(nitro == "L")
nitro_alto <- datos_trigo %>% filter(nitro == "H")

ggqqplot(nitro_bajo$yield)

Code

ggqqplot(nitro_alto$yield)

Pruebas de shapiro wilk:

Code

shapiro.test(nitro_alto$yield)


    Shapiro-Wilk normality test

data:  nitro_alto$yield
W = 0.90953, p-value = 2.064e-05

Code

shapiro.test(nitro_bajo$yield)


    Shapiro-Wilk normality test

data:  nitro_bajo$yield
W = 0.9039, p-value = 1.153e-05

Igualdad de varianzas

Como el valor p ( 0.001958 ) es menor que el nivel de significancia, existe evidencia para manifestar que las varianzas son diferentes para el grupo de nitrógeno alto y bajo.

Code

leveneTest(datos_trigo$yield ~ datos_trigo$nitro)

Solución con R

Usamos la función t.test() con los siguientes argumentos:
- formula: y ~ x. En este caso “y” es la producción y “x” es el nivel de fertilización con nitrógeno (bajo, alto)
- alternative: tipo de hipótesis alternativa. En este caso es bilateral.
- conf.level: nivel de confianza. En este caso es del 95%
- var.equal: toma valores TRUE o FALSE para cuando las varianzas son iguales o diferentes, respectivamente.
Nota: los supuestos matemáticos necesarios no se cumplen, por tal motivo se debe optar por algún método no paramétrico.

Code

t.test(datos_trigo$yield ~ datos_trigo$nitro,
       alternative = "two.sided",
       conf.level = 0.95,
       var.equal = FALSE)


    Welch Two Sample t-test

data:  datos_trigo$yield by datos_trigo$nitro
t = 2.9051, df = 143.36, p-value = 0.004254
alternative hypothesis: true difference in means between group H and group L is not equal to 0
95 percent confidence interval:
 17.44086 91.70200
sample estimates:
mean in group H mean in group L 
       517.4762        462.9048

Conclusión:
- Como el valor p (0.004254) es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que el promedio de producción de bajo nitrógeno es diferente al promedio de producción de alto nitrógeno.
- Como el intervalo de confianza no contiene al cero, entonces existe evidencia para rechazar la hipótesis nula. Además, como el intervalo de confianza está a la derecha del cero, quiere decir que en promedio el grupo de alta dosificación estará por encima del promedio de baja dosificación en nitrógeno.

Inferencia sobre p1-p2

Hipótesis

\[H_0: p1 = p2\] \[H_1: p1 \neq p2\]

O de forma equivalente:

\[H_0: p1 - p2 = 0\]

\[H_1: p1 - p2 \neq 0\]

Nivel de significancia

En este caso usaremos un nivel de significancia del 5% (0.05)

Calculando proporciones

En este caso queremos comaparar la proporción de estudiantes que trabajan vs los que no trabajan.

Code

table(datos$trabajo)


No Sí 
19 11

Solución con R

Usamos la función prop.test() con los siguientes argumentos:
- x: el número de éxitos para el evento de interés. En este caso son estudiantes que trabajan (n = 11) y no trabajan (n = 19)
- n: total de ensayos (observaciones). En este caso equivale a 30.
- alternative: tipo de hipótesis alternativa. En este caso es “two.sided”
- conf.level: nivel de confianza. En este caso es 1 - 0.05 = 0.95

Code

prop.test(x = c(11, 19),
          n = c(30, 30),
          alternative = "two.sided",
          conf.level = 0.95)


    2-sample test for equality of proportions with continuity correction

data:  c(11, 19) out of c(30, 30)
X-squared = 3.2667, df = 1, p-value = 0.0707
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.5438677  0.0105344
sample estimates:
   prop 1    prop 2 
0.3666667 0.6333333

Conclusión:
- Como el valor p (0.0707) es mayor que el nivel de significancia (0.05), no existe evidencia para rechazar la hipótesis nula, es decir, que la proporción de estudiantes que trabajan y no trabajan es la misma.
- Como el intervalo de confianza para la diferencia de proporciones ([LI: -0.5438677 , LS: 0.0105344]) contiene al cero, entonces no existe evidencia para rechazar la hipótesis nula.

T-student medias pareadas

Datos

Code

calificaciones <- read_excel("datos/datos_parciales.xlsx")
calificaciones

Juego de hipótesis

\[H_0: \mu_{post} - \mu_{pre} = 0\]

\[H_1: \mu_{post} - \mu_{pre} \neq 0\]

Nivel de significancia

En este caso vamos a usar un nivel de significancia del 5% (0.05)

Diferencia de medias

Code

diferencia <- calificaciones$Post - calificaciones$Pre
diferencia

 [1]  4  4  1  2 -3  5  3  2 -4  2  1 -1  2  7  0  4  6  3  4 -1

Normalidad

Gráfico:

Code

ggqqplot(data = diferencia)

Shapiro Wilk: la normalidad se cumple (0.725)

Code

shapiro.test(x = diferencia)


    Shapiro-Wilk normality test

data:  diferencia
W = 0.9686, p-value = 0.725

Igualdad de varianzas

Existe igualdad de varianzas (0.2795)

Code

var.test(x = calificaciones$Pre,
         y = calificaciones$Post,
         ratio = 1,
         alternative = "two.sided",
         conf.level = 0.95)


    F test to compare two variances

data:  calificaciones$Pre and calificaciones$Post
F = 0.60329, num df = 19, denom df = 19, p-value = 0.2795
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.238790 1.524186
sample estimates:
ratio of variances 
         0.6032913

Prueba t-student

Code

t.test(x = calificaciones$Pre,
       y = calificaciones$Post,
       alternative = "two.sided",
       conf.level = 0.95,
       paired = TRUE,
       var.equal = TRUE)


    Paired t-test

data:  calificaciones$Pre and calificaciones$Post
t = -3.2313, df = 19, p-value = 0.004395
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -3.3778749 -0.7221251
sample estimates:
mean difference 
          -2.05

Conclusión:
- Como el valor p (0.004395) es menor que el nivel de significancia, existe evidencia para rechazar la hipótesis nula, es decir, que el promedio del “post” es diferente al promedio para el “pre”.
- Como el intervalo de confianza está a la izquierda del cero ([-3.3778749, -0.7221251]), existe evidencia para rechazar la hipótesis nula. Además, el promedio del “pre” podría estar desde 0.72 hasta 3.37 puntos por debajo del “post”.

Inferencia sobre \(\mu_1 - \mu_2\) - Método no paramétrico

Code

library(infer)

set.seed(2024)
remuestreo_dif_medias <- datos_trigo |>
  specify(formula = yield ~ nitro) |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "diff in means", order = c("H", "L"))

remuestreo_dif_medias |> 
  visualize()

Calculamos el intervalo de confianza:

Code

ic_dif_medias <-
  remuestreo_dif_medias |> 
  get_confidence_interval(level = 0.95,
                          type = "percentile")

ic_dif_medias