figures.Rmd

---
title: "R Notebook"
---

# Load libraries

```{r}

library(tidyverse)
library(magrittr)
library(lubridate)
library(ggalluvial)
library(ggrepel)
library(corrr) # for correlation matrix
library(patchwork) # for creating nice plot grids
library(colorspace) # better color palettes

```

# Theme options

```{r}

# Spruce up the default figure formatting
theme_set(theme_minimal() +
            theme(text = element_text(size = 10),
                  axis.title.x = element_text(size = 10,
                                              margin = margin(5, 0, 5, 0)),
                  axis.text.x = element_text(size = 8),
                  axis.title.y = element_text(size = 10,
                                              margin = margin(0, 5, 0, 5)),
                  axis.text.y = element_text(size = 8),
                  plot.title = element_text(size = 10),
                  panel.border = element_rect(color = "#E0E0E0", 
                                              size = 0.5, 
                                              fill = NA),
                  plot.margin = margin(5,5,5,5),
                  legend.key.size = unit(0.5, "cm"),
                  legend.text = element_text(size = 8)))

# Make sure axis labels (e.g. month names) are in English
Sys.setlocale("LC_ALL", "English")

# Create custom date scales for time-series plots
custom_week_scale <- function(start_date, end_date, format = "%d-%b") {
  scale_x_date(date_breaks = "1 week",
               date_minor_breaks = "1 day",
               date_labels = format,
               expand = c(0, 0),
               limits = c(ymd(start_date), ymd(end_date)))
}
custom_month_scale <- function(start_date, end_date, format = "%d-%b") {
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = format,
               expand = c(0, 0),
               limits = c(ymd(start_date), ymd(end_date)))
}

```

# Helper variables and functions

```{r}

analysis_start <- ymd("2020-01-01")
analysis_end <- ymd("2020-10-31")

assign_covid_preprint <- function(d) {
  d %>%
    mutate(covid_preprint = case_when(
      covid_preprint == T ~ "COVID-19 preprints",
      covid_preprint == F ~ "non-COVID-19 preprints"
    ))
}

# Function for basic cleaning of author names - convert to lowercase, remove whitespace, remove periods after initials
clean_names <- function(d) {
  d %>% 
    mutate(author_corresponding = str_to_lower(str_trim(str_remove(author_corresponding, "\\."))))
}

```


# Load datasets

```{r}

# Dataset of COVID-19 cases and deaths. Data is aggregated from the repository
# maintained by Johns Hopkins (https://github.com/CSSEGISandData/COVID-19)
covid <- read_csv("https://raw.githubusercontent.com/datasets/covid-19/master/data/worldwide-aggregate.csv")

# Full dataset of preprints from the launch of bioRxiv (Nov 2013) to the end of
# our study period (October 2020).
preprints_all <- read_csv("data/preprints_basic_20131101_20201031.csv")

# Abstracts of all preprints. Annual files are stored separately due to file 
# size limits of github - here we read them all and reconstruct a single data frame
preprint_abstracts <- map_dfr(2013:2020, ~ read_csv(str_c("data/preprint_abstracts_", .x, ".csv")))

# More detailed preprint metadata covering period Sept 2019 - April 2020
preprints <- read_csv("data/preprints_full_20190101_20201031.csv",
                      col_types = cols(n_refs = "n", n_words = "n"))

# Preprint usage data (PDF downloads, abstract views)
preprint_usage <- read_csv("data/preprint_usage_20190101_20201031.csv")

# Preprint citation data
preprint_citations <- read_csv("data/preprint_citations_20190101_20201031.csv")

# Preprint altmetrics data
preprint_altmetrics <- read_csv("data/preprint_altmetrics_20190101_20201031.csv")

# Data of preprint timing (i.e. timing of first preprint per country - see Figure 2)
preprint_timing <- read_csv("data/preprint_earliest_per_country_20200101_20201031.csv")

# Preprint comments (counts) data
preprint_comments <- read_csv("data/preprint_comments_20190101_20201031.csv")

# Full dataset of preprints from other servers (note: dataset is generated
# externally and updated weekly)
preprints_other <- read_csv("https://raw.githubusercontent.com/nicholasmfraser/covid19_preprints/master/data/covid19_preprints.csv", col_types = "cccccc")

# Usage data from selected other servers (data originally provided by respective servers)
other_server_dloads <- read_csv("data/other_preprint_servers_20200101_20201031.csv")

# Dataset of journal articles on COVID-19 from Dimensions
journal_articles <- read_csv("data/journal_articles_20200101_20201031.csv")

# Policy documents
policy_docs <- read_csv("data/policy_docs.csv")

```

# Figure 1: Development of COVID-19, journal articles and preprints over time

```{r}

# Panel A: Development of COVID-19 cases and deaths

# Define  coordinates for annotations
x_1 <- ymd("2020-01-30")
x_2 <- ymd("2020-03-11")
y <- 50000000

p1A <- covid %>%
  filter(
    Date >= analysis_start,
    Date <= analysis_end) %>%
  select(Date, Confirmed, Deaths) %>%
  rename(Cases = Confirmed) %>%
  pivot_longer(Cases:Deaths) %>%
  ggplot(aes(x = Date, y = value, color = name)) +
  geom_line(size = 1) +
  geom_point(size = 1, shape = 21, fill = "white", stroke = 1) +
  geom_vline(aes(xintercept = x_1), size = 0.5, linetype = "dotted") +
  geom_vline(aes(xintercept = x_2), size = 0.5, linetype = "dotted") +
  geom_text(aes(x = x_1, y = y), label = " (i)", hjust = 0, color = "black") +
  geom_text(aes(x = x_2, y = y), label = " (ii)", hjust = 0, color = "black") +
  labs(x = "Date", y = "Confirmed Cases & Deaths", color = "") +
  scale_y_log10(labels = scales::comma) +
  custom_month_scale(analysis_start, analysis_end + 1) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Dynamic")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_1A.png", width = 10, height = 4)

# Panel B: Number of journal articles and preprints on COVID-19

p1B <- preprints_other %>%
  mutate(posted_date = ymd(posted_date), 
         type = "Preprints",
         clean_title = str_squish(tolower(title))) %>%
  # Remove duplicates occuring in the same source by exact title matching
  group_by(source, clean_title) %>%
  slice_min(order_by = posted_date, n = 1, with_ties = F) %>%
  ungroup() %>%
  bind_rows(journal_articles %>% 
              mutate(posted_date = cr_created,
                     type = "Journal Articles",
                     clean_title = str_squish(tolower(title))) %>%
              # Remove duplicates occuring in the same source by exact title matching
              group_by(clean_title) %>%
              slice_min(order_by = posted_date, n = 1, with_ties = F) %>%
              ungroup()) %>%
  count(type, posted_date) %>%
  complete(posted_date, nesting(type), fill = list(n = 0)) %>%
  group_by(type) %>%
  mutate(cumulative_n = cumsum(n),
         total_n = sum(n)) %>%
  ggplot() +
  geom_area(aes(x = posted_date, y = cumulative_n, fill = reorder(type, -total_n)), 
            color = "grey50", size = 0.5, alpha = 0.75) +
  labs(x = "Publication date", y = "Articles & Preprints", fill = "") +
  scale_y_continuous(labels = scales::comma) +
  custom_month_scale(analysis_start, analysis_end + 1) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Cold")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_1B.png", width = 10, height = 4)

# Panel C: Preprints on COVID-19 from multiple sources

# Find the repositories that have > X preprints
n_min <- 160

top_repositories <- preprints_other %>%
  filter(
    posted_date >= analysis_start,
    posted_date <= analysis_end) %>%
  count(source) %>%
  filter(n > n_min) %>%
  arrange(desc(n)) %>%
  pull(source)

fct_order <- c(top_repositories, "Other")

p1C <- preprints_other %>%
  mutate(posted_date = ymd(posted_date),
         source = case_when(
           source %in% top_repositories ~ source,
           T ~ "Other"),
         source = factor(source, levels = fct_order),
         clean_title = str_squish(tolower(title))) %>%
  # Remove duplicates occuring in the same source by exact title matching
  group_by(source, clean_title) %>%
  slice_min(order_by = posted_date, n = 1, with_ties = F) %>%
  ungroup() %>%
  count(source, posted_date) %>%
  complete(posted_date, nesting(source), fill = list(n = 0)) %>%
  group_by(source) %>%
  mutate(cumulative_n = cumsum(n),
         total_n = sum(n)) %>%
  ggplot() +
  geom_area(aes(x = posted_date, y = cumulative_n, fill = source), 
            color = "grey50", size = 0.5) +
  labs(x = "Publication date", y = "Preprints", fill = "") +
  scale_y_continuous(labels = scales::comma) +
  custom_month_scale(analysis_start, analysis_end + 1) +
  scale_fill_manual(values = qualitative_hcl(16, palette = "Set2")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_1C.png", width = 10, height = 4)

# Patchwork
p1A + p1B + p1C + plot_layout(ncol = 1) +
  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_1.png", width = 10, height = 10)

```

