Tidytuesday, KMeans y Tidymodels

Análisis de Datos Económicos: EEUU 📈

En este blog post analizaremos los datos correspondientes a la semana 9 de 2021 de TidyTuesday donde se ven los salarios medianos semanales y el número de personas empleadas por sexo, raza y edad a lo largo del tiempo (en concreto se analiza el período 2010-2020). Tenemos dos opciones para descargar los datos: la primera es utilizar el paquete tidytuesdayR pasando como argumento a la función principal la semana de los juegos de datos que deseamos descargar, mientras que la segunda opción es descargar los datos directamente desde el enlace en el que están alojados:

library(tidytuesdayR)
library(tidyverse)
library(broom)
library(gt)

#Option 1
tidyweek <- '2021-02-23'
tuesdata <- tidytuesdayR::tt_load(tidyweek)

## 
##  Downloading file 1 of 2: `earn.csv`
##  Downloading file 2 of 2: `employed.csv`

earn <- tuesdata$earn

#Option 2
earn <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-02-23/earn.csv')

earn %>% 
  head() %>%
  gt() %>%
  tab_header(title = md("Datos de Salarios Medianos en EEUU"))

Vamos a utilizar el paquete skimr para una primera aproximación al juego de datos. A través de la función skim() obtendremos información muy útil sobre las variables dividas por su tipo:

skimr::skim(earn)

Variable type: character

Variable type: numeric

Podemos observar que las variables sex, race y ethnic_origin tienen un número limitado de niveles (valores únicos) mientras que la variable edad está dividida en más niveles (12). También se puede observar que la variable número de personas contratadas tiene una alta variabilidad, lo cuál nos da una pista de que muy posiblemente tendremos que tratar esta variable si queremos utilizarla de manera efectiva. Vamos a comprobar de manera rápida cuales son los valores de las variables categóricas:

earn %>%
  select_if(is_character) %>%
  map(unique)

## $sex
## [1] "Both Sexes" "Men"        "Women"     
## 
## $race
## [1] "All Races"                 "White"                    
## [3] "Black or African American" "Asian"                    
## 
## $ethnic_origin
## [1] "All Origins"        "Hispanic or Latino"
## 
## $age
##  [1] "16 years and over" "16 to 24 years"    "16 to 19 years"   
##  [4] "20 to 24 years"    "25 years and over" "25 to 54 years"   
##  [7] "25 to 34 years"    "35 to 44 years"    "45 to 54 years"   
## [10] "55 years and over" "55 to 64 years"    "65 years and over"

Lo que observamos es que los valores de la edad se solapan, por lo que no parece claro que lógica siguen los niveles en esa variable. En el siguiente apartado de visualización de datos profundizaremos en esta variable para ver cómo se distribuyen estos valores en función de las otras variables.

Visualización de los Datos 🔎📊

Vamos a ver cómo varía la edad en función del sexo y la raza para ver si podemos ver que se estén utilizando diferentes escalas que no se solapen para diferentes niveles de otras variables. Comprobémoslo:

earn %>%
  filter(sex != "Both Sexes") %>%
  ggplot(aes(race, age, color = sex)) +
  geom_count() +
  facet_wrap(~sex) +
  coord_flip() +
  tidyquant::theme_tq() + 
  ylab(NULL) +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "right") 

Desafortunadamente, vemos que los niveles se siguen solapando para los niveles de las variables race y sex, por lo que esta variable no nos será muy útil.

A continuación vamos a analizar como varía el salario mediano semanal en función del sexo y de la raza:

earn %>%
  filter(sex != "Both Sexes") %>%
  ggplot(aes(median_weekly_earn, fill = sex)) +
  geom_histogram(color = "black") +
  scale_fill_manual(values = c("#FFB612", "#97233F")) +
  tidyquant::theme_tq() +
  facet_wrap(~race, nrow = 2, scales = "free") +
  theme(legend.position = "right")

Aquí ya tenemos las primeras conclusiones interesantes: Podemos ver como los hombres cobran más que las mujeres independientemente de la raza (como puede verse al crecer las barras de color rojo de los últimos tramos de las distribuciones respecto a las azules). También es interesante ver como las personas negras o africanas tienen unos rangos menores que las otras razas (su máximo está en torno a 900, mientras que el máximo de los blancos está por encima de los 1250 y el de los asiáticos por encima de los 1500). También es interesante ver como la zona intermedia de la distribución es mayoritariamente de las mujeres en todas las razas mientras que en los primeros valores vemos un mayor equilibrio entro sexos.

