exploration.Rmd

---
title: "Some preliminary data exploration"
author: "Anne Hobert"
output: github_document
---

```{r, echo = FALSE, message = FALSE, warning = FALSE}
knitr::opts_chunk$set(
  comment = "#>",
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = TRUE,
  fig.width = 6,
  fig.asp = 0.618,
  out.width = "70%",
  fig.align = "center",
  dpi = 300
)
```

```{r}
library(tidyverse)
library(urltools)
library(cowplot)
```


Here, we perform some preliminary exploration of the generated datasets. We obtain tibbles `pubs_wos`, `upw_evidence`, the matched data frame `pubs_wos_upw`, and the final dataframe `pubs_cat` as described in [data_gathering.Rmd](data_gathering.Rmd).

## How many publications are there? How many can be matched to Unpaywall?

There are `r pubs_cat %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n` different articles in our sample (differentiated by pk_items). `r pubs_cat %>% filter(is.na(DOI)) %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n` of them, or `r round(pubs_cat %>% filter(is.na(DOI)) %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n/pubs_wos %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n, 2)*100 ` %, do not have any DOI information. `r pubs_cat %>% group_by(DOI) %>%   summarise(n = n_distinct(PK_ITEMS)) %>% filter(n==2) %>% count() %>% .$n` DOIs appear twice, that is, they are associated with two different `PK_ITEMS`. The remaining DOIs are linked to a unique `PK_ITEMS`. In total, there are `r pubs_cat %>% filter(!is.na(DOI)) %>% summarise(n=n_distinct(DOI)) %>% .$n` distinct DOIs within the sample. Of these, `r pubs_cat %>% filter(upw_matched) %>% summarise(n = n_distinct(DOI)) %>% .$n` are matched to Unpaywall records, whereas `r pubs_cat %>% filter(!is.na(DOI), upw_matched == FALSE) %>% summarise(n = n_distinct(DOI)) %>% .$n`, or `r round(pubs_cat %>% filter(!is.na(DOI), upw_matched == FALSE) %>% summarise(n = n_distinct(DOI)) %>% .$n/(pubs_cat %>% summarise(n = n_distinct(DOI)) %>% .$n), 2)*100` %, did not correspond to any of the DOIs in our Unpaywall dataset. In other words, `r round(pubs_cat %>% filter(upw_matched == FALSE, is.na(DOI)) %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n/pubs_cat %>% filter(upw_matched == FALSE) %>% summarise(n = n_distinct(PK_ITEMS)) %>% .$n, 2)*100` of the unmatched items do not have any DOI information in WoS.


Briefly evaluate the matching procedure of matching WoS and Unpaywall: The following table shows the number of articles from WoS that have no DOI information (`no_doi`), articles that have DOI information in WoS but could not be matched to any record in the Unpwaywall dataset we used (`unmatched`), articles matched to Unpaywall for which Unpaywall did not find any open version (`upw_closed`) and articles with OA evidence in Unpaywall (`OA`).
```{r}
pubs_wos_upw %>% 
  mutate(matched_oa = case_when(
    is.na(DOI) ~ 'no_doi',
    is.na(upw_doi) ~ 'unmatched',
    !is_oa ~ 'upw_closed',
    TRUE ~ 'OA'
  )) %>% 
  group_by(matched_oa) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  # filter(!matched_oa %in% c('unmatched', 'no_doi')) %>% 
  mutate(prop = n / sum(n))%>% 
  knitr::kable()
```

Extract unmatched articles with DOI for checking against API:

```{r}
pubs_wos_upw %>% 
  filter(!is.na(DOI), is.na(upw_doi)) %>% 
#  summarise(n = n_distinct(DOI))
  select(DOI) %>% 
  distinct() %>% 
  readr::write_csv("data/unmatched_wos_dois.csv")
```


## Fulltext links occuring for multiple dois

How many locations in total? 

```{r}
upw_evidence %>% 
  filter(is_oa == TRUE) %>% 
  group_by(url_for_landing_page, url_for_pdf) %>% 
  summarise(n = n_distinct(upw_doi)) %>% 
  arrange(desc(n))

```

## Number of Full-Text links per OA Article 

How many OA DOIs have more than one OA location?

```{r}
upw_evidence %>% 
  filter(is_oa == TRUE) %>% 
  group_by(upw_doi) %>% 
  summarise(n_oa_loc = n()) %>% 
  group_by(n_oa_loc > 1) %>% 
  count() %>% 
  mutate(prop = n / sum(.$n))
```

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
excl_from_analysis <- readr::read_csv("data/exclude_from_analysis.csv")
pubs_cat_details <- pubs_cat_details %>% 
  anti_join(excl_from_analysis, by = c("INST_NAME" = "Universitäten"))
rm(excl_from_analysis)

pubs_cat_details %>% 
  filter(!is.na(upw_doi)) %>% 
  group_by(upw_doi) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  filter(n == 2) %>% 
  count()

  
