Ejercicio Contraste Kolmogorov-Smirnov
2023-08-01
Compraciación de distribuciones estadísticas
#=============================================
1.- KOLMOGOROV SMIRNOV
#=============================================
El estadístico de Kolmogorov–Smirnov, también conocido como distancia Kolmogorov–Smirnov (K-S), se define como la distancia vertical máxima entre las funciones de distribución acumulada empíricas de dos muestras, o entre una función de distribución empírica y una función de distribución acumulada teórica de referencia. La ventaja principal de este estadístico es que es sensible a diferencias tanto en la localización como en la forma de la función de distribución acumulada.
EJEMPLO
El set de datos Snmesp del paquete plm contiene una muestra de los salarios (en escala logarítmica) pagados en España durante los años 1983 a 1990 (783 observaciones por año). Se pretende determinar si la distribución de los salarios cambió entre 1989 y 1993. Para ello se calcula la distancia de Kolmogorov–Smirnov.
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'readr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'forcats' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorsdata(Snmesp, package = "plm")
Snmesp <- Snmesp %>%
filter(year %in% c(1989, 1990)) %>%
mutate(year = as.factor(year)) %>%
select(year, salario = w)
Snmesp %>%
ggplot(aes(x = salario, fill = year)) +
geom_density(alpha = 0.3) +
scale_fill_manual(values = c("lightblue", "red")) +
labs(title = "Salarios en España entre 1989 y 1990",
subtitle = "Escala logarítmica, 738 observaciones por año",
fill = "Año") +
theme_bw() +
theme(legend.position = "bottom")
Cálculo de la función de distribución acumulada empírica
La función ecdf() permite ajustar la función de distribución acumulada empírica a partir de una muestra. El resultado de esta función es un objeto ecdf que se comporta de forma similar a un modelo predictivo, recibe un vector de observaciones y devuelve su probabilidad acumulada.
# Se ajustan las funciones ecdf con cada muestra.
ecdf_1989 <- ecdf(Snmesp %>% filter(year == 1989) %>% pull(salario))
ecdf_1990 <- ecdf(Snmesp %>% filter(year == 1990) %>% pull(salario))
# Se calcula la probabilidad acumulada de cada valor de salario observado con cada
# una de las funciones ecdf.
grid_salario <- unique(Snmesp %>% pull(salario))
prob_acumulada_ecdf_1989 <- ecdf_1989(v = grid_salario)
prob_acumulada_ecdf_1990 <- ecdf_1990(v = grid_salario)# Se unen los valores calculados en un dataframe.
df_ecdf <- data.frame(
salario = grid_salario,
ecdf_1989 = prob_acumulada_ecdf_1989,
ecdf_1990 = prob_acumulada_ecdf_1990
) %>%
pivot_longer(
cols = c(ecdf_1989, ecdf_1990),
names_to = "year",
values_to = "ecdf"
)
grafico_ecdf <- ggplot(data = df_ecdf,
aes(x = salario, y = ecdf, color = year)) +
geom_line(size = 0.6) +
scale_color_manual(values = c("lightblue", "red")) +
labs(
title = "Función de distribución acumulada empírica salarios",
subtitle = "Escala logarítmica, 738 observaciones por año",
color = "Año",
y = "Probabilidad acumulada"
) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(size = 12))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.grafico_ecdf
Snmesp %>%
ggplot(aes(x=salario, color = year))+
stat_ecdf(geom = "step")+
labs(title = "Fx",
subtitle = "escala logarítmica",
color = "Año")+
theme_bw()+
theme(legend.position = "bottom",
plot.title = element_text(size=12))
#### Diferencias en distribuciones
# Se calcula la diferencia absoluta entre las probabilidades acumuladas de cada
# función.
abs_dif <- abs(prob_acumulada_ecdf_1989 - prob_acumulada_ecdf_1990)
# La distancia Kolmogorov–Smirnov es el máximo de las distancias absolutas.
distancia_ks <- max(abs_dif)
paste("Distancia Kolmogorov–Smirnov:", distancia_ks)## [1] "Distancia Kolmogorov–Smirnov: 0.105691056910569"indice_ks <- which.max(abs_dif)
grafico_ecdf +
geom_segment(aes(
x = grid_salario[indice_ks],
xend = grid_salario[indice_ks],
y = prob_acumulada_ecdf_1989[indice_ks],
yend = prob_acumulada_ecdf_1990[indice_ks]
),
arrow = arrow(ends = "both", length = unit(0.2,"cm")),
color = "blue")
test <- ks.test(
x = Snmesp %>% filter(year == 1989) %>% pull(salario),
y = Snmesp %>% filter(year == 1990) %>% pull(salario)
)## Warning in ks.test.default(x = Snmesp %>% filter(year == 1989) %>%
## pull(salario), : p-value will be approximate in the presence of tiestest$statistic## D
## 0.1056911test$p.value## [1] 0.0005257129El estadístico del contraste asciende a 0.10 con un p valor de 0.0005<0.05 Tenemos evidencia empírica para rechazar la hipótesis nula de que la distribución de salarios es la misma en los dos años considerados
##2.- Ejemplo PRO #==================================================
# Distribuciones teóricas
p1 <- ggplot(data = data.frame(x=c(-15,15)), aes(x = x)) +
stat_function(fun=dnorm, args=list(mean = 0, sd = 3),
geom = "line", color = "blue") +
stat_function(fun=dnorm, args=list(mean = 0, sd = 3),
geom = "area", fill = "blue", alpha = 0.5) +
stat_function(fun=dnorm, args=list(mean = 2, sd = 3),
geom = "line", color = "red") +
stat_function(fun=dnorm, args=list(mean = 2, sd = 3),
geom = "area", fill = "red", alpha = 0.5) +
labs(title = "Distribuciones normales",
subtitle = "(mean_1 = 0, sd_1 = 1) y (mean_2 = 0, sd_2 = 5)") +
theme_bw()
# Simulación de datos extraídos de dos distribuciones
set.seed(1234)
sample_normal_1 <- rnorm(n = 1000000, mean = 0, sd = 3)
sample_normal_2 <- rnorm(n = 1000000, mean = 2, sd = 3)
# Distribuciones de las muestras simuladas
p2 <- ggplot() +
geom_histogram(data = tibble(valor = sample_normal_1),
aes(x = valor, stat(density)),
fill = "blue",
alpha = 0.5) +
geom_histogram(data = tibble(valor = sample_normal_2),
aes(x = valor, stat(density)),
fill = "red",
alpha = 0.5) +
labs(title = "Distribuciones empíricas",
subtitle = paste(
"mean_1 = ", round(mean(sample_normal_1), 3),
"sd_1 =", round(sd(sample_normal_1), 3),
"mean_2 = ", round(mean(sample_normal_2), 3),
"sd_2 =", round(sd(sample_normal_2), 3)
)
) +
theme_bw()
ggpubr::ggarrange(p1, p2, nrow = 1)## Warning: `stat(density)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
set.seed(1234)
test_ks_boot <- Matching::ks.boot(
Tr = sample_normal_1,
Co = sample_normal_2,
alternative = "two.sided",
nboots = 1000
)
test_ks_boot$ks.boot.pvalue## [1] 0El p valor es cero por tanto las distribuciones de sloarios son diferentes en los dos años.
# Distribuciones teóricas
p1 <- ggplot(data = data.frame(x=c(0,1)), aes(x = x)) +
stat_function(fun=dbeta, args=list(shape1 = 2, shape2 = 5),
geom = "line", color = "gray50") +
stat_function(fun=dbeta, args=list(shape1 = 2, shape2 = 5),
geom = "area", fill = "gray50", alpha = 0.5) +
stat_function(fun=dbeta, args=list(shape1 = 4, shape2 = 10),
geom = "line", color = "orangered2") +
stat_function(fun=dbeta, args=list(shape1 = 4, shape2 = 10),
geom = "area", fill = "orangered2", alpha = 0.5) +
labs(title = "Distribuciones beta",
subtitle = "(shape1 = 2, shape2 = 5) y (shape1 = 4, shape2 = 10)") +
theme_bw()
# Simulación de datos extraídos de dos distribuciones
set.seed(123)
sample_beta_1 <- rbeta(n = 1000, shape1 = 2, shape2 = 5)
sample_beta_2 <- rbeta(n = 1000, shape1 = 4, shape2 = 10)
# Distribuciones de las muestras simuladas
p2 <- ggplot() +
geom_histogram(data = tibble(valor = sample_beta_1),
aes(x = valor, stat(density)),
fill = "gray50",
alpha = 0.5) +
geom_histogram(data = tibble(valor = sample_beta_2),
aes(x = valor, stat(density)),
fill = "orangered2",
alpha = 0.5) +
labs(title = "Distribuciones empíricas",
subtitle = paste(
"mean_1 = ", round(mean(sample_beta_1), 3),
"sd_1 =", round(sd(sample_beta_1), 3),
"mean_2 = ", round(mean(sample_beta_2), 3),
"sd_2 =", round(sd(sample_beta_2), 3)
)
) +
theme_bw()
ggpubr::ggarrange(p1, p2, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
set.seed(123)
test_ks_boot <- Matching::ks.boot(
Tr = sample_beta_1,
Co = sample_beta_2,
alternative = "two.sided",
nboots = 1000
)
test_ks_boot$ks.boot.pvalue## [1] 0El p valor es cero por tanto las distribuciones de salarios son diferentes en los dos años.
# Funciones auxiliares que calculan el estadístico para cada permutación (split)
#-------------------------------------------------------------------------------
calcular_dif_media <- function(split) {
grupo_1 <- rsample::analysis(x = split)
grupo_2 <- rsample::assessment(x = split)
dif <- mean(grupo_1[[1]]) - mean(grupo_2[[1]])
dif_abs <- abs(dif)
return(dif_abs)
}
calcular_dif_mediana <- function(split) {
grupo_1 <- rsample::analysis(x = split)
grupo_2 <- rsample::assessment(x = split)
dif <- median(grupo_1[[1]]) - median(grupo_2[[1]])
dif_abs <- abs(dif)
return(dif_abs)
}
calcular_dif_sd <- function(split) {
grupo_1 <- rsample::analysis(x = split)
grupo_2 <- rsample::assessment(x = split)
dif <- sd(grupo_1[[1]]) - sd(grupo_2[[1]])
dif_abs <- abs(dif)
return(dif_abs)
}
calcular_dif_iqr <- function(split) {
grupo_1 <- rsample::analysis(x = split)
grupo_2 <- rsample::assessment(x = split)
dif <- IQR(grupo_1[[1]]) - IQR(grupo_2[[1]])
dif_abs <- abs(dif)
return(dif_abs)
}
# Función para calcular la diferencia observada y su significancia estadística
#-------------------------------------------------------------------------------
calcular_diferencias_permutacion <- function(sample_1,
sample_2,
n_permut = 1000,
seed = NULL) {
# Se calcula la longitud de cada sample para que el posterior reparto sea proporcional.
len_sample_1 <- length(sample_1)
len_sample_2 <- length(sample_2)
len_pull_samples <- len_sample_1 + len_sample_2
# Se almacenan las observaciones en un dataframe
df <- tibble(
sample = rep(c("sample_1", "sample_2"),
times = c(len_sample_1, len_sample_2)),
valor = c(sample_1, sample_2)
)
# Cálculo del valor observado de cada estadístico y su diferencia.
valor_observado <- df %>%
group_by(sample) %>%
summarise(
media = mean(valor, na.rm = TRUE),
mediana = median(valor, na.rm = TRUE),
sd = sd(valor, na.rm = TRUE),
iqr = IQR(valor, na.rm = TRUE)
) %>%
pivot_longer(
cols = -sample,
names_to = "estadistico",
values_to = "valor"
) %>%
pivot_wider(
names_from = sample,
values_from = valor) %>%
mutate(
diferencia = sample_1 - sample_2,
diferencia_abs = abs(diferencia)
)
# Se generan las permutaciones con la función rsample::mc_cv().
if (!is.null(seed)) {
set.seed(seed = seed)
}
resamples <- rsample::mc_cv(
data = df %>% select(valor),
prop = len_sample_1/len_pull_samples,
times = n_permut
)
# Se paraleliza el cálculo de cada permutación.
future::plan(strategy = "multiprocess")
permut_media <- furrr::future_map_dbl(.x = resamples$splits,
.f = calcular_dif_media)
permut_mediana <- furrr::future_map_dbl(.x = resamples$splits,
.f = calcular_dif_mediana)
permut_sd <- furrr::future_map_dbl(.x = resamples$splits,
.f = calcular_dif_sd)
permut_iqr <- furrr::future_map_dbl(.x = resamples$splits,
.f = calcular_dif_iqr)
# Cálculo de significancia
dif_observada_abs <- valor_observado %>% pull(diferencia_abs)
names(dif_observada_abs) <- valor_observado %>% pull(estadistico)
pvalue_media <- mean(permut_media > dif_observada_abs[["media"]])
pvalue_mediana <- mean(permut_mediana > dif_observada_abs[["mediana"]])
pvalue_sd <- mean(permut_sd > dif_observada_abs[["sd"]])
pvalue_iqr <- mean(permut_iqr > dif_observada_abs[["iqr"]])
resultados <- valor_observado %>%
mutate(p_value = c(pvalue_media, pvalue_mediana,
pvalue_sd,pvalue_iqr)
)
return(resultados)
}