Inferencia sobre dos poblaciones

Autor/a

Edimer David Jaramillo

Fecha de publicación

30 de mayo de 2025

Bibliotecas

Código 
library(tidyverse)
library(readxl)
library(infer)

theme_set(theme_bw())

Inferencia sobre 2 poblaciones

Inferencia sobre \(\sigma^2\)

Datos

Código 
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:
Código 
library(ggpubr)
ggqqplot(data = datos_cebada$yield)

  • Prueba de shapiro wilk:
Código 
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)
Código 
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.
Código 
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:
Código 
library(car)
leveneTest(datos_cebada$yield ~ datos_cebada$year)

Inferencia sobre \(\mu_1 - \mu_2\)

Datos

Código 
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.
Código 
nitro_bajo <- datos_trigo |> filter(nitro == "L")
nitro_alto <- datos_trigo |> filter(nitro == "H")

ggqqplot(nitro_bajo$yield)

Código 
ggqqplot(nitro_alto$yield)

  • Pruebas de shapiro wilk:
Código 
shapiro.test(nitro_alto$yield)

    Shapiro-Wilk normality test

data:  nitro_alto$yield
W = 0.90953, p-value = 2.064e-05
Código 
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.
Código 
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.
Código 
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

Datos

Código 
datos <- read_excel("../datos/datos-encuestas-historia.xlsx")
datos

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.
Código 
table(datos$trabaja)

No Sí 
71 22

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 = 22) y no trabajan (n = 71)
    • n: total de ensayos (observaciones). En este caso equivale a 93.
    • 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
Código 
prop.test(x = c(22, 71),
          n = c(93, 93),
          alternative = "two.sided",
          conf.level = 0.95)

    2-sample test for equality of proportions with continuity correction

data:  c(22, 71) out of c(93, 93)
X-squared = 49.548, df = 1, p-value = 1.935e-12
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.6597804 -0.3939831
sample estimates:
   prop 1    prop 2 
0.2365591 0.7634409
  • Conclusión:
    • Como el valor p (1.935e-12) es menor que el nivel de significancia (0.05), existe evidencia para rechazar la hipótesis nula, es decir, que la proporción de estudiantes que trabajan y no trabajan no es la misma.
    • Como el intervalo de confianza para la diferencia de proporciones ([LI: -0.6597804 , LS: -0.3939831]) no contiene al cero, entonces existe evidencia para rechazar la hipótesis nula.

T-student medias pareadas

Datos

Código 
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

Código 
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:
Código 
ggqqplot(data = diferencia)

  • Shapiro Wilk: la normalidad se cumple (0.725)
Código 
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)
Código 
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

Código 
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 con Bootstrapping

\(\mu_1 - \mu_2\)

Código 
set.seed(2025)
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:
Código 
ic_dif_medias <-
  remuestreo_dif_medias |> 
  get_confidence_interval(level = 0.95,
                          type = "percentile")

ic_dif_medias

Medias pareadas

Código 
set.seed(2025)
remuestreo_mu_pareadas <- calificaciones |>
  mutate(diff = Post - Pre) |> 
  specify(response = diff) |>
  generate(reps = 1000, type = "bootstrap") |> 
  calculate(stat = "mean")

remuestreo_mu_pareadas |> 
  visualize()

  • Calculamos el intervalo de confianza:
Código 
ic_dif_pareadas <-
  remuestreo_mu_pareadas |> 
  get_confidence_interval(level = 0.95,
                          type = "percentile")

ic_dif_pareadas