pubs_cat_details %>% 
  filter(is_oa == TRUE) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  select(n_total, PK_ITEMS, url_for_landing_page, url_for_pdf, evidence) %>% 
  distinct() %>% 
  group_by(PK_ITEMS, n_total) %>% 
  summarise(n_loc = n()) %>%  
  group_by(n_loc > 1, n_total) %>% 
  summarise(n_items = n()) %>% 
  mutate(prop = n_items/n_total)

pubs_cat_details %>% 
  filter(oa_category =="not_oa") %>% 
  group_by(PK_ITEMS) %>% 
  mutate(n_loc = n()) %>% 
  filter(n_loc > 1)
```


How many OA DOIs have more than one OA category according to our schema?

```{r}
pubs_cat %>% 
  filter(oa_category != "not_oa") %>% 
  group_by(DOI) %>% 
  summarise(n_oa_cat = n_distinct(oa_category)) %>% 
  group_by(n_oa_cat > 1) %>% 
  count() %>% 
  mutate(prop = n / sum(.$n))
```


```{r}
oa_fxt_count <- upw_evidence %>% 
  filter(!is.na(host_type)) %>%
  group_by(upw_doi) %>%
  mutate(n_oa_versions = n()) %>%
  ungroup() %>%
  group_by(DOI, n_oa_versions, host_type) %>%
  summarise(ftxt_n = n()) %>%
  ungroup()
```

```{r}
oa_fxt_count %>%
  distinct(DOI, n_oa_versions) %>%
  ggplot(aes(x = n_oa_versions)) +
  geom_histogram(binwidth = 1, fill="#56B4E9", alpha=0.5, color = "#56B4E9") +
  geom_vline(aes(xintercept = mean(n_oa_versions, na.rm = T)),
             colour = "#E69F00", linetype ="dashed", size = .8) +
  geom_vline(aes(xintercept = median(n_oa_versions, na.rm = T)),
             colour = "red", linetype ="dashed", size = .8) +
  theme_minimal_hgrid(12) +
  labs(x = "Number of Full-Text Links",
       y = "German OA Articles")
```

### per host type

```{r}
# compute densities for full text numbers lengths
oa_fxt_count %>%
  #filter(!is.na(host_type)) %>% 
  #distinct() %>% 
  ggplot(aes(x = ftxt_n, fill = host_type, color = host_type)) +
  geom_histogram(binwidth = 1, alpha=0.5) +
  facet_grid(~ host_type, scales = "free_x") +
  theme_minimal_hgrid(12) +
  labs(x = "Number of Full-Text Links",
       y = "German OA Articles") 
```

## Number of OA Articles per Host Type

```{r}
oa_fxt_count %>% 
  group_by(DOI) %>%
  mutate(n_host_type = n()) %>%
  mutate(host_type = ifelse(n_host_type == 2, "publisher & repository", host_type)) %>%
  ungroup() %>%
  group_by(host_type) %>%
  summarise(n = n_distinct(DOI)) %>%
  mutate(prop = n / sum(n)) %>%
  knitr::kable()
```

### Repository analysis using OpenDOAR

```{r}
excl_from_analysis <- readr::read_csv("data/exclude_from_analysis.csv")
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
 repo_df <- pubs_cat_details %>% 
  filter(host_type == "repository") %>%
   ## exclude some inst
  anti_join(excl_from_analysis, by = c("INST_NAME" = "Universitäten"))
```

```{r}
opendoar <- readr::read_csv("data/opendoar_data_tidier.csv") %>%
  mutate(repo_domain = domain(repo_url)) %>%
  mutate(repo_domain = gsub("www.", "", repo_domain))
```


#### Number of self-archived full-texts by repository type

```{r}
repo_df %>% 
  distinct(DOI, pdf_domain, landing_domain) %>% 
  left_join(opendoar,  by = c("pdf_domain" = "repo_domain", "landing_domain" = "repo_domain")) %>% 
  ungroup() %>% 
  group_by(repo_type) %>% 
  summarise(n = n()) %>% 
  arrange(desc(n)) %>% 
  mutate(prop = n / sum(n)) %>%
  knitr::kable()
```

without "pdfs.semanticscholar.org".

```{r}
repo_df %>% 
  distinct(DOI, pdf_domain, landing_domain) %>% 
  filter(pdf_domain != "pdfs.semanticscholar.org") %>%
  left_join(opendoar,  by = c("pdf_domain" = "repo_domain", "landing_domain" = "repo_domain")) %>% 
  ungroup() %>% 
  group_by(repo_type) %>% 
  summarise(n = n()) %>% 
  arrange(desc(n)) %>% 
  mutate(prop = n / sum(n)) %>%
  knitr::kable()