# Supplementary Figure 1: Number of preprints in relation to previous epidemics (Zika, Ebola)

```{r}

# Build search strings containing terms related to other epidemics
search_string_ebola <- "ebola|ZEBOV"
search_string_zika <- "zika|ZIKV"

ebola_start <- "2014-01-24"
ebola_end <- "2016-06-09"  # Outbreak declared ended: https://apps.who.int/iris/bitstream/handle/10665/208883/ebolasitrep_10Jun2016_eng.pdf

zika_start <- "2015-03-02"
zika_end <- "2016-11-18"  # Public health emergency declared ended: https://www.who.int/en/news-room/detail/18-11-2016-fifth-meeting-of-the-emergency-committee-under-the-international-health-regulations-(2005)-regarding-microcephaly-other-neurological-disorders-and-zika-virus

# Determine preprint counts of previous epidemics
covid_counts <- preprints_all %>% 
  filter(posted_date >= analysis_start & posted_date <= analysis_end) %>% 
  count(covid_preprint) %>%
  mutate(epidemic = "COVID-19") %>%
  rename(epi_preprint = covid_preprint)

ebola_counts <- preprints_all %>%
  inner_join(preprint_abstracts, by = "doi") %>%
  filter(posted_date >= ebola_start & posted_date <= ebola_end) %>%
  mutate(epi_preprint = case_when(
    str_detect(title, regex(search_string_ebola, ignore_case = TRUE)) ~ T,
    str_detect(abstract, regex(search_string_ebola, ignore_case = TRUE)) ~ T,
    T ~ F)) %>% 
  count(epi_preprint) %>%
  mutate(epidemic = "Western Africa Ebola virus")

zika_counts <- preprints_all %>% 
  inner_join(preprint_abstracts, by = "doi") %>%
  filter(posted_date >= zika_start & posted_date <= zika_end) %>%
  mutate(epi_preprint = case_when(
    str_detect(title, regex(search_string_zika, ignore_case = TRUE)) ~ T,
    str_detect(abstract, regex(search_string_ebola, ignore_case = TRUE)) ~ T,
    T ~ F)) %>% 
  count(epi_preprint) %>%
  mutate(epidemic = "Zika virus")

# Generate plot
pS1A <- bind_rows(covid_counts, ebola_counts, zika_counts) %>% 
  mutate(epi_preprint = case_when(
    epi_preprint == T ~ "Epidemic preprints",
    epi_preprint == F ~ "Non-epidemic preprints"),
    epidemic = factor(epidemic, 
                      levels = c("Western Africa Ebola virus", 
                                 "Zika virus", 
                                 "COVID-19"),
                      labels = c("Western Africa Ebola virus\n(2014/04/24 - 2016/06/09)", 
                                 "Zika virus\n(2015/03/02 - 2016/11/18)", 
                                 "COVID-19\n(2020/01/01 - 2020/10/31)"))) %>%
  ggplot(aes(x = epi_preprint, y = n, colour = epi_preprint, fill = epi_preprint)) +
  geom_col(color = "grey50", size = 0.25) +
  labs(x = "", y = "Total preprints", fill = "", color = "") +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  facet_wrap(epidemic ~ ., strip.position = "bottom") +
  theme(axis.text.x = element_blank(),
        plot.margin = margin(0, 0, 20, 0)) +
  ggsave("outputs/figures/partials/Figure_S1A.png", width = 10, height = 4)

pS1B <-journal_articles %>% 
  mutate(clean_title = str_squish(tolower(title))) %>%
  # Remove duplicates occuring in the same source by exact title matching
  group_by(clean_title) %>%
  slice_min(order_by = cr_created, n = 1, with_ties = F) %>%
  ungroup() %>%
  count(article_is_oa, cr_created) %>%
  complete(cr_created, nesting(article_is_oa), fill = list(n = 0)) %>%
  group_by(article_is_oa) %>%
  mutate(cumulative_n = cumsum(n),
         total_n = sum(n)) %>%
  ungroup() %>%
  mutate(article_is_oa = case_when(
    article_is_oa == T ~ "Open Access",
    T  ~ "Closed Access")) %>%
  ggplot() +
  geom_area(aes(x = cr_created, y = cumulative_n, 
                fill = reorder(article_is_oa, -total_n)), 
            color = "grey50", size = 0.5) +
  labs(x = "Publication Date", y = "Published articles", fill = "") +
  scale_y_continuous(labels = scales::comma) +
  custom_month_scale("2020-01-01", "2020-11-01", format = "%d-%b") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Dynamic")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_S1B.png", width = 10, height = 4)

pS1A + pS1B + plot_layout(ncol = 1) +
  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_S1.png", width = 10, height = 6)

```

# Figure 2: Preprint properties

