12.analyze-affiliation.Rmd

---
title: 'Affiliation analysis'
---

```{r setup, include=FALSE}
library(tidyverse)
library(lubridate)
library(rnaturalearth)
library(epitools)
source('utils/r-utils.R')
library(DT)
theme_set(
  theme_bw() + 
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major.y = element_blank(),
      legend.title = element_blank()
    ))
```

## Country level analysis

Observed vs. expected

```{r}
load('Rdata/raws.Rdata')

iscb_aff_country <- 
  keynotes %>% 
  separate_rows(afflcountries, sep = '\\|') %>% 
  filter(!is.na(afflcountries)) %>% 
  add_count(year, conference, full_name, name = 'num_affls') %>%
  mutate(probabilities = 1 / num_affls,
         publication_date = ymd(year, truncated = 2),
         year = ymd(year, truncated = 2)) %>%
  left_join(nat_to_reg, by = c('afflcountries' = 'country_name'))

pubmed_aff_country <- corr_authors %>%
  filter(!is.na(countries)) %>% 
  add_count(pmid, name = 'num_corr_authors') %>% # number of corresponding authors per pmid
  select(pmid, journal, publication_date, year, countries, fore_name_simple, last_name_simple, num_corr_authors) %>%
  separate_rows(countries) %>%
  add_count(pmid, fore_name_simple, last_name_simple, name = 'num_affls') %>% 
  mutate(probabilities = 1 / num_corr_authors / num_affls)
  
country_rep <- iscb_aff_country %>% 
  group_by(countries) %>% 
  summarise(Observed = sum(probabilities)) %>% 
  arrange(desc(Observed)) 

num_papers <- corr_authors %>% 
  filter(!is.na(countries)) %>% 
  pull(pmid) %>% 
  unique() %>% 
  length()

num_papers

num_honorees <- sum(iscb_aff_country$probabilities)

obs_vs_exp_all <- pubmed_aff_country %>%
  group_by(countries) %>%
  summarise(num_authors = sum(probabilities)) %>%
  ungroup() %>% 
  mutate(freq_affi = num_authors / sum(num_authors)) %>%
  arrange(desc(num_authors)) %>%
  mutate(Expected = freq_affi * num_honorees) %>%
  left_join(country_rep, by = 'countries') %>%
  left_join(nat_to_reg, by = 'countries') %>%
  select(country_name, everything()) %>%
  select(-c(region, countries)) %>%
  mutate(
    Observed = replace_na(Observed, 0),
    over_rep = Observed - Expected,
    other_honorees = num_honorees - Observed,
    other_authors = num_papers - num_authors
  )
```


### Fisher's exact test and 
#### Table of representations
Null hypothesis: For each country, the proportion of honorees affiliated with an institution/company from that country is similar to the proportion of authors affiliated with an institution/company from that country.

```{r warning=FALSE, fig.width = 7*1.5, fig.height = 2.5*1.5}
my_fish <- function(df) {
  res <- df %>%
    unlist() %>%
    matrix(ncol = 2, byrow = TRUE) %>%
    my_riskratio(correction = TRUE)
  
  res_fish <- df %>%
    unlist() %>%
    matrix(ncol = 2) %>%
    fisher.test()
  
  or <- res_fish$estimate
  l_or <- res_fish$conf.int[1]
  u_or <- res_fish$conf.int[2]
  p0 <- df[2]/(df[2] + df[1])
  fish_rr <- or/(1 - p0 + p0*or)
  fish_rr_lower <- l_or/(1 - p0 + p0*l_or)
  fish_rr_upper <- u_or/(1 - p0 + p0*u_or)
  
  res$measure[2, 1:3] %>% 
    c(p_value = res$p.value[2,2]) %>% 
    as.matrix() %>% t() %>% data.frame()
}

nested_obs_exp <- obs_vs_exp_all %>%
  filter(!is.na(country_name)) %>% 
  select(country_name, other_authors, num_authors, other_honorees, Observed) %>% 
  group_by(country_name) %>% 
  nest() 

fish_obs_exp <- nested_obs_exp %>% 
  mutate(fish = map(data, my_fish)) %>% 
  dplyr::select(-data) %>% 
  unnest()

```

```{r eval=FALSE, include=FALSE}
# manual check
x = nested_obs_exp %>% 
  filter(country_name == 'Finland') %>% 
  ungroup() %>% 
  select(data) %>% 
  unlist() %>%
  matrix(ncol = 2, byrow = TRUE) 

a0 <- x[1, 2]
b0 <- x[1, 1]
a1 <- x[2, 2]
b1 <- x[2, 1]
n1 <- a1 + b1
n0 <- a0 + b0
m0 <- b0 + b1
m1 <- a0 + a1

log2(n0/n1*(a1+1)/(a0)*qf(1-0.025, 2*(a1+1), 2*a0))
log2(n0/n1*(a1)/(a0 + 1)/qf(1-0.025, 2*(a0 + 1), 2*a1))
```


### Country enrichment table {#enrichment_tab}
```{r}
fish_obs_exp %>% 
  mutate(ci = paste0('(', round(log2(lower), 1), ', ', 
                     round(log2(upper), 1), ')'),
         lestimate = log2(estimate)) %>% 
  left_join(obs_vs_exp_all, by = 'country_name') %>% 
  select(country_name, freq_affi, Observed, Expected, over_rep,
         # log2fc, 
         estimate, lestimate, ci) %>% 
  rename('Country' = 'country_name',
         'Author proportion' = 'freq_affi',
         'Observed - Expected' = 'over_rep',
         'Enrichment' = 'estimate',
         'Log2(enrichment)' = 'lestimate',
         '95% Confidence interval' = 'ci') %>% 
  datatable(rownames = FALSE) %>% 
  formatPercentage('Author proportion', 2) %>% 
  formatRound(c('Observed', 'Expected', 'Observed - Expected', 
                'Enrichment', 'Log2(enrichment)'), 1)
```


