KNN vs. XGBoost Rivalry: Women Employment in Management

Finding a high-profile job position has been very hard for women especially those living in countries with few opportunities for related acquires. This problem can be stemmed from many reasons like contextual factors and accessibility dimensions. This article will examine females’ high-profile job positions rate in the working environment worldwide from the aforementioned perspective.

I prefer the female share of employment in senior and middle management (%) indicator for our target variable modeled on it. The explanatory variables I chose and their explanations link are in the below code block. Besides those, the region and income variables have been included.

library(tidyverse)library(WDI)#Female share of employment in senior and middle management (%)df_profile <-   WDI(indicator = "SL.EMP.SMGT.FE.ZS", extra = TRUE) %>%   as_tibble() %>%   rename(senior_middle_management = SL.EMP.SMGT.FE.ZS) %>%   crosstable::remove_labels()#Employment in industry, female (% of female employment)#(https://data.worldbank.org/indicator/SL.IND.EMPL.FE.ZS)df_industry <-   WDI(indicator = "SL.IND.EMPL.FE.ZS", extra = TRUE) %>%   as_tibble() %>%   rename(industry_employment = SL.IND.EMPL.FE.ZS) %>%   crosstable::remove_labels()#Women Business and the Law Index Score (scale 1-100)#(https://data.worldbank.org/indicator/SG.LAW.INDX)df_wbl <-   WDI(indicator = "SG.LAW.INDX", extra = TRUE) %>%   as_tibble() %>%   rename(wbl = SG.LAW.INDX) %>%   crosstable::remove_labels()#Proportion of seats held by women in national parliaments (%)#(https://data.worldbank.org/indicator/SG.GEN.PARL.ZS)df_seats <-   WDI(indicator = "SG.GEN.PARL.ZS", extra = TRUE) %>%   as_tibble() %>%   rename(seats = SG.GEN.PARL.ZS) %>%   crosstable::remove_labels()  #Merging all the datasetsdf_merged <-   reduce(list(df_profile,              df_industry,              df_wbl,              df_seats),         inner_join,         by = c("country","year","region","income")) %>%   select(region,          income,         wbl,         seats,         senior_middle_management,         industry_employment) %>%  mutate(across(wbl:industry_employment,  ~ round(., digits = 2))) %>%   drop_na()

Before the modeling, we will examine the distributions of the target variable by region and income level to understand its characteristics.

#Exploratory Data Analysis (EDA)library(kableExtra)#Female management by regionfemale_management_region <-   split(df_merged$senior_middle_management,         df_merged$region)inline_plot_region <-   data.frame(var = c(df_merged$region %>% unique()),              Boxplot = "",              Histogram = "")#Female management by incomefemale_management_income <-   split(df_merged$senior_middle_management,         df_merged$income)inline_plot_income <-   data.frame(var = c(df_merged$income %>% unique()),              Boxplot = "",              Histogram = "")#Kable table of the distributions of female share of employment #in senior and middle management rate by region and incomeplot_df <-   inline_plot_region %>%   rbind(inline_plot_income)plot_df %>%   kbl(    col.names = c("", "Boxplot", "Histogram"),    caption = "<center><b>The Distributions of Female Management (%)<br>by Region and Income</b></center>") %>%  kable_paper("striped", full_width = F) %>%  pack_rows("Region", 1, 7) %>%  pack_rows("Income", 8, 11) %>%   column_spec(2, image = c(spec_boxplot(female_management_region),                           spec_boxplot(female_management_income))) %>%   column_spec(3, image = c(spec_hist(female_management_region),                            spec_hist(female_management_income))) %>%   kable_classic(full_width = F,                 html_font = "Bricolage Grotesque") %>%   kable_styling(position = "center")

When we look at the High income, we can see that it is right-skewed, indicating the majority of the observations are below the average. The same is applied to Latin America & Caribbean.

Now, we can pass the modeling phase; while doing that we will compare the algorithms below.