```{r}

# Panel A: COVID vs non-COVID preprints deposited for each server
p2A <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  assign_covid_preprint() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  count(covid_preprint, posted_month) %>%
  complete(posted_month, nesting(covid_preprint), fill = list(n = 0)) %>%
  ggplot(aes(x = posted_month, y = n, fill = covid_preprint)) +
  geom_col(color = "grey50", size = 0.25,
           position = position_dodge(preserve = "single")) +
  labs(x = "Posted month", y = "Preprints deposited", fill = "") +
  scale_y_continuous(labels = scales::comma) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  ggsave("outputs/figures/partials/Figure_2A.png", width = 6, height = 4)

# Panel B: Screening time
p2B <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(doi_date = gsub("\\.", "-", substr(doi, 9, 18)),
         days = as.numeric(ymd(posted_date) - ymd(doi_date)),
         days = case_when(
           days > 7 ~ "8+",
           T ~ as.character(days)),
         days = factor(days,
                       levels = c("0", "1", "2", "3", "4", "5", "6", "7", "8+"),
                       labels = c("0", "1", "2", "3", "4", "5", "6", "7", "8+"))
  ) %>%
  count(covid_preprint, days) %>%
  group_by(covid_preprint) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = days, y = proportion,  fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), 
           width = 0.75, size = 0.25) +    
  labs(x = "Screening time (days)", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  guides(fill = FALSE) +
  ggsave("outputs/figures/partials/Figure_2B.png", width = 5, height = 4)

# Panel C: Versions
p2C <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  count(covid_preprint, n_versions) %>%
  group_by(covid_preprint) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = n_versions, y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), width = 0.75, size = 0.25) +
  labs(x = "Number of versions", y = "% of preprints", fill = "") +
  scale_x_continuous(breaks = seq(1:10)) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  guides(fill = FALSE) +
  ggsave("outputs/figures/partials/Figure_2C.png", width = 5, height = 4)

# Panel D: Licenses
p2D <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(license = case_when(
    str_detect(license, "cc0") ~ "cc0",
    T ~ license
  )) %>%
  count(covid_preprint, license) %>%
  na.omit() %>%
  group_by(covid_preprint) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  mutate(license = factor(license,
                          levels = c("cc0", "cc_by", "cc_by_nc", "cc_by_nd", 
                                     "cc_by_nc_nd", "cc_no", "cc0_ng"),
                          labels = c("CC0", "CC BY", "CC BY-NC", "CC BY-ND", 
                                     "CC BY-NC-ND", "None", "CC0 NG?"))) %>%
  ggplot(aes(x = license, y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), width = 0.75, size = 0.25) +
  labs(x = "License Type", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_discrete(label = function(x) stringr::str_wrap(x, 6)) +
  guides(fill = FALSE) +
  ggsave("outputs/figures/partials/Figure_2D.png", width = 5, height = 4)

# Panel E: Word counts
p2E <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = posted_month, y = n_words, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1, 
              position = position_jitterdodge(jitter.width = 10, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "Word counts", 
       fill = "", color = "") +
  scale_y_continuous(labels = scales::comma, expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 15000)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_2E.png", width = 10, height = 4)

# Panel F: Reference counts
p2F <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = posted_month, y = n_refs, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1,  
              position = position_jitterdodge(jitter.width = 10, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "Reference counts", 
       fill = "", color = "") +
  scale_y_continuous(labels = scales::comma, expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 150)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_2F.png", width = 10, height = 4)

p <- p2A + p2B + p2C + p2D + p2E + p2F +
  plot_layout(ncol = 2, guides = 'collect') + 
  plot_annotation(tag_levels = "A") & theme(legend.position = 'bottom')

p + ggsave("outputs/figures/Figure_2.png", width = 10, height = 10) 

```

# Supplementary Figure 2: Preprint properties

```{r}

# Panel A: COVID vs non-COVID preprints for each server
pS2A <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  assign_covid_preprint() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month"),
         source = factor(source, 
                         levels = c("biorxiv", "medrxiv"),
                         labels = c("bioRxiv", "medRxiv"))) %>%
  count(covid_preprint, source, posted_month) %>%
  complete(posted_month, nesting(covid_preprint, source), fill = list(n = 0)) %>%
  ggplot(aes(x = posted_month, y = n, fill = covid_preprint)) +
  geom_col(color = "grey50", size = 0.25,
           position = position_dodge(preserve = "single")) +
  labs(x = "Posted month", y = "Preprints deposited", fill = "") +
  scale_y_continuous(labels = scales::comma) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  facet_wrap(. ~ source, nrow = 2, strip.position = "right") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_rect(color = "grey90", fill = "grey95"),
        legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_S2A.png", width = 6, height = 4)

# Panel B: Screening time
pS2B <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(doi_date = gsub("\\.", "-", substr(doi, 9, 18)),
         days = as.numeric(ymd(posted_date) - ymd(doi_date)),
         days = case_when(
           days > 7 ~ "8+",
           T ~ as.character(days)),
         days = factor(days,
                       levels = c("0", "1", "2", "3", "4", "5", "6", "7", "8+"),
                       labels = c("0", "1", "2", "3", "4", "5", "6", "7", "8+"))
  ) %>%
  count(covid_preprint, source, days) %>%
  group_by(covid_preprint, source) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = days, y = proportion,  fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), 
           width = 0.75, size = 0.25) +    
  labs(x = "Screening time (days)", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  guides(fill = FALSE) +
  facet_wrap(. ~ source, nrow = 2, strip.position = "right") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_rect(color = "grey90", fill = "grey95"),
        legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_S2B.png", width = 5, height = 4)

# Panel C: Versions
pS2C <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  count(covid_preprint, source, n_versions) %>%
  group_by(covid_preprint, source) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = n_versions, y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), width = 0.75, size = 0.25) +
  labs(x = "Number of versions", y = "% of preprints", fill = "") +
  scale_x_continuous(breaks = seq(1:10)) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  guides(fill = FALSE) +
  facet_wrap(. ~ source, nrow = 2, strip.position = "right") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_rect(color = "grey90", fill = "grey95"),
        legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_S2C.png", width = 5, height = 4)

# Panel D: Licences
pS2D <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(license = case_when(
    str_detect(license, "cc0") ~ "cc0",
    T ~ license
  )) %>%
  count(covid_preprint, source, license) %>%
  na.omit() %>%
  group_by(covid_preprint, source) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  mutate(license = factor(license,
                          levels = c("cc0", "cc_by", "cc_by_nc", "cc_by_nd", 
                                     "cc_by_nc_nd", "cc_no", "cc0_ng"),
                          labels = c("CC0", "CC BY", "CC BY-NC", "CC BY-ND", 
                                     "CC BY-NC-ND", "None", "CC0 NG?"))) %>%
  ggplot(aes(x = license, y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), width = 0.75, size = 0.25) +
  labs(x = "License Type", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_discrete(label = function(x) stringr::str_wrap(x, 6)) +
  guides(fill = FALSE) +
  facet_wrap(. ~ source, nrow = 2, strip.position = "right") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        strip.background = element_rect(color = "grey90", fill = "grey95"),
        legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_S2D.png", width = 5, height = 4)

pS2 <- pS2A + pS2B + pS2C + pS2D +
  plot_layout(ncol = 2, guides = 'collect') + 
  plot_annotation(tag_levels = "A") & theme(legend.position = 'bottom')

pS2 + ggsave("outputs/figures/Figure_S2.png", width = 10, height = 7) 

```

