Plotting ROC curves

Name origin prediction method performance

library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0

## ── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

# still need to install caret for the calibration function because tidymodels's 
# probably hasn't published this yet
library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

source('utils/r-utils.R')
theme_set(theme_bw())

roc_df <- read_tsv('https://raw.githubusercontent.com/greenelab/wiki-nationality-estimate/7c22d0a5f661ce5aeb785215095deda40973ff17/models/NamePrism_roc_curves.tsv') %>%
  rename('region' = category) %>%
  # recode_region_letter() %>%
  recode_region() %>% 
  group_by(region) %>%
  mutate(Sensitivity = tpr, Specificity = 1-fpr, dSens = c(abs(diff(1-tpr)), 0)) %>%
  ungroup()

## 
## ── Column specification ───────────────────────────────────────────────────────────────────────────────
## cols(
##   fpr = col_double(),
##   tpr = col_double(),
##   threshold = col_double(),
##   category = col_character()
## )

## Warning: Problem with `mutate()` input `region`.
## ℹ Unknown levels in `f`: OtherCategories
## ℹ Input `region` is `fct_recode(...)`.

## Warning: Unknown levels in `f`: OtherCategories

auc_df <- roc_df %>%
  group_by(region) %>%
  # add_count() %>%
  summarise(auc = sum((1 - fpr) * dSens),
            n = n()) %>%
  arrange(desc(auc)) %>%
  mutate(auc_pct = 100 * auc,
         reg_auc = paste0(region, ', AUC = ', round(auc_pct, 1), '%'))

## `summarise()` ungrouping output (override with `.groups` argument)

# region_levels <- c('Celtic English', 'European', 'East Asian', 'Hispanic', 'South Asian', 'Muslim', 'Israeli', 'African')
region_levels <- paste(c('Celtic/English', 'European', 'East Asian', 'Hispanic', 'South Asian', 'Arabic', 'Hebrew', 'African', 'Nordic', 'Greek'), 'names')
region_levels_let <- toupper(letters[1:8])
region_cols <- c('#b3de69', '#fdb462',  '#bc80bd', '#8dd3c7', '#fccde5', '#ffffb3', '#ccebc5', '#bebada', '#80b1d3', '#fb8072')

fig_3a <- roc_df %>%
  left_join(auc_df, by = 'region') %>%
  ggplot(aes(x = Sensitivity, y = Specificity, color = fct_relevel(reg_auc, as.character(auc_df$reg_auc)))) +
  scale_color_manual(values = region_cols) +
  geom_step(size = 1, alpha = 0.8) +
  coord_fixed() +
  scale_x_reverse(breaks = seq(1, 0, -0.2), labels = scales::percent) +
  scale_y_continuous(breaks = seq(0, 1, 0.2), labels = scales::percent, limits = c(NA, 1.05)) +
  theme(legend.position = c(0.62, 0.42),
        legend.title = element_blank(),
        legend.text.align = 1,
        legend.text = element_text(size = 7),
        legend.margin = margin(-0.2, 0.2, 0.2, 0, unit='cm'))

predictions_df <- read_tsv('https://raw.githubusercontent.com/greenelab/wiki-nationality-estimate/7c22d0a5f661ce5aeb785215095deda40973ff17/data/NamePrism_results_test.tsv') %>%
  mutate(y_true = as.factor(truth)) %>%
  select(-truth)

## 
## ── Column specification ───────────────────────────────────────────────────────────────────────────────
## cols(
##   African = col_double(),
##   CelticEnglish = col_double(),
##   EastAsian = col_double(),
##   European = col_double(),
##   Greek = col_double(),
##   Hispanic = col_double(),
##   Jewish = col_double(),
##   Muslim = col_double(),
##   Nordic = col_double(),
##   SouthAsian = col_double(),
##   truth = col_character()
## )

regs <- predictions_df %>% select(African:SouthAsian) %>% colnames()
cal_dfs <- list()
for (reg in regs) {
  pred_reg <- predictions_df %>%
    mutate(y_true_bin = as.factor((y_true == reg))) %>%
    rename(prob = reg) %>%
    select(y_true_bin, prob)

  cal_dfs[[reg]] <- calibration(y_true_bin ~ prob,
                                data = pred_reg,
                                cuts = 11,
                                class = 'TRUE')$data %>%
    mutate(region = reg)
}

## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(reg)` instead of `reg` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

cal_dfs$EastAsian