```

### Top Repository Domains

```{r}
repo_df %>%
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(domains) %>%
  summarise(n = n_distinct(DOI)) %>%
  mutate(prop = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
#  head(20) %>%
#  knitr::kable()
  readr::write_csv("data/export_tables/repo_domain_counts.csv")
```

per category:

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
pubs_cat_details %>%
    anti_join(excl_from_analysis, by = c("INST_NAME" = "Universitäten")) %>% 
  filter(!oa_category %in% c("not_oa", "full_oa_journal", "other_oa_journal")) %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(oa_category, domains) %>%
  summarise(n = n()) %>%
  ungroup() %>% 
  group_by(oa_category) %>% 
#  mutate(prop = round(n / sum(n) * 100)) %>%
  arrange(oa_category, desc(n)) %>%
  pivot_wider(names_from = oa_category, values_from = n, values_fill = list(n = 0)) %>% 
#  head(20) %>%
#  knitr::kable()
  readr::write_csv("data/export_tables/repo_domain_category_counts.csv")
```

Development of disciplinary repos:

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
top_disc_repos <- pubs_cat_details %>%
  filter(oa_category == "opendoar_subject") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(domains) %>%
  summarise(n = n_distinct(DOI)) %>%
  arrange(desc(n)) %>% 
  head(20) %>% 
  select(domains)
pubs_cat_details %>%
  filter(oa_category == "opendoar_subject") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  inner_join(top_disc_repos) %>% 
  group_by(PUBYEAR, domains) %>%
  summarise(n = n_distinct(DOI)) %>%
  ggplot(aes(x = PUBYEAR, y = n)) +
  geom_line() + 
  facet_wrap(~ domains)
```


```{r}
repo_df %>% 
  distinct(pdf_domain, landing_domain) %>% 
  left_join(opendoar,  by = c("pdf_domain" = "repo_domain", "landing_domain" = "repo_domain"))

repo_df %>% 
  distinct(FK_ITEMS, DOI, pdf_domain, landing_domain) %>% 
  left_join(opendoar,  by = c("pdf_domain" = "repo_domain", "landing_domain" = "repo_domain")) %>% 
  ungroup() %>% 
  write_csv("data/german_unis_url_matching.csv")
```
## Special role of MPG

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
pubs_cat_details %>% 
  filter(sector == "Max-Planck-Gesellschaft",
#         host_type == "repository") %>% 
         oa_category == "opendoar_subject") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(domains) %>%
  summarise(n = n()) %>%
  mutate(prop = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
 head(20) %>%
 knitr::kable()
```

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
pubs_cat_details %>% 
  filter(sector == "Max-Planck-Gesellschaft",
#         host_type == "repository") %>% 
         oa_category == "other_oa_journal") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(pdf_domain) %>%
  summarise(n = n()) %>%
  mutate(prop = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
 head(20) %>%
 knitr::kable()
```

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
pubs_cat_details %>% 
  filter(sector == "Max-Planck-Gesellschaft",
#         host_type == "repository") %>% 
         oa_category == "opendoar_other") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(pdf_domain) %>%
  summarise(n = n()) %>%
  mutate(prop = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
 head(20) %>%
 knitr::kable()
```

```{r}
library(bigrquery)
con_bq <- dbConnect(
  bigrquery::bigquery(),
  project = "api-project-764811344545",
  dataset = "oadoi_full"
)
pubs_cat <- readr::read_csv("data/pubs_cat.csv")
mpg_dois <- pubs_cat %>% 
  filter(oa_category == "other_oa_journal") %>% 
  select(DOI) %>% 
  distinct()
schema <- data_frame("DOI" = character())
bq_dois <- bigrquery::bq_table(project = "api-project-764811344545",
                            dataset = "oadoi_full",
                            table = "mpg_dois")
## publishers
bq_mpg_dois <- bq_table_create(
  bq_dois,
  fields = schema
)
bq_table_upload(bq_mpg_dois, mpg_dois)
cr_publisher_mpg_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT publisher, c.doi FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.mpg_dois` u on c.doi = u.DOI", destination_table = "api-project-764811344545.oadoi_full.cr_publisher_mpg")
cr_top_publisher_mpg_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT publisher, count(distinct c.doi) as n FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.mpg_dois` u on c.doi = u.DOI group by publisher order by n desc LIMIT 20", destination_table = "api-project-764811344545.oadoi_full.cr_top_publisher_mpg")
cr_dois_top_publisher_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT c.publisher, doi FROM `oadoi_full.cr_publisher_mpg` c inner join `oadoi_full.cr_top_publisher_mpg` t on c.publisher = t.publisher")
cr_publisher_mpg <- bq_table_download(cr_dois_top_publisher_bq)
bq_table_delete(cr_top_publisher_mpg_bq)
bq_table_delete(cr_publisher_mpg_bq)

pubs_cat %>% 
  inner_join(cr_publisher_mpg, by = c("DOI" = "doi")) %>% 
  group_by(publisher) %>% 
  summarise(n = n_distinct(DOI)) %>%
  arrange(desc(n)) %>% 
  readr::write_csv("data/export_tables/publishers_mpg_other_oa_journal.csv")

## journals
cr_journal_mpg_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT container_title, publisher, c.doi FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.mpg_dois` u on c.doi = u.DOI", destination_table = "api-project-764811344545.oadoi_full.cr_journal_mpg")
cr_top_journal_mpg_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT container_title, publisher, count(distinct c.doi) as n FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.mpg_dois` u on c.doi = u.DOI group by container_title, publisher order by n desc LIMIT 30", destination_table = "api-project-764811344545.oadoi_full.cr_top_journal_mpg")
cr_dois_top_journal_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT c.container_title, c.publisher, doi FROM `oadoi_full.cr_journal_mpg` c inner join `oadoi_full.cr_top_journal_mpg` t on c.container_title = t.container_title and c.publisher = t.publisher")
cr_journal_mpg <- bq_table_download(cr_dois_top_journal_bq)
bq_table_delete(cr_top_journal_mpg_bq)
bq_table_delete(cr_journal_mpg_bq)

pubs_cat %>% 
  inner_join(cr_journal_mpg, by = c("DOI" = "doi")) %>% 
  group_by(container_title, publisher) %>% 
  summarise(n = n_distinct(DOI)) %>%
  arrange(desc(n)) %>% 
  readr::write_csv("data/export_tables/journals_mpg_other_oa_journal.csv")
```

## Drop in other repo

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")

most_prev_doms <- pubs_cat_details %>% 
  filter(oa_category == "other_repo") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  group_by(domains) %>%
  summarise(n = n()) %>%
  # mutate(prop = n / sum(n) * 100) %>%
  arrange(desc(n)) %>% 
  head(25) %>% 
  select(domains)

pubs_cat_details %>% 
  filter(oa_category == "other_repo") %>% 
  mutate(domains = case_when(
    is.na(pdf_domain) && is.na(landing_domain) ~ NA_character_,
    is.na(pdf_domain) ~ landing_domain,
    TRUE ~ pdf_domain
  )) %>% 
  inner_join(most_prev_doms) %>% 
  group_by(domains, PUBYEAR) %>%
  summarise(n = n()) %>% 
  ggplot(aes(x = PUBYEAR, y = n)) +
  geom_line()+
  facet_wrap(~ domains)
  
```
## rise in full OA journal


```{r}
library(bigrquery)
con_bq <- dbConnect(
  bigrquery::bigquery(),
  project = "api-project-764811344545",
  dataset = "oadoi_full"
)
pubs_cat <- readr::read_csv("data/pubs_cat.csv")
excl_from_analysis <- readr::read_csv("data/exclude_from_analysis.csv")
pubs_cat <- pubs_cat %>% 
  anti_join(excl_from_analysis, by = c("INST_NAME" = "Universitäten"))
rm(excl_from_analysis)
full_dois <- pubs_cat %>% 
  filter(oa_category == "full_oa_journal") %>% 
  select(DOI) %>% 
  distinct()
schema <- data_frame("DOI" = character())
bq_dois <- bigrquery::bq_table(project = "api-project-764811344545",
                            dataset = "oadoi_full",
                            table = "full_dois")
bq_full_dois <- bq_table_create(
  bq_dois,
  fields = schema
)
bq_table_upload(bq_full_dois, full_dois)
cr_publisher_full_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT publisher, c.doi FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.full_dois` u on c.doi = u.DOI", destination_table = "api-project-764811344545.oadoi_full.cr_publisher_full")
cr_top_publisher_full_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT publisher, count(distinct c.doi) as n FROM `cr_dump_march_20.cr_08_20` c inner join `oadoi_full.full_dois` u on c.doi = u.DOI group by publisher order by n desc LIMIT 20", destination_table = "api-project-764811344545.oadoi_full.cr_top_publisher_full")
cr_dois_top_publisher_bq <- bigrquery::bq_project_query("api-project-764811344545", query = "SELECT c.publisher, doi FROM `oadoi_full.cr_publisher_full` c inner join `oadoi_full.cr_top_publisher_full` t on c.publisher = t.publisher")
cr_publisher_full <- bq_table_download(cr_dois_top_publisher_bq)
bq_table_delete(cr_top_publisher_full_bq)
bq_table_delete(cr_publisher_full_bq)
bq_table_delete(bq_full_dois)

pubs_cat %>% 
  inner_join(cr_publisher_full, by = c("DOI" = "doi")) %>% 
  group_by(PUBYEAR, publisher) %>% 
  summarise(n = n_distinct(DOI)) %>% 
  ggplot(aes(x=PUBYEAR, y = n)) +
  geom_line() +
  facet_wrap(~ publisher)

pubs_cat %>% 
  inner_join(cr_publisher_full, by = c("DOI" = "doi")) %>% 
  group_by(publisher) %>% 
  summarise(n = n_distinct(DOI)) %>% 
  arrange(desc(n))

pubs_cat %>% 
  inner_join(cr_publisher_full, by = c("DOI" = "doi")) %>% 
  filter(publisher == "Public Library of Science (PLoS)") %>% 
  group_by(PUBYEAR) %>% 
  summarise(n = n_distinct(DOI))

```
## Comparison Gold-list, Unpaywall

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
excl_from_analysis <- readr::read_csv("data/exclude_from_analysis.csv")
pubs_cat_details <- pubs_cat_details %>% 
  anti_join(excl_from_analysis, by = c("INST_NAME" = "Universitäten"))
rm(excl_from_analysis)

# pubs in both
pubs_cat_details %>%
  filter(host_type == "publisher") %>% 
  filter(issn_gold ==TRUE, journal_is_oa == TRUE) %>% 
  select(`ISSN-L`) %>% 
  distinct() -> both

pubs_cat_details %>%
  anti_join(both) %>%
  filter(issn_gold == TRUE & (host_type == "publisher"| is.na(host_type))) %>%
  select(`ISSN-L`) %>% 
  distinct() -> only_list
pubs_cat_details %>%
  anti_join(both) %>%
  filter(journal_is_oa == TRUE, host_type == "publisher") %>%
  select(`ISSN-L`) %>% 
  distinct() -> only_upw
```


## Institutional details

## Subject Uni wuppertal

```{r}
pubs_cat_details <- readr::read_csv("data/pubs_cat_details.csv")
pubs_cat_details %>% 
  filter(oa_category != "not_oa", INST_NAME == "Bergische Universität Wuppertal") %>% 
  mutate(oa_host = fct_collapse(
    oa_category,
    journal = c("full_oa_journal", "other_oa_journal"),
    repo = c(
      "opendoar_inst",
      "opendoar_subject",
      "opendoar_other",
      "other_repo"
    )
  )) %>%
  group_by(oa_host, pdf_domain) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  arrange(desc(n)) %>% View()
```


### Outliers with all OA publications

Investigate, which institutions are the outliers with only open access articles.

```{r}
outliers <- pubs_cat %>%
  filter(INST_NAME %in% c("Heidelberg Graduate School of Fundamental Physics", "Schloss Dagstuhl - Leibniz-Zentrum für Informatik GmbH (LZI)"))
outliers %>%
  group_by(INST_NAME) %>% 
  summarise(n = n_distinct(PK_ITEMS))
outliers %>% 
  group_by(INST_NAME, oa_category) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  arrange(INST_NAME, desc(n))
outliers %>%
  group_by(INST_NAME, url_domain, oa_category) %>%
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  arrange(desc(n))
```
Schloss Dagstuhl with just one publication, that is OA. 2nd one is Heidelberg Graduate School of Fundamental Physics with 56 publications, 55 of which are published in the journal Astronomy and Astrophysics. Most of these publications are tagged as journal based OA but are also deposited in various repositories.

### Highest / lowest OA shares

Which institutions have the highest OA shares (apart from the two just discussed ones)?
```{r}
oa_shares %>% 
  select(INST_NAME, oa_share) %>%
  top_n(15) %>% 
  arrange(desc(oa_share))
```

Mostly (astro-)physics-related MP institutes, but also Leibniz-Institut für Astrophysik, some computer science / biology related institutes.

Lowest shares?

```{r}
oa_shares %>% 
  select(INST_NAME, n_total, oa_share) %>% 
  filter(oa_share == 0)
```
OA share zero mostly for institutions with few publications in data set (less then 10 in total), exceptions: "Staatliche Hochschule für Musik Freiburg im Breisgau" (39 pubs), "Herder-Institut für historische Ostmitteleuropaforschung"(14 pubs). 

Distribution of OA shares?

```{r}
oa_shares %>% 
  ggplot(aes(x = oa_share)) +
    geom_histogram(binwidth = 1) +
    geom_vline(aes(xintercept = 7), col = "red")
```

### Total publication output

Publication output of sectors/institutions?
```{r}
pubs_cat %>% 
  group_by(sector) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  ggplot(aes(x = sector, y = n)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
  geom_col()
pubs_cat %>% 
  group_by(INST_NAME) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  filter(n < 10000) %>% 
  ggplot(aes(x =n)) +
  geom_histogram(bins = 150) +
  # geom_vline(xintercept = 35000, color = "red")+
  # geom_vline(xintercept = 21000, color = "red")+
  # geom_vline(xintercept = 10000, color = "red")+
   geom_vline(xintercept = 6000, color = "red")
```

## Validation

Extract sample DOIs of the matched records in order to manually check the OA categorisation:

```{r}
sample_rows <- pubs_cat_update %>% 
  filter(upw_matched == TRUE) %>% 
  select(upw_doi) %>% 
  distinct() %>% 
  sample_n(100)
sample_rows <- sample_rows %>% 
  select(upw_doi) %>% 
  distinct() %>% 
  left_join(pubs_cat_update, by = "upw_doi")
sample_rows %>% group_by(oa_category, oa_category_old) %>% summarise(n = n_distinct(upw_doi)) %>% arrange(desc(n))
sample_rows %>% 
  select(upw_doi) %>% 
  distinct() %>% 
  readr::write_csv("data/sample_dois_manual.csv")
sample_rows %>% 
  select(upw_doi, oa_category, evidence, host_type, repository_institution, url_for_landing_page, url_for_pdf) %>% 
  distinct() %>% 
  readr::write_csv("data/sample_manual_check_2.csv")
```

Explore registration agencies of the unmatched DOIs:
```{r}
unmatched_dois <- readr::read_csv("data/unmatched_dois_agency.csv")
unmatched_dois %>% 
  ## update after exclusion of theological etc. universities
  inner_join(distinct(select(pubs_cat, DOI)), by = c("doi" = "DOI")) %>% 
  group_by(doi_agency_label) %>% 
  summarise(n = n_distinct(doi)) %>% 
  arrange(desc(n)) %>% 
  mutate(doi_agency_label = replace_na(doi_agency_label, "No agency found")) %>% 
  readr::write_csv("data/unmatched_dois_agency_counts.csv")
```


## Comparison of matching strategies

After the first validation we updated the matching strategy regarding OpenDOAR information. Here, we want to briefly investigate the effects of this change in procedure.

```{r}
pubs_cat_update %>% 
  mutate(new_cat = ifelse(oa_category == oa_category_old, FALSE, TRUE)) %>% 
  group_by(new_cat) %>% 
  summarise(n = n_distinct(PK_ITEMS))
pubs_cat_update %>% 
  mutate(new_cat = ifelse(oa_category == oa_category_old, FALSE, TRUE)) %>%
  filter(new_cat == TRUE) %>% 
  group_by(oa_category, oa_category_old) %>% 
  summarise(n = n_distinct(PK_ITEMS))
pubs_cat_update %>% 
  mutate(new_cat = ifelse(oa_category == oa_category_old, FALSE, TRUE)) %>%
  filter(new_cat == TRUE) %>% 
  group_by(repository_institution, pdf_domain, landing_domain, oa_category, oa_category_old) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  arrange(desc(n)) %>% View()
pubs_cat_update %>% 
  mutate(new_cat = ifelse(oa_category == oa_category_old, FALSE, TRUE)) %>%
  filter(new_cat == TRUE) %>% 
  group_by(pdf_domain == landing_domain) %>% 
  summarise(n = n_distinct(PK_ITEMS)) %>% 
  arrange(desc(n))
```


## Exemplary analytics on national, sectoral and institutional level

### National level

First, we look at the highest aggregation level, namely the national level. Here, publications of all institutions of the selected sectors (5+2) are counted per year and category. 

#### by host type

The first figure shows the development over time of the total number of publications per host type, that is, we collate all articles where access is provided via a journal based platform as `OA (journal)`, all repository based OA is collated as `OA (repository)`, and all remaining articles are labelled as `Not OA`.

```{r}
pubs_cat %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo")
  )) %>% 
  group_by(PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    #geom_line(aes(color = oa_category)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    # geom_bar(aes(fill = oa_category), stat = "identity", position = "dodge") +
    #scale_color_viridis(end = 0.85, discrete = TRUE, direction = -1) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA Categories", labels = c("OA (journal)", "OA (repository)", "Not OA")) +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA host type (national level)") +
    theme(legend.position = "bottom")
```

In comparison to all articles published:

```{r}
all_articles <- pubs_cat%>% 
  group_by(PUBYEAR) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS))