# Figure 3: Authorships

```{r}

# Panel A: Number of authors per preprint
p3A <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  mutate(n_authors = case_when(
    n_authors > 10 ~ "11+",
    T ~ as.character(n_authors)),
    n_authors = factor(n_authors,
                       levels = c(as.character(1:10), "11+"),
                       labels = c(as.character(1:10), "11+"))) %>%
  count(covid_preprint, n_authors) %>%
  group_by(covid_preprint) %>%
  mutate(proportion = (n/sum(n)) * 100) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = n_authors, y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), 
           width = 0.75, size = 0.25) +
  labs(x = "Number of authors", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  theme(legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_3A.png", width = 5, height = 4)

# Panel B: Top countries of authors
top_n_countries <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end,
         !is.na(institution_match_country_code)) %>%
  count(institution_match_country_code) %>%
  slice_max(order_by = n, n = 15) %>% # select top 15
  pull(institution_match_country_code)

p3B <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end,
         institution_match_country_code %in% top_n_countries) %>%
  count(institution_match_country_code, covid_preprint) %>%
  group_by(covid_preprint) %>%
  mutate(total_n = sum(n),
         proportion = (n/sum(n))*100) %>%
  ungroup() %>%
  mutate(continent = countrycode::countrycode(institution_match_country_code, 
                                              origin = 'iso2c', 
                                              destination = 'continent')) %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = reorder(institution_match_country_code, -n), y = proportion, fill = covid_preprint)) +
  geom_col(color = "grey50", position = position_dodge(preserve = "single"), 
           width = 0.75, size = 0.25) +
  labs(x = "Country of corresponding author", y = "% of preprints", fill = "") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  theme(legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_3B.png", width = 5, height = 4)

# Panel C - % senior authors vs country

# Authors of COVID-19 preprints
covid_authors <- preprints %>% 
  filter(covid_preprint == T) %>%
  clean_names() %>%
  pull(author_corresponding)

# Control authors (authors not in COVID-19 group)
control_authors <- preprints %>% 
  filter(covid_preprint == F & !(author_corresponding %in% covid_authors)) %>%
  clean_names() %>%
  pull(author_corresponding)

# Determine preprinter history
preprinter_history <- preprints_all %>%
  filter(posted_date < analysis_start,
         author_corresponding != "? ?",
         !is.na(author_corresponding)) %>%
  clean_names() %>%
  mutate(author_corresponding = factor(author_corresponding, 
                                       levels = unique(c(covid_authors, control_authors)))) %>%
  group_by(author_corresponding, .drop = F) %>%
  summarize(n = n()) %>%
  mutate(
    author_group = case_when(
      author_corresponding %in% covid_authors ~ "COVID-19 preprints",
      author_corresponding %in% control_authors ~ "non-COVID-19 preprints"
    ),
    preprinter_status = case_when(
      n == 0 ~ "First-time", 
      T ~ "Previous")
  ) %>%
  distinct(author_corresponding, author_group, preprinter_status) %>%
  filter(author_corresponding %in% (preprints %>% 
                                      filter(posted_date >= analysis_start,
                                             posted_date <= analysis_end) %>% 
                                      clean_names() %>%
                                      pull(author_corresponding))) # Limit to authors during the study period only

# Append author locations and calculate proportions for plot
preprinter_history_tabular <- preprints %>%
  clean_names() %>%
  inner_join(preprinter_history, by = "author_corresponding") %>%
  distinct(author_corresponding, author_group, preprinter_status, institution_match_country_code) %>%
  count(author_group, preprinter_status, institution_match_country_code) %>% 
  complete(author_group, preprinter_status, institution_match_country_code, 
           fill = list(n = 0))

countries_to_plot <- preprinter_history_tabular %>% 
  filter(!is.na(institution_match_country_code)) %>% 
  group_by(institution_match_country_code) %>% 
  summarise(n = sum(n)) %>% 
  top_n(15, n) %>% 
  arrange(-n) %>%
  pull(institution_match_country_code)

preprinter_history_tabular <- preprinter_history_tabular %>% filter(institution_match_country_code %in% countries_to_plot)

preprinter_history_proportions <- preprinter_history_tabular %>%
  group_by(author_group, institution_match_country_code) %>%
  mutate(prop = n*100/sum(n)) 

p3C <- preprinter_history_proportions %>%
  ggplot() +
  geom_bar(mapping = aes(x = institution_match_country_code, y = prop, 
                         fill = author_group, alpha = preprinter_status), 
           stat = "identity", position = "stack", color = "grey50", size = 0.25) + 
  coord_flip() +
  scale_x_discrete(limits = rev(countries_to_plot)) +
  scale_y_continuous(expand = c(0,0)) +
  facet_wrap(facets = ~author_group, scales = "free") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_alpha_manual(values = c(1.0, 0.4),
                     breaks = c("Previous", "First-time")) +
  labs(x = "Country", y = "Proportion corresponding authors") +
  guides(fill = FALSE) +
  theme(legend.title = element_blank(),
        legend.margin=margin(c(0,0,0,0)),
        legend.position = "top",
        plot.margin = margin(0,10,0,0),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  ggsave("outputs/figures/partials/Figure_3C.png", width = 8, height = 5)

# Panel D
p3D <- preprint_timing %>%
  mutate(continent = countrycode::countrycode(institution_match_country_code, 
                                              origin = 'iso2c', destination = 'continent')) %>% 
  filter(!is.na(first_preprint) & !is.na(first_case) & !is.na(continent)) %>%
  filter(first_case <= "2020-04-30",
         first_preprint <= "2020-04-30") %>%
  ggplot(aes(x = first_case, y = first_preprint, label = institution_match_country_code)) +
  geom_point(aes(color = continent), size=2) +
  geom_text_repel(force = 5,
                  segment.color = "grey50",
                  segment.size = 0.2,
                  segment.alpha = 0.8,
                  color = "grey25", 
                  size = 2.5) +
  geom_abline(intercept = 0, slope = 1, color = "grey50", lty = "dashed", alpha = 0.5) +
  labs(x = "Date of first case",
       y = "Date of first preprint",
       color = "") +
  scale_color_manual(values = qualitative_hcl(5, palette = "Set2")) +
  coord_cartesian(ylim = c(ymd("2020-01-01", "2020-04-30")),
                  xlim = c(ymd("2020-01-01", "2020-04-30"))) +
  theme(legend.margin = margin(0, 0, 0, 0),
        legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_3D.png", width = 5, height = 4)

# Panel E: Science focus of senior author

# For each author, select their chronologically earliest COVID-19 preprint
preprint_author_covid <- preprints %>%
  clean_names() %>%
  filter(covid_preprint == T) %>% 
  group_by(author_corresponding) %>%
  arrange(posted_date) %>%
  slice(1) %>%
  select(author_corresponding, category, posted_date) %>%
  rename(category_covid = category,
         posted_covid = posted_date)

# For each author, select all their non-COVID-19 preprints
preprint_author_non_covid <- preprints_all %>%
  clean_names() %>%
  filter(covid_preprint == F) %>% 
  select(author_corresponding, category, posted_date) %>%
  rename(category_previous = category,
         posted_previous = posted_date)

# Filter to non-COVID-19 preprints uploaded before the latest COVID-19 preprint and select the chronologically latest
preprint_author_last <- preprint_author_covid %>%
  inner_join(preprint_author_non_covid, by = "author_corresponding") %>%
  filter(posted_covid >= posted_previous) %>%
  group_by(author_corresponding) %>%
  arrange(posted_previous) %>%
  slice(n()) %>%
  mutate(category_covid = tolower(category_covid),
         category_previous = tolower(category_previous))

total_cats <- c(preprint_author_last$category_covid, preprint_author_last$category_previous) %>% unique %>% sort

# Plot alluvial diagram
p3E <- preprint_author_last %>%
  mutate(category_covid = factor(category_covid, levels = total_cats),
         category_previous = factor(category_previous, levels = total_cats)) %>%
  group_by(category_covid, category_previous) %>%
  summarize(freq = n()) %>%
  ungroup() %>%
  filter(freq >= 10) %>%
  ggplot(aes(axis1 = category_previous, axis2 = category_covid, y = freq)) +
  scale_x_discrete(limits = c("Previous", "COVID-19"), 
                   expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  geom_alluvium(aes(fill = category_previous), reverse = TRUE) +
  geom_stratum(reverse = TRUE, color = "grey50") + 
  geom_text(aes(label = after_stat(stratum)), size = 2.5, stat = "stratum", 
            reverse = TRUE, color = "grey25") +
  guides(fill = FALSE) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  ggsave("outputs/figures/partials/Figure_3E.png", width = 5, height = 8)

layout <- "
AB
CE
DE
"

p <- p3A + p3B + p3C + p3D + p3E +
  plot_layout(design = layout, heights = c(0.5, 0.75, 1))

p + plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_3.png", width = 10, height = 12) 

```

