library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
library(umap)
## Warning: package 'umap' was built under R version 4.4.3
library(plotly)
## Warning: package 'plotly' was built under R version 4.4.3
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
employed <- read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2021/2021-02-23/employed.csv")
## Rows: 8184 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): industry, major_occupation, minor_occupation, race_gender
## dbl (3): industry_total, employ_n, year
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#1: Convert to standrized form

employed_grouped <- employed %>%
    filter(!is.na(employ_n)) %>%
    group_by(occupation = paste(industry, minor_occupation), race_gender) %>%
    summarise(n = sum(employ_n)) %>%
    ungroup()
## `summarise()` has grouped output by 'occupation'. You can override using the
## `.groups` argument.
employed_tidy <- employed_grouped %>%
    #remove total
    filter(race_gender != "TOTAL") %>%
    
    #add total column
    left_join(employed_grouped %>%
                  filter(race_gender == "TOTAL") %>%
                  select(occupation, total = n)) %>%
    #Get pct in total
    mutate(pct = n / total) %>%
    
    #Remove Outliers
    filter(total > 1000) %>%
    select(-n)
## Joining with `by = join_by(occupation)`
employed_tidy
## # A tibble: 1,160 × 4
##    occupation                                          race_gender  total    pct
##    <chr>                                               <chr>        <dbl>  <dbl>
##  1 Agriculture and related Construction and extractio… Asian       7.3 e4 0.0274
##  2 Agriculture and related Construction and extractio… Black or A… 7.3 e4 0.0822
##  3 Agriculture and related Construction and extractio… Men         7.3 e4 0.973 
##  4 Agriculture and related Construction and extractio… White       7.3 e4 0.863 
##  5 Agriculture and related Construction and extractio… Women       7.3 e4 0.0274
##  6 Agriculture and related Farming, fishing, and fore… Asian       5.74e6 0.0139
##  7 Agriculture and related Farming, fishing, and fore… Black or A… 5.74e6 0.0342
##  8 Agriculture and related Farming, fishing, and fore… Men         5.74e6 0.789 
##  9 Agriculture and related Farming, fishing, and fore… White       5.74e6 0.911 
## 10 Agriculture and related Farming, fishing, and fore… Women       5.74e6 0.211 
## # ℹ 1,150 more rows
employed_standard <- employed_tidy %>%
    

    #Standardize
    group_by(race_gender) %>%
    mutate(pct = pct %>% scale() %>% as.numeric()) %>%
    ungroup() %>%
    mutate(total = total %>% log() %>% scale() %>% as.numeric()) 

employed_standard
## # A tibble: 1,160 × 4
##    occupation                                          race_gender  total    pct
##    <chr>                                               <chr>        <dbl>  <dbl>
##  1 Agriculture and related Construction and extractio… Asian       -1.30  -0.539
##  2 Agriculture and related Construction and extractio… Black or A… -1.30  -0.405
##  3 Agriculture and related Construction and extractio… Men         -1.30   1.31 
##  4 Agriculture and related Construction and extractio… White       -1.30   0.725
##  5 Agriculture and related Construction and extractio… Women       -1.30  -1.30 
##  6 Agriculture and related Farming, fishing, and fore… Asian        0.819 -0.928
##  7 Agriculture and related Farming, fishing, and fore… Black or A…  0.819 -1.21 
##  8 Agriculture and related Farming, fishing, and fore… Men          0.819  0.510
##  9 Agriculture and related Farming, fishing, and fore… White        0.819  1.38 
## 10 Agriculture and related Farming, fishing, and fore… Women        0.819 -0.503
## # ℹ 1,150 more rows

#2: Spread to object-characteristics

occupation_demo_tbl <- employed_tidy %>%
    pivot_wider(names_from = race_gender, values_from = pct) %>%
    janitor::clean_names()

#3: Perform k-means clustering

occupation_cluster <- kmeans(occupation_demo_tbl %>% select(-occupation), centers = 3, nstart = 20)
summary(occupation_cluster)
##              Length Class  Mode   
## cluster      232    -none- numeric
## centers       18    -none- numeric
## totss          1    -none- numeric
## withinss       3    -none- numeric
## tot.withinss   1    -none- numeric
## betweenss      1    -none- numeric
## size           3    -none- numeric
## iter           1    -none- numeric
## ifault         1    -none- numeric
tidy(occupation_cluster)
## # A tibble: 3 × 9
##     total  asian black_or_african_ame…¹   men white women  size withinss cluster
##     <dbl>  <dbl>                  <dbl> <dbl> <dbl> <dbl> <int>    <dbl> <fct>  
## 1  2.85e6 0.0449                  0.105 0.679 0.812 0.319   219  3.28e15 1      
## 2  3.58e7 0.0650                  0.127 0.566 0.774 0.434    12  1.92e15 2      
## 3  1.16e8 0.0720                  0.117 0.268 0.785 0.732     1  0       3      
## # ℹ abbreviated name: ¹​black_or_african_american
glance(occupation_cluster)
## # A tibble: 1 × 4
##     totss tot.withinss betweenss  iter
##     <dbl>        <dbl>     <dbl> <int>
## 1 3.00e16      5.20e15   2.48e16     2
augment(occupation_cluster, occupation_demo_tbl) %>%
    
    ggplot(aes(total, asian, color = .cluster)) +
    geom_point()

