---
title: "Code for Session 9"
author: "Dr. Richard M. Crowley"
format:
  html:
    code-tools: true
    code-link: true
    embed-resources: true
---

Note that the directories used to store data are likely different on your computer, and such references will need to be changed before using any such code.

```{r helpers, warning=FALSE, message=F}
library(knitr)
library(kableExtra)

#' Pretty formatted HTML tables
#' 
#' @param df The dataframe to format as a table
#' @param names Names to use for columns if not using the df column names. Vector
#' @param highlight_cols Columns to highlight. Default is c(1) for highlighting the first row. Use c() to remove.
#' @param highlight_rows Rows to highlight. Default is c() to not highlight any rows.
#' @param color_rows Rows to highlight in color. Default is c() to not highlight any rows.
#' @param highlight_end_rows Rows to highlight. Uses a reverse indexing where 1 is the last row, 2 is second to last, etc.
#' @param full Use to span the full width of a slide. FALSE by default.
#' @param fixed_header Whether the first row should be treated as a header row. TRUE by default.
html_df <- function(df, names=NULL, highlight_cols=c(1), highlight_rows=c(), color_rows=c(), highlight_end_rows=c(), full=F, fixed_header=T) {
  if(!length(names)) {
    names=colnames(df)
  }
  kbl(df,"html", col.names = names, align = c("l",rep('c',length(cols)-1))) %>%
    kable_paper(c("striped","hover"), full_width=full, fixed_thead=F, html_font = "\"Source Sans Pro\", Helvetica, sans-serif") %>%
    {if(length(highlight_cols)) column_spec(., highlight_cols,bold=T, background = "#ffffff99") else .} %>%
    {if(length(highlight_rows)) row_spec(., highlight_rows,bold=T, background = "#ffffff99") else .} %>%
    {if(length(color_rows)) row_spec(., color_rows,bold=T, background = "#aaaaff99") else .} %>%
    {if(length(highlight_end_rows)) row_spec(., nrow(df) + 1 - highlight_end_rows,bold=T, background = "#ffffff99") else .}
}

library(tidyverse)
```

```{r}
#| include: false
options(downlit.attached = c("knitr", "kableExtra", "tidyverse"))
```

```{r, message=F, warnings=F, message=F, results=F}
df <- readRDS('../../Data/Session_6_models.rds')
df2 <- df %>% mutate(AAER = ifelse(AAER==0,0,1))
library(xgboost)

# Prep data
train_x <- model.matrix(AAER ~ ., data=df2[df2$Test==0,-1])[,-1]
train_y <- model.frame(AAER ~ ., data=df2[df2$Test==0,])[,"AAER"]
test_x <- model.matrix(AAER ~ ., data=df2[df2$Test==1,-1])[,-1]
test_y <- model.frame(AAER ~ ., data=df2[df2$Test==1,])[,"AAER"]

set.seed(468435)  #for reproducibility
xgbCV <- xgb.cv(max_depth=5, eta=0.10, gamma=5, min_child_weight = 4,
                subsample = 0.5, objective = "binary:logistic", data=train_x,
                label=train_y, nrounds=100, eval_metric="auc", nfold=10,
                stratified=TRUE, verbose=0)
fit_ens <- xgboost(params=xgbCV$params, data = train_x, label = train_y,
                   nrounds = which.max(xgbCV$evaluation_log$test_auc_mean),
                   verbose = 0)
```

```{r, warning=F, message=F, fig.height=3.5}
library(yardstick)
all_x <- model.matrix(AAER ~ ., data=df[,-1])[,-1]
df$pred_ens <- predict(fit_ens, all_x, type="response")

df_test <- df %>% filter(Test==1)

aucs <- c(df_test %>% roc_auc(AAER, pred_ens, event_level='second') %>% pull(.estimate),
          df_test %>% roc_auc(AAER, pred_BCE, event_level='second') %>% pull(.estimate),
          df_test %>% roc_auc(AAER, pred_lmin, event_level='second') %>% pull(.estimate),
          df_test %>% roc_auc(AAER, pred_xgb, event_level='second') %>% pull(.estimate))
names(aucs) <- c("Ensemble", "Logit (BCE)", "Lasso (lambda.min)", "XGBoost")


curve_out_ens <- df_test %>% roc_curve(AAER, pred_ens, event_level='second')
curve_out_BCE <- df_test %>% roc_curve(AAER, pred_BCE, event_level='second')
curve_out_lmin <- df_test %>% roc_curve(AAER, pred_lmin, event_level='second')
curve_out_xgb <- df_test %>% roc_curve(AAER, pred_xgb, event_level='second')

ggplot() +
  geom_line(data=curve_out_BCE, aes(y=sensitivity, x=1-specificity, color="BCE")) + 
  geom_line(data=curve_out_lmin, aes(y=sensitivity, x=1-specificity, color="LASSO, lambda.min")) + 
  geom_line(data=curve_out_xgb, aes(y=sensitivity, x=1-specificity, color="XGBoost")) +
  geom_line(data=curve_out_ens, aes(y=sensitivity, x=1-specificity, color="Ensemble")) +
  geom_abline(slope=1)
```