# Figure 4: Publication outcomes

```{r}

# Panel A: proportions of preprints published over time
p4A <- preprints %>% 
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  assign_covid_preprint() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month"),
         is_published = !is.na(published_doi)) %>%
  count(posted_month, covid_preprint, is_published) %>%
  group_by(posted_month, covid_preprint) %>%
  mutate(prop = n*100 / sum(n)) %>%
  ungroup() %>%
  filter(is_published == T) %>%
  ggplot(aes(x = posted_month, y = prop, fill = covid_preprint)) +
  geom_col(color = "grey50", size = 0.25, position = "dodge") +
  labs(x = "", y = "% of preprints published", fill ="") +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  theme(legend.position = "top") +
  ggsave("outputs/figures/partials/Figure_4A.png", width = 5, height = 4)

# Panel B: Journal locations of published papers (top 10)
p4B <- preprints %>%
  filter(covid_preprint == T) %>%
  count(published_journal) %>%
  na.omit() %>%
  arrange(-n) %>%
  slice(1:10) %>%
  ggplot(aes(x = reorder(published_journal, -n), y = n)) +
  geom_col(color = "grey50", fill = qualitative_hcl(n = 1, palette = "Set2"),
           position = "dodge", size = 0.25) +
  labs(x = "", y = str_wrap("Number of COVID-19 preprints published", 20), fill = "") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 25)) +
  theme(axis.text.x = element_text(size = 8, angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_4B.png", width = 5, height = 4) 

# Panel C: time to publication
p4C <- preprints %>%
  mutate(covid_preprint = case_when(
    (covid_preprint == T & posted_date >= analysis_start & posted_date <= analysis_end) ~ "COVID-19 preprints",
    (covid_preprint == F & posted_date >= analysis_start & posted_date <= analysis_end) ~ "non-COVID-19 preprints",
    T ~ "preprints 2019")) %>%
  mutate(pub_bracket = cut(as.numeric(delay_in_days), 
                           seq(0, 360, by = 10), 
                           labels=seq(0, 350, by = 10))) %>%
  group_by(covid_preprint) %>%
  count(pub_bracket) %>%
  mutate(prop = n*100 / sum(n)) %>%
  filter(!is.na(pub_bracket)) %>%
  ggplot(aes(x = pub_bracket, y = prop, fill=covid_preprint, color=covid_preprint)) +
  geom_bar(alpha = 0.25, width = 1, size = 0.25, 
           stat = "identity", position="identity") +
  scale_x_discrete(breaks = as.character(seq(0, 350, by = 50))) +
  labs(x = "Time from preprint posting to publication (days)",
       y = "% of all preprints") +
  scale_color_manual(values = c(qualitative_hcl(2, palette = "Set2"), "grey60")) +
  scale_fill_manual(values = c(qualitative_hcl(2, palette = "Set2"), "grey60")) +
  theme(legend.title = element_blank(),
        legend.position = "top",
        legend.text = element_text(size = 8)) +
  ggsave("outputs/figures/partials/Figure_4C.png", width = 5, height = 4)

# Panel D: Time to publication by publisher
top_publishers <- preprints %>%
  filter(covid_preprint == T,
         posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  count(published_publisher) %>%
  na.omit() %>%
  top_n(10, n) %>%
  pull(published_publisher)

p4D <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end,
         published_publisher %in% top_publishers) %>%
  assign_covid_preprint() %>%
  mutate(published_publisher = factor(published_publisher,
                                      levels = c("American Association for the Advancement of Science (AAAS)",
                                                 "American Society for Microbiology",
                                                 "BMJ",
                                                 "Elsevier BV",
                                                 "Frontiers Media SA",
                                                 "MDPI AG",
                                                 "Oxford University Press (OUP)",
                                                 "Public Library of Science (PLoS)",
                                                 "Springer Science and Business Media LLC",
                                                 "Wiley"),
                                      labels = c("AAAS", "ASM", "BMJ", "Elsevier", "Frontiers", 
                                                 "MDPI", "OUP", "PLoS", "Springer", "Wiley"))) %>%
  filter(delay_in_days > 0) %>%
  ggplot(aes(x = published_publisher, y = delay_in_days, color = covid_preprint)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.2, alpha = 0.2,
              position = position_jitterdodge(jitter.width = 0.5)) + 
  labs(x = "", y = str_wrap("Time from preprint posting to publication (days)", 30), 
       fill = "", color = "") +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        legend.position = "top") +
  guides(fill = FALSE, colour = FALSE) +
  ggsave("outputs/figures/partials/Figure_4D.png", width = 5, height = 4) 


p <- p4A + p4B + p4C + p4D +
  plot_layout(ncol = 2)

p + plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_4.png", width = 10, height = 7) 

```