#4: Select optimal number of clusters

kclusts <- tibble(k = 1:9) %>%
    mutate(kclust = map(.x = k, .f = ~kmeans(occupation_demo_tbl %>% select(-occupation), centers = 3, nstart = 20)),
           glaced = map(.x = kclust, .f = glance))

kclusts %>%
    unnest(glaced) %>%
    ggplot(aes(k, tot.withinss)) +
    geom_point() +
    geom_line()

final_cluster <- kmeans(occupation_demo_tbl %>% select(-occupation), centers = 5, nstart = 20)
augment(final_cluster, occupation_demo_tbl) %>%
    
    ggplot(aes(total, asian, color = .cluster)) +
    geom_point()

#5: Reduce Dimension using UMAP

umap_results <- occupation_demo_tbl %>%
    select(-occupation) %>%
    umap()

umap_results_tble <- umap_results$layout %>%
    as.tibble() %>%
    bind_cols(occupation_demo_tbl %>% select(occupation))
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tibble package.
##   Please report the issue at <https://github.com/tidyverse/tibble/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
umap_results_tble
## # A tibble: 232 × 3
##         V1    V2 occupation                                                     
##      <dbl> <dbl> <chr>                                                          
##  1 -8.10   -6.50 Agriculture and related Construction and extraction occupations
##  2 -0.803   6.24 Agriculture and related Farming, fishing, and forestry occupat…
##  3 -1.21   -7.29 Agriculture and related Installation, maintenance, and repair …
##  4  8.43   -1.46 Agriculture and related Manage-ment, business, and financial o…
##  5  0.0444  6.10 Agriculture and related Management, business, and financial op…
##  6  4.10   -5.63 Agriculture and related Office and administrative support occu…
##  7 -1.25   -7.14 Agriculture and related Production occupations                 
##  8  1.46   -6.22 Agriculture and related Professional and related occupations   
##  9 -7.96   -6.20 Agriculture and related Protective service occupations         
## 10 -7.83   -6.06 Agriculture and related Sales and related occupations          
## # ℹ 222 more rows
umap_results_tble %>%
    ggplot(aes(V1, V2)) +
    geom_point()

#6: Visualize clusters by adding k-means results

kmeans_umap_tbl <- final_cluster %>%
    augment(occupation_demo_tbl) %>%
    select(occupation, .cluster) %>%
    
    #Add umap results
    left_join(umap_results_tble) %>%
    
    # Add employement info
    left_join(employed_tidy %>% 
                  select(-total) %>%
                  pivot_wider(names_from = race_gender, values_from = pct) %>%
                  janitor::clean_names())
## Joining with `by = join_by(occupation)`
## Joining with `by = join_by(occupation)`
kmeans_umap_tbl
## # A tibble: 232 × 9
##    occupation  .cluster      V1    V2   asian black_or_african_ame…¹   men white
##    <chr>       <fct>      <dbl> <dbl>   <dbl>                  <dbl> <dbl> <dbl>
##  1 Agricultur… 1        -8.10   -6.50 0.0274                 0.0822  0.973 0.863
##  2 Agricultur… 5        -0.803   6.24 0.0139                 0.0342  0.789 0.911
##  3 Agricultur… 1        -1.21   -7.29 0.0155                 0.0309  0.985 0.918
##  4 Agricultur… 1         8.43   -1.46 0.00992                0.00794 0.739 0.967
##  5 Agricultur… 5         0.0444  6.10 0.00997                0.00882 0.741 0.962
##  6 Agricultur… 1         4.10   -5.63 0.0233                 0.0155  0.159 0.938
##  7 Agricultur… 1        -1.25   -7.14 0.0332                 0.104   0.815 0.820
##  8 Agricultur… 1         1.46   -6.22 0.0339                 0.0373  0.675 0.902
##  9 Agricultur… 1        -7.96   -6.20 0                      0.0682  0.864 0.875
## 10 Agricultur… 1        -7.83   -6.06 0                      0.0213  0.585 0.968
## # ℹ 222 more rows
## # ℹ abbreviated name: ¹​black_or_african_american
## # ℹ 1 more variable: women <dbl>
g <- kmeans_umap_tbl %>%
    
    #Create text label
    mutate(text_label = str_glue("Occupation: {occupation}
                                 Cluster: {.cluster}
                                 Asian: {asian %>% scales::percent(1)}
                                 Women: {women%>% scales::percent(1)}")) %>%
    #Plot
    ggplot(aes(V1, V2, color = .cluster, text = text_label)) +
    geom_point()

g %>% ggplotly(tooltip = "text")