##    calibModelVar            bin   Percent      Lower     Upper Count  midpoint    region
## 1           prob     [0,0.0909]  0.973038  0.9061138  1.043559   777  4.545455 EastAsian
## 2           prob (0.0909,0.182] 12.715105 10.7555376 14.887108   133 13.636364 EastAsian
## 3           prob  (0.182,0.273] 20.620843 16.9791523 24.652952    93 22.727273 EastAsian
## 4           prob  (0.273,0.364] 29.924242 24.4643714 35.841394    79 31.818182 EastAsian
## 5           prob  (0.364,0.455] 35.897436 29.7515540 42.405681    84 40.909091 EastAsian
## 6           prob  (0.455,0.545] 38.536585 31.8402892 45.569554    79 50.000000 EastAsian
## 7           prob  (0.545,0.636] 45.637584 37.4635833 53.988516    68 59.090909 EastAsian
## 8           prob  (0.636,0.727] 56.953642 48.6544756 64.974492    86 68.181818 EastAsian
## 9           prob  (0.727,0.818] 61.421320 54.2394760 68.253900   121 77.272727 EastAsian
## 10          prob  (0.818,0.909] 71.764706 66.6571343 76.488532   244 86.363636 EastAsian
## 11          prob      (0.909,1] 97.209555 96.8524649 97.536348  8953 95.454545 EastAsian

fig_3b <- bind_rows(cal_dfs) %>%
  recode_region() %>%
  ggplot(aes(x = midpoint/100, y = Percent/100, color = fct_relevel(region, as.character(auc_df$region)))) +
  geom_abline(slope = 1, linetype = 2, alpha = 0.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 20L), breaks = seq(0, 1, 0.2), limits = c(-0.005, 1.045)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 20L), breaks = seq(0, 1, 0.2), limits = c(0, 1)) +
  coord_fixed() +
  geom_point() +
  geom_line() +
  scale_color_manual(values = region_cols) +
  theme(legend.position = 'None') +
  labs(x = 'Predicted probability', y = 'Fraction of names')

## Warning: Problem with `mutate()` input `region`.
## ℹ Unknown levels in `f`: OtherCategories
## ℹ Input `region` is `fct_recode(...)`.

## Warning: Unknown levels in `f`: OtherCategories

n_obs <- sum(auc_df$n)
short_regs <- auc_df$region %>% 
  as.character() %>% 
  gsub(' names', '', .)

heat_dat <- predictions_df %>%
  group_by(y_true) %>%
  summarise_if(is.numeric, mean, na.rm = T) %>%
  ungroup() %>%
  pivot_longer(- y_true, names_to = 'region', values_to = 'pred_prob') %>%
  recode_region() %>%
  rename('reg_hat' = region, 'region' = y_true) %>% 
  recode_region() %>%
  rename('y_true' = region, 'region' = reg_hat) %>% 
  left_join(auc_df, by = 'region') %>%
  mutate(scale_pred_prob = log2((pred_prob)/(n/n_obs)),
         region = region %>% gsub(' names', '', .) %>% fct_relevel(short_regs),
         y_true = y_true %>% gsub(' names', '', .) %>% fct_relevel(short_regs))

## Warning: Problem with `mutate()` input `region`.
## ℹ Unknown levels in `f`: OtherCategories
## ℹ Input `region` is `fct_recode(...)`.

## Warning: Unknown levels in `f`: OtherCategories

## Warning: Problem with `mutate()` input `region`.
## ℹ Unknown levels in `f`: OtherCategories
## ℹ Input `region` is `fct_recode(...)`.

## Warning: Unknown levels in `f`: OtherCategories

fig_3c <- ggplot(heat_dat, aes(y_true, region,
                               fill = scale_pred_prob)) +
  geom_tile() +
  scale_fill_gradientn(
    colours = c("#3CBC75FF","white","#440154FF"),
    values = scales::rescale(
      c(min(heat_dat$scale_pred_prob),
        0,
        max(heat_dat$scale_pred_prob)))
  ) +
  coord_fixed() +
  labs(x = 'True region', y = 'Predicted region', fill = bquote(log[2]~'FC')) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        legend.position = 'top',
        legend.key.height = unit(0.2, 'cm'),
        legend.title = element_text(vjust = 1),
        legend.margin = margin(0, 0,0, -1, unit='cm'),
        axis.title.x = element_text(margin = margin(t = 27, r = 0, b = 0, l = 0)),
        axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0)))

fig_3 <- cowplot::plot_grid(fig_3a, fig_3b, fig_3c, labels = 'AUTO', nrow = 1,
                            rel_widths = c(2,2,1.6))
fig_3

# ggsave('figs/fig_3.png', fig_3, height = 4, width = 10)