pubs_cat %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  ))  %>% 
  group_by(PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
#  .$oa_category %>% levels()
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_col(data = all_articles, aes(color = "All Articles"), fill = "#b3b3b3", size = 0.05) +
    scale_color_manual(values = "#ffffff", name = "", labels = "All Articles") +
    geom_col(aes(fill = oa_category)) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "", labels = c("OA (journal)", "OA (repository)", "Not OA"))+
    facet_wrap(~ oa_category, nrow = 1) +
    theme(strip.text.x = element_blank()) +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA host type (national level)") +
    theme(legend.position = "bottom")
```


Looking at proportions:

```{r}
pubs_cat %>% 
  group_by(PUBYEAR) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo")
  )) %>% 
  group_by(PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    #geom_line(aes(color = oa_category)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    # geom_bar(aes(fill = oa_category), stat = "identity", position = "dodge") +
    #scale_color_viridis(end = 0.85, discrete = TRUE, direction = -1) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, labels = c("OA (journal)", "OA (repository)", "Not OA"), direction = -1) +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA host type", fill = "OA Category") +
    theme(legend.position = "right") 
    # theme(axis.title = element_text(size = 16),
    #     plot.title = element_text(size = 24),
    #     axis.text = element_text(size = 14),
    #     legend.text = element_text(size = 14),
    #     legend.title = element_text(size = 16))
```

As faceted graph:

```{r}
pubs_cat %>% 
  group_by(PUBYEAR) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo")
  )) %>% 
  group_by(PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    #geom_line(aes(color = oa_category)) +
    geom_col(aes(fill = oa_category)) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, labels = c("OA (journal)", "OA (repository)", "Not OA"), direction = -1) +
    facet_wrap(~ oa_category, nrow = 1) +
    theme(strip.text.x = element_blank()) +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA host type", fill = "OA Category") +
    theme(legend.position = "bottom")