• Bayesian additive regression trees (BART)
• Boosted trees via xgboost
• K-nearest neighbors (KNN)
• Generalized additive models (GAMs)
• Linear regression via keras/tensorflow

#Modelinglibrary(tidymodels)##Split into a train and test setset.seed(12345)df_split <- initial_split(df_merged,                           prop = 0.9,                          strata = region)df_train <- training(df_split)df_test <- testing(df_split)#Resampling for comparing many modelsset.seed(12345)df_folds <-   vfold_cv(df_train, strata = region)#Bayesian additive regression trees (BART)spec_bart <-   parsnip::bart(trees = 20) %>%  set_mode("regression") %>%   set_engine("dbarts")#Boosted trees via xgboostspec_boost <-   boost_tree(trees = 20) %>%  set_mode("regression") %>%   set_engine("xgboost")#K-nearest neighbors via kknnspec_knn <-   nearest_neighbor(neighbors = 5,                    weight_func = "triangular") %>%  set_mode("regression") %>%  set_engine("kknn") #Generalized additive models via mgcvspec_gen_add <-   gen_additive_mod(select_features = FALSE,                    adjust_deg_free = 10) %>%   set_mode("regression") %>%   set_engine("mgcv")#Linear regression via keras/tensorflowspec_linreg <-   linear_reg(penalty = 0.1) %>%   set_engine("keras")#Workflow sets#Workflow set With only a basic formulano_pre_proc <-   workflow_set(    preproc = list("formula" =                      senior_middle_management ~                      region + income + wbl + seats + industry_employment),    models = list(BART = spec_bart)  )#Preprocessing with featuresrec_features<-   recipe(senior_middle_management ~ ., data = df_train) %>%   step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%   step_zv(all_predictors())#Workflow sets With featureswith_features <-   workflow_set(    preproc = list(dummy = rec_features),    models = list(      Linear = spec_linreg,      XGBoost = spec_boost,      KNN = spec_knn    )  )#Workflow set for GAMwflwset_gam <-   workflow() %>%   add_model(spec_gen_add,             formula = senior_middle_management ~               region + income + wbl + seats + industry_employment) %>%   add_formula(senior_middle_management ~                 region +                 income +                wbl +                 seats +                 industry_employment) %>%   as_workflow_set(GAM = .)#Combining all workflow setsall_wflws <-   bind_rows(no_pre_proc,            with_features,            wflwset_gam) %>%   mutate(wflow_id = gsub("(formula_)|(dummy_)", "", wflow_id))#Evaluating the modelsresamples_ctrl <-  control_grid(    save_pred = TRUE,    parallel_over = "everything",    save_workflow = TRUE  )mods_results <-  all_wflws %>%  workflow_map(    seed = 98765,    resamples = df_folds,    control = resamples_ctrl  )#Accuracy ranking of the modelsdf_rmse <-   mods_results %>%   rank_results() %>%  select(wflow_id, .metric, mean) %>%   filter(.metric == "rmse") %>%   select(Models = wflow_id, RMSE = mean)  df_rsq <-   mods_results %>%   rank_results() %>%  select(wflow_id, .metric, mean) %>%   filter(.metric == "rsq") %>%   select(Models = wflow_id, RSQ = mean)#Accuracy tablelibrary(gt)df_acc <-   df_rmse %>%   full_join(df_rsq, by = "Models")df_acc %>%   #rounding decimal digits down to 2 for all numeric variables  mutate(across(where(is.numeric), ~ round(., 2))) %>%   gt() %>%   data_color(    method = "numeric",    palette = c("red", "green")  ) %>%   cols_align(align = "center", columns = -Models) %>%   #separating alignment of column names from cells-alignment  tab_style(    style = cell_text(align = "left"),    locations = cells_body(      columns = Models    )) %>%   #separating cell body from each other  tab_style(    style = cell_borders(sides = "all",                          color = "white",                         weight = px(12),                          style = "solid"),    locations = cells_body(columns = everything())) %>%   tab_header(title = "Accuracy")%>%   opt_table_font(font = "Bricolage Grotesque")