# Figure 5: Preprint access

```{r}

# Prepare data
d <- preprint_usage %>%
  inner_join(preprints, by = "doi") %>%
  group_by(doi, posted_date, covid_preprint) %>%
  summarize(full_text_views = sum(full_text_views),
            abstract_views = sum(abstract_views),
            pdf_downloads = sum(pdf_downloads)) %>%
  ungroup() %>%
  assign_covid_preprint() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  filter(posted_month >= analysis_start,
         posted_month <= analysis_end)

# Abstract views
p5A <- d %>%
  ggplot(aes(x = posted_month, y = abstract_views, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1, 
              position = position_jitterdodge(jitter.width = 15, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "Total abstract views", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e6), expand = c(0, 0)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_5A.png", width = 10, height = 4)

# PDF downloads
p5B <- d %>%
  ggplot(aes(x = posted_month, y = pdf_downloads, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1, 
              position = position_jitterdodge(jitter.width = 15, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "Total PDF downloads", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e6), expand = c(0, 0)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_5B.png", width = 10, height = 4)

# Patchwork
p5A + p5B + plot_layout(ncol = 1) +
  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_5.png", width = 10, height = 8)

```

# Supplementary Figure 3: Preprint access in time prior to our study period

```{r}

# Panel A: Abstract views in first month after posting
d <- preprint_usage %>%
  mutate(
    collection_month = str_sub(collection_date, 1, 3),
    collection_year = str_c("20", str_sub(collection_date, 
                                          nchar(collection_date) - 1,
                                          nchar(collection_date))),
    collection_date = my(str_c(collection_month, collection_year))) %>%
  inner_join(preprints, by = "doi") %>%
  assign_covid_preprint() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  filter(posted_month >= analysis_start,
         posted_month <= analysis_end,
         posted_month == collection_date)

pS3A <- d %>%
  ggplot(aes(x = posted_month, y = abstract_views, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1, 
              position = position_jitterdodge(jitter.width = 15, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "Abstract views in posting month", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e6), expand = c(0, 0)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_S3A.png", width = 10, height = 4)

# Panel B: PDF downloads
pS3B <- d %>%
  ggplot(aes(x = posted_month, y = pdf_downloads, color = covid_preprint,
             group = interaction(posted_month, covid_preprint))) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.25, alpha = 0.1, 
              position = position_jitterdodge(jitter.width = 15, dodge.width = 21)) +
  geom_boxplot(width = 20, outlier.shape = NA, fill = NA,
               position = position_dodge(width = 21)) +
  labs(x = "Posted Month", y = "PDF downloads in posting month", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e6), expand = c(0, 0)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b") +
  ggsave("outputs/figures/partials/Figure_S3B.png", width = 10, height = 4)

# Panel C: Abstract views (bioRxiv + medRxiv, non-COVID)
d <- preprint_usage %>%
  inner_join(preprints, by = "doi") %>%
  filter(covid_preprint != T) %>%
  group_by(doi, posted_date) %>%
  summarize(full_text_views = sum(full_text_views),
            abstract_views = sum(abstract_views),
            pdf_downloads = sum(pdf_downloads)) %>%
  ungroup() %>%
  mutate(posted_month = floor_date(posted_date, unit = "month")) %>%
  filter(posted_date >= "2019-01-01",
         posted_date <= analysis_end)

pS3C <- d %>%
  ggplot(aes(x = posted_month, y = abstract_views, group = posted_month)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(shape = 21, size = 0.2, alpha = 0.1) +
  annotate("rect", xmin = ymd("2019-12-15"), xmax = ymd("2020-10-15"),
           ymin = 1, ymax = 1e5, fill = "red", alpha = 0.1) +
  labs(x = "Posted Month", y = "Total abstract views", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e5), expand = c(0, 0)) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b-%y",
               expand = c(0.01,0.01)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_S3C.png", width = 10, height = 4)

# Panel D: PDF downloads (bioRxiv + medRxiv, non-COVID)
pS3D <- d %>%
  ggplot(aes(x = posted_month, y = pdf_downloads, group = posted_month)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(shape = 21, size = 0.2, alpha = 0.1) +
  annotate("rect", xmin = ymd("2019-12-15"), xmax = ymd("2020-10-15"),
           ymin = 1, ymax = 1e5, fill = "red", alpha = 0.1) +
  labs(x = "Posted Month", y = "Total PDF downloads", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e5), expand = c(0, 0)) +
  scale_x_date(date_breaks = "1 month",
               date_minor_breaks = "1 month",
               date_labels = "%b-%y",
               expand = c(0.01,0.01)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_S3D.png", width = 10, height = 4)

# Panel E: PDF downloads for additional preprint servers
# Bind into one data frame
all_server_dloads <- bind_rows(preprint_usage %>%
                                 inner_join(preprints, by = c("doi", "source")) %>%
                                 group_by(source, doi, posted_date, covid_preprint) %>%
                                 summarize(pdf_downloads = sum(pdf_downloads)) %>%
                                 ungroup() %>%
                                 select(-doi),
                               other_server_dloads %>%
                                 mutate(posted_date = as.Date(posted_date, "%d/%m/%Y"))) %>%
  filter(posted_date >= analysis_start & posted_date <= analysis_end)

# Generate plot
pS3E <- all_server_dloads %>%
  assign_covid_preprint() %>%
  ggplot(aes(x = source, y = pdf_downloads, color = covid_preprint)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(fill = factor(covid_preprint)),
              shape = 21, size = 0.2, alpha = 0.2,
              position = position_jitterdodge(jitter.width = 0.5)) + 
  labs(x = "Preprint server", y = "Total PDF downloads", 
       fill = "", color = "") +
  scale_y_log10(labels = scales::comma, limits = c(1, 1e7), expand = c(0, 0)) +
  scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
  scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
  ggsave("outputs/figures/partials/Figure_S3E.png", width = 10, height = 4)

# Patchwork
pS3A + pS3B + pS3C + pS3D + pS3E +
  plot_layout(nrow = 5) +
  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_S3.png", width = 10, height = 14)

```

# Figure 6: Preprint usage and sharing