```{r}
aucs  # Out of sample
```

```{r, fig.height=4.5}
xgb.train.data = xgb.DMatrix(train_x, label = train_y, missing = NA)
col_names = attr(xgb.train.data, ".Dimnames")[[2]]
imp = xgb.importance(col_names, fit_ens)
# Variable importance
xgb.plot.importance(imp)
```

```{r, message=F, warning=F}
library(gender)
df <- read_csv('../../Data/Session_9-City_of_Chicago_Salary_2019.04.csv')
names(df) <- make.names(names(df),unique = TRUE)
df <- df %>%
  mutate(Full.Time = ifelse(Full.or.Part.Time == 'F',1 ,0),
         Salaried = ifelse(Salary.or.Hourly == 'Salary', 1, 0),
         Salary = ifelse(is.na(Annual.Salary), Typical.Hours * Hourly.Rate, Annual.Salary),
         Obs = n()) %>%
  group_by(Job.Titles) %>%
  mutate(Job.Titles.Freq = n()/Obs) %>%
  ungroup() %>%
  group_by(Department) %>%
  mutate(Department.Freq = n()/Obs) %>%
  ungroup() %>%
  mutate(Job.Titles = ifelse(Job.Titles.Freq > 0.03, Job.Titles, 'Other'),
         Department = ifelse(Department.Freq > 0.05, Department, 'Other'),
         Job.Titles = factor(Job.Titles),
         Department = factor(Department),
         first_name = str_extract(df$Name, '(?<=,[:space:]{2})[:alpha:]+'))

genders <- gender(df$first_name) %>% distinct(name, .keep_all = T)
genders <- genders %>%
  mutate(Female = ifelse(proportion_female > proportion_male,1,0))
df <- left_join(df, genders[, c('name', 'Female')], by = c("first_name" = "name"))

train_x <- model.matrix(Salary ~ Job.Titles + Department + Full.Time + Salaried + Female, data=df)[,-1]
train_y <- model.frame(Salary ~ Job.Titles + Department + Full.Time + Salaried + Female, data=df)[, "Salary"]

set.seed(654687)  #for reproducibility
xgbCV <- xgb.cv(max_depth=6, eta=0.30, gamma=0, min_child_weight = 1,
                subsample = 1, booster = 'gbtree', data=train_x,
                label=train_y, nrounds=100, eval_metric="rmse", nfold=10,
                verbose=0)

fit_xgb <- xgboost(params=xgbCV$params, data = train_x, label = train_y,
                   nrounds = which.min(xgbCV$evaluation_log$test_rmse_mean))
```

```{r, warning=F}
library('SHAPforxgboost')
#https://liuyanguu.github.io/post/2019/07/18/visualization-of-shap-for-xgboost/



vals <- shap.values(xgb_model = fit_xgb, X_train = train_x)
vals$mean_shap_score
```

```{r}
#| eval: true
#| echo: true

vals$shap_score[train_x[, 'Female'] == 1] %>%
  colMeans()

```

```{r}
#| eval: true
#| echo: true

vals$shap_score[train_x[, 'Female'] == 0] %>%
  colMeans()
```

```{r, fig.height = 4.5, fig.width = 8, message=F}
shap_long <- shap.prep(xgb_model = fit_xgb, X_train = train_x)
shap.plot.summary(shap_long)
```

```{r}
#| eval: true
#| echo: true
shap.plot.dependence(data_long = shap_long, x = 'Female', y = 'Female',
                           color_feature = 'Job.TitlesPOLICE OFFICER', jitter_width = 0.05)

```

```{r, fig.height=4.5, fig.width=6, message=F}
plot_data <- shap.prep.stack.data(shap_contrib = vals$shap_score,
                                  top_n = 5, n_groups = 6)
shap.plot.force_plot_bygroup(plot_data)
```

```{r}
g1 <- readRDS('../../Data/Session_9_wmt_1.rds')
g2 <- readRDS('../../Data/Session_9_wmt_2.rds')
g3 <- readRDS('../../Data/Session_9_wmt_3.rds')
g4 <- readRDS('../../Data/Session_9_wmt_4.rds')
g5 <- readRDS('../../Data/Session_9_wmt_5.rds')
g6 <- readRDS('../../Data/Session_9_wmt_6.rds')
```

```{r, fig.height=5, fig.width=6}
g5
```

```{r, fig.height=5, fig.width=6}
g6
```

```{r, fig.height=5, fig.width=6}
g3
```

```{r, fig.height=5, fig.width=6}
g4
```

```{r, fig.height=5, fig.width=6}
g2
```

```{r, fig.height=6, fig.width=10}
g1
```