```{r include = F, eval = F, fig.width = 7*1.5, fig.height = 2.5*1.5}
# world map of enrichment
enrich_df <- world %>% 
  left_join(
    fish_obs_exp  %>% 
      mutate(log_lower = case_when(
        lower > 1 ~ log2(lower),
        upper < 1 ~ log2(upper),
        TRUE ~ 0)
      ), by = c('name' = 'country_name'))

enrichment_map <- enrich_df %>% 
  ggplot() +
  geom_sf(aes(fill = log_lower)) +
  scale_fill_gradientn(
    colours = c('#3CBC75FF','white','#440154FF'),
    na.value = NA,
    values = scales::rescale(
      c(min(enrich_df$log_lower, na.rm = T),
        0,
        max(enrich_df$log_lower, na.rm = T)))
  ) +
  coord_sf(crs = '+proj=eqearth +wktext')
enrichment_map
# ggsave('figs/enrichment-map.png', enrichment_map, width = 7, height = 2.5)
# equal earth map projection:
# http://equal-earth.com/equal-earth-projection.html
```

#### Compute enrichment from proportion comparisons
Adapted from `epitools::riskratio()`.

Presentation enrichment/depletion of 20 countries that have the most publications.

```{r fig.width = 7, fig.height = 3.5}
filtered_obs_exp <- obs_vs_exp_all %>% 
  left_join(fish_obs_exp, by = 'country_name') %>% 
  top_n(25, num_authors) %>% 
  mutate(
    distance_to_null = case_when(
      lower > 1 ~ lower - 1,
      TRUE ~ upper - 2
    ),
    presentation = case_when(
      lower > 1 ~ '#377EB8', 
      upper < 1 ~ '#E41A1C',
      TRUE ~ 'grey20'
    ),
    country_name =  as.factor(country_name) %>% fct_reorder(num_authors)) %>% 
  arrange(desc(num_authors))

plot_obs_exp <- filtered_obs_exp %>%
  mutate(lestimate = log2(estimate),
         llower = log2(lower), 
         lupper = log2(upper)) %>% 
  select(country_name, Expected, Observed, lestimate, llower, lupper, presentation, over_rep) %>% 
  pivot_longer(- c(country_name, presentation, over_rep), names_to = 'type') %>% 
  mutate(subtype = ifelse(type == 'Expected' | type == 'Observed', 'Sqrt(number of honorees)', 'Log2 enrichment, 95% CI')) 

```

```{r}
save(plot_obs_exp, file = 'Rdata/affiliations.Rdata')
```

### Log enrichment figure {#enrichment_plot}
```{r}
plot_obs_exp_right <- plot_obs_exp %>% filter(subtype == 'Sqrt(number of honorees)')
plot_obs_exp_left <- plot_obs_exp %>% filter(subtype != 'Sqrt(number of honorees)')

enrichment_plot_left <- plot_obs_exp_left %>%
  ggplot(aes(x = country_name)) +
  coord_flip() +
  labs(x = NULL, y = bquote(Log[2] ~ 'enrichment, 95% CI')) +
  theme(
    legend.position = c(0.9, 0.3),
    axis.title = element_text(size = 9),
    plot.margin = margin(5.5, 2, 5.5, 5.5, unit = 'pt')
  ) +
  scale_color_brewer(type = 'qual', palette = 'Set1') +
  geom_pointrange(
    data = plot_obs_exp_left %>%
      pivot_wider(names_from = type),
    aes(y = lestimate,
        ymin = llower,
        ymax = lupper,
        group = country_name),
    color = filtered_obs_exp$presentation,
    stroke = 0.3, fatten = 2
  ) +
  scale_x_discrete(position = 'top', labels = NULL) +
  scale_y_reverse() +
  geom_hline(data = plot_obs_exp_left, aes(yintercept = 0), linetype = 2)

overrep_countries <- plot_obs_exp_right %>% 
  filter(over_rep > 0) %>% 
  pull(country_name)

enrichment_plot_right <- plot_obs_exp_right %>%
  ggplot(aes(x = country_name, y = value)) +
  geom_line(aes(group = country_name),
            color = rev(plot_obs_exp_right$presentation)) +
  geom_point(aes(shape = type), color = 'grey20') +
  labs(x = NULL, y = 'Number of honorees') +
  theme(
    axis.title = element_text(size = 9),
    legend.position = c(0.75 , 0.3),
    axis.text.y = element_text(
      color = rev(filtered_obs_exp$presentation),
      hjust = 0.5),
    plot.margin = margin(5.5, 5.5, 8.5, 2, unit = 'pt')
  ) +
  scale_y_sqrt(breaks = c(0, 1, 4, 16, 50, 120, 225)) +
  # scale_color_brewer(type = 'qual', palette = 'Set1') +
  coord_flip() +
  geom_text(
    data = . %>% 
      filter(type == 'Expected', !(country_name %in% overrep_countries)),
    nudge_y = 1.2, aes(label = round(over_rep, 1)), size = 2.5) +
  geom_text(
    data = . %>% 
      filter(type == 'Expected', country_name %in% overrep_countries), 
    nudge_y = -1.2, aes(label = round(over_rep, 1)), size = 2.5)

enrichment_plot <- cowplot::plot_grid(enrichment_plot_left, enrichment_plot_right,
                                rel_widths = c(1, 1.3))
enrichment_plot
ggsave('figs/enrichment-plot.png', enrichment_plot, width = 5.5, height = 3.5, dpi = 600)

```


```{r}
sessionInfo()
```