earn %>%
  filter(sex != "Both Sexes") %>%
  ggplot(aes(median_weekly_earn, fill = sex)) +
  geom_histogram(color = "black", position = "dodge") +
  tidyquant::theme_tq() +
  scale_fill_manual(values = c("#FFB612", "#97233F")) +
  scale_x_continuous(limits = c(300, 1800)) +
  facet_wrap(~ethnic_origin, nrow = 2, scales = "free") +
  theme(legend.position = "right")

Se observa el mismo patrón que vimos en el gráfico anterior, donde vemos como por lo general los hombres cobran más que las mujeres. También se observa como los hispanos tienen unos sueldos menores a la media (tanto hombres como mujeres) como puede verse en los rangos que alcanzan las distribuciones (el valor máximo para los hispanos no alcanza los 900 dólares semanales de mediana)

earn %>%
  filter(sex != "Both Sexes") %>%
  filter(ethnic_origin == "All Origins" & !is.na(age)) %>%
  mutate(yearq = as.Date(paste(year, case_when(nchar(quarter*3) == 1 ~ paste("0", quarter*3, sep = ""),
                                       TRUE ~ as.character(quarter*3)), "01", sep = "-"))) %>%
  group_by(sex, race, yearq) %>%
  summarise(median_weekly_earn = median(median_weekly_earn),
            n_persons          = sum(n_persons)) %>%
  ungroup() %>%
  mutate(agrupar = paste(sex, race, sep = "-"),
         agrupar = case_when(agrupar == "Men-Black or African American" ~ "Men-Black/African",
                             agrupar == "Women-Black or African American" ~ "Women-Black/African",
                             TRUE ~ agrupar)) %>%
  group_by(agrupar) %>%
  timetk::plot_time_series(
    .date_var    = yearq,
    .value       = median_weekly_earn,
    .color_var   = agrupar, 
    .smooth      = FALSE,
    .facet_ncol  = 3,
    .line_size   = 0.8, 
    .interactive = FALSE,
    .legend_show = FALSE
  ) +
  ggtitle("Evolución Salario Mediano por Raza y Sexo") +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20)))

Claramente podemos ver una tendencia creciente para todas las razas y sexos entre los años 2010 y 2020, con un incremento superior del salario para las personas asiáticas sobre el resto de grupos estudiados (para ambos sexos, tanto hombres como mujeres). También es interesante comentar que este crecimiento tiene un comportamiento prácticamente lineal para las personas de raza blanca (tanto hombres como mujeres) mientras que para las personas tanto asiáticas como africanas/negras se observan mayores fluctuaciones en el crecimiento, no siendo tan lineal.

earn %>%
  filter(sex != "Both Sexes") %>%
  filter(ethnic_origin == "All Origins" & !is.na(age)) %>%
  mutate(yearq = as.Date(paste(year, case_when(nchar(quarter*3) == 1 ~ paste("0", quarter*3, sep = ""),
                                       TRUE ~ as.character(quarter*3)), "01", sep = "-"))) %>%
  group_by(sex, race, yearq) %>%
  summarise(median_weekly_earn = median(median_weekly_earn),
            n_persons          = sum(n_persons)) %>%
  ungroup() %>%
  mutate(agrupar = paste(sex, race, sep = "-"),
         agrupar = case_when(agrupar == "Men-Black or African American" ~ "Men-Black/African",
                             agrupar == "Women-Black or African American" ~ "Women-Black/African",
                             TRUE ~ agrupar)) %>%
  group_by(agrupar) %>%
  timetk::plot_time_series(
    .date_var    = yearq,
    .value       = n_persons,
    .color_var   = agrupar, 
    .smooth      = FALSE,
    .facet_ncol  = 3,
    .line_size   = 0.8, 
    .interactive = FALSE,
    .legend_show = FALSE
  ) +
  ggtitle("Evolución Personas Empleadas por Raza y Sexo") +
  theme(plot.title = element_text(hjust = 0.5, margin = margin(b = 20)))

Podemos ver claramente el efecto de la pandemia COVID en el año 2020 donde aparecen las caídas abruptas para todos los grupos analizados debido al shock inicial y se ve posteriormente el inicio de la recuperación. También resulta curioso ver como las caídas asiáticas han sido menos marcadas que las del resto de grupos (un ejercicio interesante sería ver cuáles son los motivos que se encuentran detrás de las diferencias en estos comportamientos).

Implementación de KMeans 🆒💎