```{r}

d <- preprints %>%
  filter(posted_date >= analysis_start,
         posted_date <= analysis_end) %>%
  inner_join(preprint_citations, by = "doi") %>%
  inner_join(preprint_altmetrics, by = "doi") %>%
  inner_join(preprint_comments, by = "doi") %>%
  mutate(posted_week = floor_date(posted_date, unit = "week", 
                                  week_start = 1)) %>%
  assign_covid_preprint()

plot_main <- function(d, title_y, limit_y) {
  
  d %>%
    group_by(covid_preprint, posted_week) %>%
    summarise(prop_cited = sum(is_cited)/n()) %>%
    ggplot(aes(x = posted_week, y = prop_cited, 
               color = covid_preprint, fill = covid_preprint)) +
    geom_line() +
    geom_point(shape = 21) +
    labs(x = "Date posted", y = title_y, fill = "", color = "") +
    scale_y_continuous(labels = scales::percent, limits = c(0,limit_y)) +
    scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
    scale_fill_manual(values = qualitative_hcl(2, palette = "Set2"))  +
    guides(fill = FALSE) +
    theme(axis.text.y = element_text(size = 8),
          axis.title.y = element_text(size = 8))
}

plot_inset <- function(d, title_y) {
  
  d %>%
    ggplot(aes(x = covid_preprint, y = metric + 1, color = covid_preprint)) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(aes(fill = factor(covid_preprint)), shape = 21, 
                size = 0.25, alpha = 0.2, width = 0.25) +
    labs(x = "", y = title_y, 
         fill = "", color = "") +
    scale_y_log10(labels = scales::comma_format(accuracy = 1), expand = c(0, 0.1)) +
    scale_color_manual(values = qualitative_hcl(2, palette = "Set2")) +
    scale_fill_manual(values = qualitative_hcl(2, palette = "Set2")) +
    scale_x_discrete(label = function(x) str_wrap(x, 12)) +
    guides(fill = FALSE, color = FALSE) +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_text(size = 7),
          axis.title.y = element_text(size = 7),
          plot.margin = margin(0,0,0,0),
          panel.background = element_rect(fill = "white"),
          panel.border = element_rect(color = "grey75"))
  
}

# Citations
p6A_main <- d %>%
  mutate(is_cited = as.numeric(citations > 0)) %>%
  plot_main(title_y = str_wrap("% of preprints with \u2265 1 citation", 30),
            limit_y = 1)
p6A_inset <- d %>%
  mutate(metric = citations) %>%
  plot_inset(title_y = "Citations")
p6A <- p6A_main + 
  annotation_custom(ggplotGrob(p6A_inset), 
                    xmin = as.Date("2020-06-15"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.4, ymax = 1.04) +
  ggsave("outputs/figures/partials/Figure_6A.png", width = 5, height = 4)

# Twitter
p6B_main <- d %>%
  mutate(is_cited = as.numeric(twitter > 1)) %>%
  plot_main(title_y = str_wrap("% of preprints with \u2265 2 tweets*", 30),
            limit_y = 1)
p6B_inset <- d %>%
  mutate(metric = twitter) %>%
  plot_inset(title_y = "Tweets")
p6B <- p6B_main + 
  annotation_custom(ggplotGrob(p6B_inset), 
                    xmin = as.Date("2020-06-01"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.1, ymax = 0.74) +
  ggsave("outputs/figures/partials/Figure_6B.png", width = 5, height = 4)

# News
p6C_main <- d %>%
  mutate(is_cited = as.numeric(news > 0)) %>%
  plot_main(title_y = str_wrap("% of preprints mentioned in \u2265 1 news article", 30),
            limit_y = 1)
p6C_inset <- d %>%
  mutate(metric = news) %>%
  plot_inset(title_y = "News")
p6C <- p6C_main + 
  annotation_custom(ggplotGrob(p6C_inset), 
                    xmin = as.Date("2020-06-01"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.4, ymax = 1.04) +
  ggsave("outputs/figures/partials/Figure_6C.png", width = 5, height = 4)

# Blogs
p6D_main <- d %>%
  mutate(is_cited = as.numeric(blogs > 0)) %>%
  plot_main(title_y = str_wrap("% of preprints mentioned in \u2265 1 blog post", 30),
            limit_y = 1)
p6D_inset <- d %>%
  mutate(metric = blogs) %>%
  plot_inset(title_y = "Blogs")
p6D <- p6D_main + 
  annotation_custom(ggplotGrob(p6D_inset), 
                    xmin = as.Date("2020-06-01"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.4, ymax = 1.04) +
  ggsave("outputs/figures/partials/Figure_6D.png", width = 5, height = 4)

# Wikipedia
p6E_main <- d %>%
  mutate(is_cited = as.numeric(wikipedia > 0)) %>%
  plot_main(title_y = str_wrap("% of preprints with \u2265 1 Wikipedia citation", 30),
            limit_y = 0.1)
p6E_inset <- d %>%
  mutate(metric = wikipedia) %>%
  plot_inset("Wikipedia")
p6E <- p6E_main + 
  annotation_custom(ggplotGrob(p6E_inset), 
                    xmin = as.Date("2020-06-01"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.04, ymax = 0.104) +
  ggsave("outputs/figures/partials/Figure_6E.png", width = 5, height = 4)

# Comments
p6F_main <- d %>%
  mutate(is_cited = as.numeric(comments_count > 0)) %>%
  plot_main(title_y = str_wrap("% of preprints with \u2265 1 comment", 30),
            limit_y = 1)
p6F_inset <- d %>%
  mutate(metric = comments_count) %>%
  plot_inset(title_y = "Comments")
p6F <- p6F_main + 
  annotation_custom(ggplotGrob(p6F_inset), 
                    xmin = as.Date("2020-06-01"), 
                    xmax = as.Date("2020-11-07"), 
                    ymin = 0.4, ymax = 1.04) +
  ggsave("outputs/figures/partials/Figure_6F.png", width = 5, height = 4)

# Policy docs
p6G <- policy_docs %>%
#  filter(type == "COVID-19") %>%
  mutate(preprint_prop = npreprints*100/nrefs,
         set = case_when(set == "ECDC" ~ "ECDC",
                         set == "National academy of medicine" ~ "NAM",
                         set == "US committee COVID" ~ "US HSSCC",
                         set == "WHO Scientific Briefs" ~ "WHO SB",
                         set == "UK post" ~ "UK POST")) %>%
  ggplot(aes(x = set, y = preprint_prop)) +
  geom_dotplot(fill = qualitative_hcl(1, palette = "Set2"), color = qualitative_hcl(1, palette = "Set2"), binaxis = "y", binwidth = 1,  stackdir = "center", binpositions="all") +
  ylab("% preprints in references") +
  guides(size = guide_legend(title="total\ncitations")) +
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggsave("outputs/figures/partials/Figure_6G.png", width = 4, height = 3)

# PCorrelation plots between citations, tweets, news articles, blogs, wikipedia, and comments
covid_usage <- d %>% 
  inner_join(preprint_usage %>%
               group_by(doi) %>%
               summarise(
                 abstract_views = sum(abstract_views),
                 pdf_downloads = sum(pdf_downloads)
               ), by = "doi") %>%
  filter(covid_preprint == "COVID-19 preprints") %>%
  rename(comments = comments_count) %>%
  select(citations, twitter, news, blogs, wikipedia, comments, abstract_views, pdf_downloads)

non_covid_usage <- d %>% 
  inner_join(preprint_usage %>%
               group_by(doi) %>%
               summarise(
                 abstract_views = sum(abstract_views),
                 pdf_downloads = sum(pdf_downloads)
               ), by = "doi") %>%
  filter(covid_preprint == "non-COVID-19 preprints") %>%
  rename(comments = comments_count) %>%
  select(citations, twitter, news, blogs, wikipedia, comments, abstract_views, pdf_downloads)

covid_corr <- correlate(covid_usage, method = "spearman") %>%
  pivot_longer(citations:pdf_downloads) %>%
  rename(rows = term,
         cols = name) %>%
  filter(str_trunc(rows, 2, ellipsis = "") < str_trunc(cols, 2, ellipsis = ""))

non_covid_corr <- correlate(non_covid_usage, method = "spearman") %>%
  pivot_longer(citations:pdf_downloads) %>%
  rename(rows = term,
         cols = name) %>%
  filter(str_trunc(rows, 2, ellipsis = "") < str_trunc(cols, 2, ellipsis = ""))

p6H <- covid_corr %>%
  ggplot(aes(x = rows, y = cols, fill = value)) +
  geom_tile() +
  geom_text(aes(label = round(value, 2)), size = 3) +
  scale_fill_continuous_sequential(palette = "Reds3", 
                                   breaks = c(0, 0.5, 1), 
                                   limits = c(0, 1)) +
  labs(x = "", y = "", title = "COVID-19 preprints",
       fill = "Spearman Correlation") +
  theme(plot.title = element_text(size = 10),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill = FALSE) +
  ggsave("outputs/figures/partials/Figure_6H.png", width = 5, height = 4)

p6I <- non_covid_corr %>%
  ggplot(aes(x = rows, y = cols, fill = value)) +
  geom_tile() +
  geom_text(aes(label = round(value, 2)), size = 3) +
  scale_fill_continuous_sequential(palette = "Reds3", 
                                   breaks = c(0, 0.5, 1), 
                                   limits = c(0, 1)) +
  labs(x = "", y = "", title = "non-COVID-19 preprints",
       fill = "Spearman Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  guides(fill = FALSE) +
  ggsave("outputs/figures/partials/Figure_6I.png", width = 5, height = 4)

layout <- "
AAABBB
CCCDDD
EEEFFF
GGHHII
"

# Patchwork. For some reason the widths are not consistent - some manual
# adjustments are made
p <- p6A + p6B + p6C + p6D + p6E + p6F + p6G + p6H + p6I +
  plot_layout(design = layout, 
              widths = c(1,1,0.55,0.55,1,1), 
              guides = 'collect') & theme(legend.position = "top")

p +  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_6.png", width = 10, height = 13)

```

