Este documento explica el paso a paso cómo realizar un test de hipótesis en R.
Para ello vamos a usar una base de datos de pingüinos.

Pero antes, ¿cómo se construye un test de hipótesis en R?

Un test de hipótesis en R se construye a través de una función estadística que hace 3 cosas:
1. Recibe los datos y los agrupa según la hipótesis.
2. Calcula un estadístico.
3. Devuelve un objeto con los resultados: el valor del estadístico, el valor p, los grados de libertad, los intervalos de confianza, la hipótesis alternativa, otras informaciones específicas.

Por convención, los tests en R devuelven un objeto tipo htest. Y cada librería tiene su tipo de test. En R base (el R predeterminado) hay varias, entre ellas:
- t.test()
- chisq.test()
- fisher.test()
- wilcox.test()
- prop.test()
- var.test()

También está: - coin: que tiene tests exactos (wilcoxon, kruskal-wallis, etc.).
- rstatix: que tiene una sintaxis tidy para tests.
- statsr: bastante amigable en general.
- DescTools: tiene test adicionales como el test de la mediana y demás giladas (na, mentira).

En general, más o menos, la lógica es esta:
resultado <- x.test(variable1 ~ agrupamiento, data = data)

Vamos con el ejemplo de la función t.test().
t.test(x, y = NULL, alternative = c(“two.sided”, “less”, “greater”), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, …) t.test(formula, data, subset, na.action, …)

Argumentos

  • x: representa un vector numérico de valores de datos.
  • y: un vector numérico de datos.
  • alternative: una cadena de caracteres que especifica la hipótesis alternativa. Puede ser: "two.sided", "greater" o "less".
  • mu: un número que indica el valor de la media bajo la hipótesis nula.
  • paired: valor lógico; si es TRUE, se hace una prueba pareada.
  • var.equal: valor lógico; si es TRUE, se asumen varianzas iguales.
  • conf.level: el nivel de confianza para el intervalo (por defecto 0.95).
  • formula: fórmula tipo respuesta ~ grupo; respuesta es numérica y grupo un factor.
  • data: el data frame que contiene las variables.
  • subset: vector opcional para filtrar observaciones.
  • na.action: función para manejar valores faltantes (na.omit, na.exclude, etc.).
#Cargamos la base de datos
library(palmerpenguins)
data("penguins")

1 Un rato de exploración

Antes de comenzar formalmente con el análisis miremos algunas propiedades.

library(dplyr)
df <- penguins %>%
  filter(!is.na(sex),
         !is.na(body_mass_g))

colnames(df) #nos da los nombres
## [1] "species"           "island"            "bill_length_mm"   
## [4] "bill_depth_mm"     "flipper_length_mm" "body_mass_g"      
## [7] "sex"               "year"
summary(df) #resumen general
##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :146   Biscoe   :163   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :123   1st Qu.:39.50   1st Qu.:15.60  
##  Gentoo   :119   Torgersen: 47   Median :44.50   Median :17.30  
##                                  Mean   :43.99   Mean   :17.16  
##                                  3rd Qu.:48.60   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172       Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190       1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197       Median :4050                Median :2008  
##  Mean   :201       Mean   :4207                Mean   :2008  
##  3rd Qu.:213       3rd Qu.:4775                3rd Qu.:2009  
##  Max.   :231       Max.   :6300                Max.   :2009
library(ggplot2)
ggplot(df, aes(x = species, y = body_mass_g, fill = sex))+
  geom_boxplot()+
  theme_light()

2 Test de Hipótesis

2.1 Pruebas Paramétricas

2.1.1 Caso 1: Diferencia de medias con varianzas conocidas (z test).

Vamos a utilizar la columna Sex como separador de grupos, así no tenemos que realizar demasiadas modificaciones.

El punto es que acá no conocemos la varianza poblacional: bue, no importa jaja. Las vamos a “asumir” de la varianza de esta base de datos.