```

#### by OA category

The following figure shows very similar data as the one before - only that here, all OA categories from the schema are shown individually, that is, journal based and repository based OA are further distinguished into subgroups.

```{r}
pubs_cat %>% 
  group_by(PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_line(aes(color = oa_category), size = 1) +
    #geom_col(aes(fill = oa_category), position = "dodge") +
    # geom_bar(aes(fill = oa_category), stat = "identity", position = "dodge") +
    #scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1) +
    scale_color_viridis(end = 0.95, discrete = TRUE, direction = -1, name= "OA category") +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA category (national level)") +
    theme(legend.position = "right")
```

In comparison to all articles published:

```{r}
all_articles <- pubs_cat%>% 
  group_by(PUBYEAR) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS))

pubs_cat %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  group_by(PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
#  .$oa_category %>% levels()
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_col(data = all_articles, aes(color = "All Articles"), fill = "#b3b3b3", size = 0.05) +
    scale_color_manual(values = "#ffffff", name = "", labels = "All Articles") +
    geom_col(aes(fill = oa_category)) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA categories")+
    facet_wrap(~ oa_category, ncol = 2) +
    theme(strip.text.x = element_blank()) +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA host type (national level)") +
    theme(legend.position = "right")
```

Looking at proportions, we get:

```{r}
pubs_cat %>% 
  group_by(PUBYEAR) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  group_by(PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    #geom_line(aes(color = oa_category), size = 1) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category") +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA category (national level)") +
    theme(legend.position = "bottom")
```

As faceted graph:

```{r}
pubs_cat %>% 
  group_by(PUBYEAR) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  group_by(PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    #geom_line(aes(color = oa_category)) +
    geom_col(aes(fill = oa_category)) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category") +
    facet_wrap(~ oa_category, ncol = 2) +
    theme(strip.text.x = element_blank()) +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA host type", fill = "OA Category") +
    theme(legend.position = "right")
```

### Level of sectors

Next, we go one level down and look at how the separate sectors developed.

Number of OA articles per sector in comparison to all articles published within the sector:

```{r, fig.asp=1}
pubs_cat %>% 
  mutate(oa_category = fct_relevel(factor(oa_category), "full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa")) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    is_oa = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  ))  %>%
  group_by(sector, PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>% 
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_area(aes(fill = fct_rev(oa_category), group = fct_rev(oa_category)),  alpha = 0.8) +
    scale_fill_manual(values = c("#b3b3b3a0", "#56B4E9"), name = "Is OA?", labels = c("Not OA", "OA"))+
    facet_wrap(~ sector, scales = "free", ncol = 2) +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of OA Articles (level of sectors)") +
    theme(legend.position = "right")
```

#### by host type

The first figure depicts the number of publications per host type category over time for each sector.

```{r, fig.asp=1}
#, fig.asp=1.5}
pubs_cat %>% 
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  )) %>% 
  group_by(sector, PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    facet_wrap(~ sector, scales = "free", ncol = 2) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category", labels = c("OA (journal)", "OA (repository)", "Not OA")) +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA host type (level of sectors)", color = "oa_category") +
    theme(legend.position = "right")