# Supplementary Figure 4: Twitter word clouds, sentiment

```{r}

# Here we generate Supplementary Figure 6. The method is a little messy, as we
# generate wordclouds and venn diagrams which do not seem to play well with 
# ggplot objects. Thus, we create each panel and save it as a separate PNG.
# Panels are then arranged using cowplot's "ggdraw" function. 

# Panel A: Word cloud

library(wordcloud2)
# The below packages are necessary for exporting wordclouds to PNG, see
# https://www.r-graph-gallery.com/196-the-wordcloud2-library.html
library("htmlwidgets")
library("webshot")
#webshot::install_phantomjs()
library(eulerr)
library(cowplot)

# Read tweet data (not shared publicly at the moment)
tweets <- read_csv("data/preprint_tweets_20190101_20201031.csv")

# First we make a wordcloud for tweets of all top-10 tweeted articles
top_preprints <- preprints %>% 
  inner_join(preprint_altmetrics, by = "doi") %>%
  filter(covid_preprint == T) %>% 
  top_n(100, twitter) %>%
  pull(doi)

# Find hashtags associated with top-25 most tweeted covid-19 preprints, only original tweets
hashtags <- tweets %>%
  filter(doi %in% top_preprints) %>% 
  filter(is_retweet == F) %>% 
  select(hashtags) %>%
  na.omit() %>%
  mutate(hashtags = strsplit(as.character(hashtags), ",")) %>% 
  unnest(hashtags) %>%
  mutate(hashtags = str_remove_all(hashtags, "[^a-zA-Z0-9]"), # remove special characters
         hashtags = str_to_lower(hashtags)) %>% # convert all hashtags to lowercase
  filter(!hashtags %in% c("", "covid19", "coronavirus", "ncov2019", "covid", "covid2019", "sarscov2", "2019ncov", "hcov19", "19", "novelcoronavirus", "corona", "coronaovirus", "coronarovirus", "coronarvirus")) %>% # remove empty hashtags and specific keywords
  count(hashtags) %>%
  filter(n > 5) %>%
  mutate(n = log(n+1))

# Generate the wordcloud
p <- wordcloud2(hashtags, color = "random-dark", size = .25)

# Save the wordcloud to a PNG
saveWidget(p, "tweet_wordcloud.html", selfcontained = F)
webshot("tweet_wordcloud.html","outputs/figures/partials/Figure_S4A.png", 
        delay = 30, vwidth = 960, vheight = 720)

# Panel B: Euler diagram of overlap between the top-10 most tweeted, most mentioned
# in news and blogs, most cited and most commented. Here we use the eulerr package

top_twitter <- preprints %>%
  inner_join(preprint_altmetrics, by = "doi") %>%
  filter(covid_preprint == T) %>%
  top_n(10, twitter) %>%
  pull(doi)

top_news <- preprints %>%
  inner_join(preprint_altmetrics, by = "doi") %>%
  filter(covid_preprint == T) %>%
  top_n(10, news) %>%
  pull(doi)

top_blogs <- preprints %>%
  inner_join(preprint_altmetrics, by = "doi") %>%
  filter(covid_preprint == T) %>%
  top_n(10, blogs) %>%
  pull(doi)

top_cited <- preprints %>%
  inner_join(preprint_citations, by = "doi") %>%
  filter(covid_preprint == T) %>%
  top_n(10, citations) %>%
  pull(doi)

top_comments <- preprints %>%
  inner_join(preprint_comments, by = "doi") %>%
  filter(covid_preprint == T) %>%
  top_n(10, comments_count) %>%
  pull(doi)

fit <- list(
  Twitter = top_twitter,
  News = top_news,
  Blogs = top_blogs,
  Citations = top_cited,
  Comments = top_comments)

# Output the Venn diagram to a PNG
png(filename = "outputs/figures/partials/Figure_S4B.png",
    width = 960, height = 720)
plot(euler(fit), 
     fills = qualitative_hcl(n = 5, palette = "Set2"),
     alpha = 0.75,
     labels = list(cex = 1),
     shape = "ellipse",
     quantities = list(cex = 1),
     width = 200)
dev.off()

# Now we fit all the PNG files together. Some manual adjustments may be necessary

pS6A <- ggdraw() +
  draw_image("outputs/figures/partials/Figure_S4A.png") +
  theme(plot.margin = margin(5, 5, 5, 5))

pS6B <- ggdraw() +
  draw_image("outputs/figures/partials/Figure_S4B.png") +
  theme(plot.margin = margin(5, 5, 5, 5))

# Patchwork

layout <- "
A
B
"

pS6A + pS6B +
  plot_layout(design = layout) +
  plot_annotation(tag_levels = "A") +
  ggsave("outputs/figures/Figure_S4.png",
         height = 10, width = 10)

```