To the above table, the K-nearest neighbors (KNN) algorithm has the best accuracy metrics. So we will build the explanatory model analysis on that model and the parameters.

#Variable importancelibrary(DALEXtra)#Fitted workflow for KNNset.seed(98765)knn_wflow_fitted <-   workflow() %>%   add_recipe(rec_features) %>%   add_model(spec_knn) %>%   fit(df_train)knn_wflow_fitted %>%   extract_fit_parsnip()#Processed data frame for variable importance calculationimp_data <-   rec_features %>%   prep() %>%   bake(new_data = NULL) #Explainer objectexplainer_knn <-   explain_tidymodels(    knn_wflow_fitted %>% extract_fit_parsnip(),     data = imp_data %>% select(-senior_middle_management),     y = imp_data$senior_middle_management,    label = "",    verbose = FALSE  )#Calculating permutation-based variable importance set.seed(1983)vip_knn <- model_parts(explainer_knn,                        loss_function = loss_root_mean_square,                       type = "difference",                       B = 100,#the number of permutations                       label = "")plot(vip_knn)#Variable importance plotvip_knn %>%   mutate(variable = str_remove_all(variable, "region_|income_")) %>%  #removes (...) and replacing with 'and' instead  mutate(variable = str_replace_all(variable, "\\.{3}"," and ")) %>%   #removes (.) and replacing with blank space  mutate(variable = str_replace_all(variable, "\\.", " ")) %>%  mutate(variable = case_when(    variable == "industry_employment" ~ "Industry Employment",    variable == "seats" ~ "Parliaments Seats",    variable == "wbl" ~ "WBL",    TRUE ~ variable  )) %>%  plot() +  labs(color = "",       x = "",       y = "",       subtitle = "Higher indicates more important",       title = "Factors Affecting Female Employment in Senior and Middle Management") +  theme_minimal(base_family = "Bricolage Grotesque",                base_size = 16) +  theme(legend.position = "none",        plot.title = element_text(hjust = 0.5,                                   size = 14,                                  face = "bold"),        plot.subtitle = element_text(hjust = 0.5, size = 12),        panel.grid.minor.x = element_blank(),        panel.grid.major.y = element_blank(),        plot.background = element_rect(fill = "#FFEBFE"))

As seen in the plot above, the rate of female participation in the industry is the most determinant factor followed by the rate of female parliamentarians and women, business, and law (WBL) score.

But, how do those factors affect our target variable? To answer this, we will calculate the partial dependence profiles.

#A function for plotting partial dependence profiles (PDP)library(rlang)plot_pdp <- function(var){    #Partial dependence profiles  set.seed(1983)  pdp_obj <- model_profile(explainer_knn,                            N = NULL, #for all observations                           variables = var)    x_title <-     var %<>%    #removes region_ or income_    gsub("region_|income_", "", .) %>%     #removes (...) and replacing with 'and' instead    gsub("\\.{3}"," and ", .) %>%     #removes (.) or (_) and replacing with blank space    gsub("[._]", " ", .)    pdp_obj$agr_profiles %>%     as_tibble() %>%    ggplot(aes(`_x_`, `_yhat_`)) +    geom_line(color = "#ffb8fb",               linewidth = 1.2) +    labs(x = glue::glue('{x_title %>% str_to_title()}'),         y = "Female Senior and Middle-Level Management (%)") +    theme_minimal(base_family = "Bricolage Grotesque")}#Combining the plotslibrary(patchwork)p_industry_employment <- plot_pdp("industry_employment")p_parliaments_seats <- plot_pdp("seats")p_wbl <- plot_pdp("wbl")p_lower_middle_income <- plot_pdp("income_Lower.middle.income")p_industry_employment +   p_parliaments_seats +  p_wbl +   p_lower_middle_income +   plot_layout(nrow = 2,               axis_titles = "collect")

According to the above graphs, the effects of those factors are mostly positive but the interesting part is that the effect of WBL score turned positive after %80; also, the lower-middle-income level countries having a positive effect is a bit surprising.