```

Looking at proportions:
```{r, fig.asp=1}
#, fig.asp=1.5}
pubs_cat %>%
  group_by(PUBYEAR, sector) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  )) %>% 
  group_by(sector, PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total *100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    facet_wrap(~ sector, scales = "free") +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category", labels = c("OA (journal)", "OA (repository)", "Not OA")) +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA host type (level of sectors)", color = "oa_category") +
    theme(legend.position = "bottom")
```

Facet grid:

```{r, fig.asp=1}
pubs_cat_abbr <- pubs_cat %>% 
  mutate(sector = case_when(
    sector == "Hochschulen" ~ "HS",
    sector == "Leibniz-Gemeinschaft" ~ "LG",
    sector == "Helmholtz-Gemeinschaft" ~ "HG",
    sector == "Max-Planck-Gesellschaft" ~ "MPG",
    sector == "Ressortforschung-Bund" ~ "R-B",
    sector == "Ressortforschung-Länder" ~ "R-L",
    sector == "Fraunhofer-Gesellschaft" ~ "FG"
  ))
pubs_cat_abbr %>%
  group_by(PUBYEAR, sector) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>%
    mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  )) %>% 
  group_by(sector, PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    geom_col(aes(fill = oa_category)) +
    facet_grid(sector ~ oa_category) +
    theme(strip.text.x = element_blank()) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category", labels = c("OA (journal)", "OA (repository)", "Not OA")) +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA host type (level of sectors)", color = "oa_category") +
    theme(legend.position = "bottom")