El primer paso será adaptar el dataset para dejarlo en un formato adecuado para poder aplicar el algoritmo KMeans:

earn_tidy <- earn %>%
  mutate(median_weekly = cut(median_weekly_earn, 50)) %>%
  group_by(sex, race,  median_weekly) %>%
  summarise(n_persons = sum(n_persons)) %>%
  select(sex, race, median_weekly, n_persons) %>%
  pivot_wider(names_from = sex:race, values_from = n_persons, id_cols = median_weekly, values_fill = 0) %>%
  janitor::clean_names() %>%
  mutate(across(where(is.numeric), ~as.numeric(scale(.))))

earn_tidy %>%
  head() %>%
  gt() %>%
  tab_header(title = md("Datos Transformados Wider")) %>%
  opt_align_table_header(align = "left")

Vamos a probar ahora con varios números de clusters para ver posteriormente según la regla del codo cual es el que mejor se adapta a nuestro juego de datos. También almacenamos en el tibble la información extraida a partir de las funciones del paquete broom tidy(), augment() y glance():

kclusts <- tibble(k = 1:9) %>%
  mutate(
    kclust = map(k, ~kmeans(select(earn_tidy, -median_weekly), .x)),
    tidied = map(kclust, tidy),
    glanced = map(kclust, glance),
    augmented = map(kclust, augment, select(earn_tidy, -median_weekly))
  )

kclusts

## # A tibble: 9 x 5
##       k kclust   tidied            glanced          augmented         
##   <int> <list>   <list>            <list>           <list>            
## 1     1 <kmeans> <tibble [1 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 2     2 <kmeans> <tibble [2 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 3     3 <kmeans> <tibble [3 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 4     4 <kmeans> <tibble [4 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 5     5 <kmeans> <tibble [5 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 6     6 <kmeans> <tibble [6 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 7     7 <kmeans> <tibble [7 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 8     8 <kmeans> <tibble [8 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>
## 9     9 <kmeans> <tibble [9 x 15]> <tibble [1 x 4]> <tibble [48 x 13]>

clusters <- kclusts %>%
  unnest(cols = c(tidied))

clusters %>%
  select(-c(kclust, glanced, augmented)) %>%
  head() %>%
  gt() %>%
  tab_header(title = md("Información a nivel Cluster")) %>%
  opt_align_table_header(align = "left")

assignments <- kclusts %>% 
  unnest(cols = c(augmented))

assignments %>%
  select(-c(kclust, tidied, glanced)) %>%
  head() %>%
  gt() %>%
  tab_header(title = md("Dataset con Cluster Asignado")) %>%
  opt_align_table_header(align = "left")

clusterings <- kclusts %>%
  unnest(cols = c(glanced))

clusterings %>%
  select(-c(kclust, tidied, augmented)) %>%
  gt() %>%
  tab_header(title = md("Información a nivel k"))

Vamos a explorar alguno de los resultados a continuación:

p1 <- ggplot(assignments, aes(x = women_all_races, y = women_asian)) +
  geom_point(aes(color = .cluster), alpha = 0.8) + 
  facet_wrap(~ k) +
  tidyquant::theme_tq() +
  theme(legend.position = "right")

p1 + geom_point(data = clusters, size = 6, shape = "x")

Podemos crear una función para explorar las combinaciones que nos interesen de manera más cómoda:

explore_kmeans_results <- function(.assignments, .clusters,  .var1, .var2){
  
  p1 <- ggplot(.assignments, aes(x = {{.var1}}, y = {{.var2}})) +
  geom_point(aes(color = .cluster), alpha = 0.8) + 
  facet_wrap(~ k) +
  tidyquant::theme_tq() +
  theme(legend.position = "right")

  p1 + geom_point(data = .clusters, size = 6, shape = "x")
  
}

explore_kmeans_results(assignments, clusters, men_all_races, men_asian)

También podemos ver cual es el número adecuado de clusters para este dataset utilizando la famosa regla del codo. Para ello necesitamos representar la varianza intraclusters por cada valor de k:

ggplot(clusterings, aes(k, tot.withinss)) +
  geom_line() +
  geom_point() +
  tidyquant::theme_tq() +
  ggtitle("Número Óptimo de Clusters")

A pesar de que en ningún momento tenemos la forma característica del codo en mi opinión posiblemente el número más óptimo de clusters en este caso sea cinco.

Fuentes 📚

• Getting started with k-means and #TidyTuesday employment status

• K-means clustering with tidy data principles

Contacto ✉

Alberto Almuiña, Linkedin, Twitter, Github, Blog.