GrupoFem <- df %>% filter(sex == "female") %>% pull(body_mass_g)
GrupoMale <- df %>% filter(sex == "male") %>% pull(body_mass_g)

x1 <- mean(GrupoMale) #Media hombres
x2 <- mean(GrupoFem) #Media mujeres
n1 <- length(GrupoMale) #N de hombres
n2 <- length(GrupoFem) #N de mujeres
sigmaMale <- var(GrupoMale) #sigma hombres
sigmaFem <- var(GrupoFem) #sigma mujeres

#Estadístico
z <- (x1 - x2) / sqrt(sigmaMale / n1 + sigmaFem / n2)

#Cálculo del p valor
pvalue <- 2 * (1 - pnorm(abs(z)))

#Resultados
library(knitr) #librería para tablas
Tabla1 <- data.frame(
  métrica = c("Media male","Media female","z","p"),
  value = c(x1,x2,z,pvalue)
)
kable(Tabla1, caption = "Test z")
Test z
métrica value
Media male 4545.684524
Media female 3862.272727
z 8.554537
p 0.000000

2.1.2 Caso 2: Diferencia de medias con varianzas conocidas (z test).

A partir de acá vamos a obviar tener que definir cada parámetro por nuestra porque R ya nos da funciones predeterminadas.

t.test(body_mass_g ~ sex, data = df)
## 
##  Welch Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -8.5545, df = 323.9, p-value = 4.794e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -840.5783 -526.2453
## sample estimates:
## mean in group female   mean in group male 
##             3862.273             4545.685

2.1.3 Caso 3: Diferencia de medias, muestras pequeñas, varianzas iguales.

Agregamos la varianza combinada (pool), el estadístico t y los grados de libertad.

t.test(body_mass_g ~ sex, data = df, var.equal = T)  
## 
##  Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -8.5417, df = 331, p-value = 4.897e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -840.8014 -526.0222
## sample estimates:
## mean in group female   mean in group male 
##             3862.273             4545.685

2.1.4 Caso 4: Muestras relacionadas (antes/después o pareadas).

En este caso la base de datos no cumple nuestros requisitos así que voy a introducir una pequeña simulación en los datos.

set.seed(43) #seteo semilla.
TimePinguis <- penguins %>%
  filter(!is.na(body_mass_g)) %>%
  slice_sample(n = 30) %>%
  mutate(
    antes = body_mass_g,
    despues = body_mass_g + rnorm(30, mean = 25, sd = 10))

t.test(TimePinguis$antes, TimePinguis$despues, paired = T)
## 
##  Paired t-test
## 
## data:  TimePinguis$antes and TimePinguis$despues
## t = -13.235, df = 29, p-value = 8.096e-14
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -25.33336 -18.55158
## sample estimates:
## mean difference 
##       -21.94247

2.1.5 Caso 5: Comparación de Varianzas (homogeneidad).

var.test(GrupoMale,GrupoFem)
## 
##  F test to compare two variances
## 
## data:  GrupoMale and GrupoFem
## F = 1.3979, num df = 167, denom df = 164, p-value = 0.03198
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.029398 1.897454
## sample estimates:
## ratio of variances 
##           1.397882

2.2 Pruebas No Paramétricas

wilcox.test(body_mass_g ~ sex, data = df, exact = T)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  body_mass_g by sex
## W = 6874.5, p-value = 1.813e-15
## alternative hypothesis: true location shift is not equal to 0

3 Chi Cuadrado

Como se esperan datos categóricos vamos a usar Species y Sex.

chicuadrado <- table(df$species, df$sex)
chisq.test(chicuadrado)
## 
##  Pearson's Chi-squared test
## 
## data:  chicuadrado
## X-squared = 0.048607, df = 2, p-value = 0.976

4 Fisher

chicuadrado <- table(df$species, df$sex)
fisher.test(chicuadrado)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  chicuadrado
## p-value = 0.979
## alternative hypothesis: two.sided