```

#### by OA category

Distinguishing all the single OA categories:

```{r, fig.asp=1}
pubs_cat %>%
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  group_by(sector, PUBYEAR, oa_category) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  ggplot(aes(x = PUBYEAR, y = number_of_articles)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    facet_wrap(~ sector, scales = "free", ncol = 2) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category") +
    labs(x = "Publication Year", y = "Number of Articles", title = "Number of Articles per OA category (level of sectors)", color = "oa_category") +
    theme(legend.position = "bottom")
```

Looking at proportions:

```{r, fig.asp=1}
pubs_cat %>%
  group_by(PUBYEAR, sector) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  group_by(sector, PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    geom_col(aes(fill = oa_category), position = "dodge") +
    facet_wrap(~ sector, scales = "free", ncol = 2) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category") +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA category (level of sectors)", color = "oa_category") +
    theme(legend.position = "bottom")
```


Facet grid:

```{r, fig.asp=1.2}
pubs_cat_abbr <- pubs_cat %>% 
  mutate(sector = case_when(
    sector == "Hochschulen" ~ "HS",
    sector == "Leibniz-Gemeinschaft" ~ "LG",
    sector == "Helmholtz-Gemeinschaft" ~ "HG",
    sector == "Max-Planck-Gesellschaft" ~ "MPG",
    sector == "Ressortforschung-Bund" ~ "R-B",
    sector == "Ressortforschung-Länder" ~ "R-L",
    sector == "Fraunhofer-Gesellschaft" ~ "FG"
  ))
pubs_cat_abbr %>%
  group_by(PUBYEAR, sector) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>% 
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>%
  group_by(sector, PUBYEAR, oa_category, n_total) %>% 
  summarise(number_of_articles = n_distinct(PK_ITEMS)) %>%
  mutate(prop = number_of_articles / n_total * 100) %>% 
  ggplot(aes(x = PUBYEAR, y = prop)) +
    geom_col(aes(fill = oa_category)) +
    facet_grid(sector ~ oa_category) +
    theme(strip.text.x = element_blank()) +
    scale_fill_viridis(end = 0.95, discrete = TRUE, direction = -1, name = "OA category") +
    labs(x = "Publication Year", y = "Proportion of Articles (in %)", title = "Proportion of Articles per OA category (level of sectors)", color = "oa_category") +
    theme(legend.position = "bottom")
```

### Institutional level

Lastly, we go down to the institutional level. Since there are more than 400 institutions in the dataset, we do not want to present figures for individual institutions. Instead, we want to show the variety of open access shares and strategies. To this end, we first calculate the individual OA shares.

```{r}
oa_shares <- pubs_cat %>%
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    is_oa = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  )) %>% 
  group_by(INST_NAME) %>% 
  mutate(n_total = n_distinct(PK_ITEMS)) %>% 
  ungroup() %>%
  group_by(INST_NAME, oa_category, n_total) %>% 
  summarise(n_cat = n_distinct(PK_ITEMS)) %>% 
  pivot_wider(names_from = oa_category,
              values_from = n_cat,
              values_fill = list(n_cat = 0)) %>% 
  mutate(oa_share = is_oa/n_total *100) %>%
  rename(n_oa = is_oa, n_not_oa = not_oa) %>% 
  ungroup()
```

We now display the variability among institutions per sector with respect to the total publication output.

```{r, fig.asp=1}
pubs_cat %>%
  mutate(oa_category = factor(oa_category, levels = c("full_oa_journal", "other_oa_journal", "opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo", "not_oa"))) %>% 
  mutate(oa_category = fct_collapse(oa_category,
    journal_oa = c("full_oa_journal", "other_oa_journal"),
    repo_oa = c("opendoar_inst", "opendoar_subject", "opendoar_other", "other_repo"),
    not_oa = "not_oa"
  )) %>% 
  left_join(oa_shares, by = 'INST_NAME') %>% 
  group_by(sector, INST_NAME, n_total, oa_share) %>%
  summarise() %>% 
  ggplot(aes(x = n_total, y = oa_share)) +
    geom_point(aes(color = sector)) + 
    scale_x_log10() + 
    scale_color_viridis(end = 0.95, direction = -1, option = "plasma", discrete = TRUE, name = "Sectors")+
    facet_wrap(~ sector, ncol = 2) +
    labs(x = "Total publication output (logarithmic scale)", y = "OA share (in %)", title = "Open Access shares of research institutions in Germany") +
    theme(legend.position = "right")
```

Sectorwise comparison of OA shares

```{r, fig.asp=0.8}
pubs_cat %>%
  left_join(oa_shares, by = 'INST_NAME') %>%
  group_by(INST_NAME, sector, oa_share, n_total) %>% 
  summarise() %>% 
  group_by(sector) %>% 
  summarise_at(vars(n_total, oa_share), funs(min, max, mean, median,sd)) %>% 
  select(sector, starts_with('oa_share'), starts_with('n_total'))

n_fun <- function(x){
  return(data.frame(y = 105,
                    label = paste ("n=", length(x))))
}

pubs_cat %>%
  group_by(sector) %>% 
  mutate(n_sector = n_distinct(PK_ITEMS),
         sector_cat = case_when(
           n_sector > 500000 ~ 'high_n_pubs',
           n_sector > 50000 ~ 'middle_n_pubs',
           TRUE ~ 'low_n_pubs'
         )) %>% 
  ungroup() %>% 
  left_join(oa_shares, by = 'INST_NAME') %>%
  group_by(INST_NAME, sector, oa_share, sector_cat) %>% 
  summarise() %>% 
  ggplot(aes(x = sector, y = oa_share)) + 
    geom_boxplot(aes(color = sector_cat), varwidth = TRUE, size = 1) +
#    geom_boxplot() +
 #   stat_summary(fun.data = n_fun, geom = "text", hjust = 0.5) +
    geom_jitter(alpha = 0.1) +
    scale_color_viridis(end = 0.8, labels = c("high", "medium", "low"), discrete = TRUE, direction = -1, option = "plasma")+
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(title = "OA shares of institutions per sector", x = "", y = "OA share (in %)", color = "Publication output",
         caption = "Widths of the boxes corresponds to the number of institutions per sector,\n gray points display oa shares for individual institutions.") 
      # theme(axis.title = element_text(size = 24),
      #   plot.title = element_text(size = 36),
      #   axis.text = element_text(size = 21),
      #   legend.text = element_text(size = 21),
      #   legend.title = element_text(size = 24))
```