---
title: "Coral Health Synthesis"
author: "Samantha Burke"
date: '2022-10-11'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Load in packages
```{r}
library(tidyverse)
library(readxl)
library(visdat)
library(maps)
library(here)
library(ggtreeExtra)
library(wordcloud)
library(wordcloud2)
library(RColorBrewer)
library(tm)
library(rAltmetric)
library(dimensionsR)
pacman::p_load(plyr,
               tibble,
               dplyr,
               stringr,
               knitr,
               forcats,
               ggplot2,
               ggraph,
               hrbrthemes,
               bibliometrix,
               patchwork,
               migest,
               circlize,
               rotl,
               ape,
               ggnewscale,
               ggtree,
               igraph,
               stringi,
               tidystringdist,
               stringdist,
               fuzzyjoin,
               ggthemes,
               viridis,
               cowplot,
               mapproj)
```


# Load in data
```{r}
# extracted data load
#systmap <- read.csv(here("data", "Systematic Map Data.csv"))
#systmap <- readRDS(here::here("data", "Systematic Map Data.rds"))
#systmap <- read.csv(here::here("data", "Final Systematic Map Data.csv"))
systmap <- readRDS(here::here("Rdata", "Full Systematic Map Data.rds"))
systematic <- readRDS(here::here("Rdata", "Systematic-Like Only Extracted Data.rds"))

# critical appraisal load
#CA <- read.csv(here("data", "Critical Appraisal Data.csv"))
#CA <- readRDS(here::here("data", "Critical Analysis Data updated.rds"))
#CA <- read.csv(here::here("data", "Final Critical Appraisal Data.csv"))
CA <- readRDS(here::here("Rdata", "Full Critical Analysis Data.rds"))

# bibliometric data load
#bib <- convert2df(here::here("Data","AcceptedBibs.csv"), dbsource = "scopus", format = "csv")
#bib <- readRDS(here::here("data", "BibliometricData.rds"))
#bibtex <- convert2df(here::here("Data","AcceptedBibs.bib"), dbsource = "scopus", format = "bibtex")
#bib1 <- readRDS(here::here("data", "BibliometricData.rds"))
#bibtex <- readRDS(here::here("data", "BibData_clean.rds"))
#all.res <- readRDS(here::here("data", "Altmetric Data.rds"))
dat_scimago <- readRDS(here::here("data", "scimag_fixed.rds"))
scimag_bib <- readRDS(here::here("data", "bib data with scimag.rds"))
#bibtex <- convert2df(here::here("Data","additional bib info.bib"), dbsource = "scopus", format = "bibtex")
#bib2 <- readRDS(here::here("data", "AdditionalBibliometricData.rds"))
bib <- readRDS(here::here("Rdata", "FullBibliometricData.rds"))
syst.bib <- readRDS(here::here("Rdata", "Systematic-Like Only Bibliometrics.rds"))
all.res <- readRDS(here::here("data", "Full Altmetric Data.rds"))
```


## Visualize Missing Data
```{r}
vis_miss(systmap)
vis_miss(CA)
# data loads in properly
```


# Data Cleaning
```{r}
# Extracted Data
# remove time stamps columns
cols <- colnames(systmap)
systmap2 <- select(systmap, cols[2:35])
#saveRDS(systmap2, here("data", "Systematic Map Data updated.rds"))
# add paper IDs
length(systmap2$Title)
length(unique(systmap2$Title)) #check for duplicates

# Identify duplicates
duplicates <- duplicated(systmap2)
# Remove duplicates
systmap <- systmap2[!duplicates, ]
systmap <- systmap[-c(27, 70, 100, 108, 109, 120, 262, 269, 295, 305), ]

# set paperID numbers
systmap$paperID <- NA
for (i in 1:length(systmap$Title)) {
  systmap$paperID[i] <- paste0("PID_",i)
}

# add in publication Year

# ID article types - they're all merged together when multiple are selected
# articles <- str_split_fixed(systmap2$Article.Type, ", ", 3)
systmap$Narrative.Review <- NA
systmap$Quant.Synthesis <- NA
systmap$Syst.Map <- NA
systmap$MetaAnalysis <- NA
systmap$Syst.Review <- NA
systmap$Critical.Review <- NA
systmap$Map <- NA
systmap$Syst.Lit.Search <- NA
for (i in 1:length(systmap$Article.Type)) {
  if_else(str_detect(systmap$Article.Type[i], "Narrative Review"), 1, 0) -> systmap$Narrative.Review[i]
  if_else(str_detect(systmap$Article.Type[i], "Data/Quantitative Synthesis"), 1, 0) -> systmap$Quant.Synthesis[i]
  if_else(str_detect(systmap$Article.Type[i], "Systematic Map"), 1, 0) -> systmap$Syst.Map[i]
  if_else(str_detect(systmap$Article.Type[i], "Meta-analysis"), 1, 0) -> systmap$MetaAnalysis[i]
  if_else(str_detect(systmap$Article.Type[i], "Systematic Review"), 1, 0) -> systmap$Syst.Review[i]
  if_else(str_detect(systmap$Article.Type[i], "Critical Review"), 1, 0) -> systmap$Critical.Review[i]
  if_else(str_detect(systmap$Article.Type[i], "^Map$"), 1, 0) -> systmap$Map[i]
  if_else(str_detect(systmap$Article.Type[i], "Systematic Literature Search"), 1, 0) -> systmap$Syst.Lit.Search[i]
}

colnames(systmap) <- c("Extractor", "Title", "First_Author", "Article_Type", "Quantitative", "Methodological", "Remote_Sensing", "Monitoring", "Conservation", "Bleaching", "Disease", "Biodiversity_Reef", "Biodiversity_Microbes", "Resilience", "Recovery", "Mortality", "Pathogen", "Microbiome", "Symbiosis", "Climate_Change", "Anthropogenic", "Habitat_Range", "SST", "UV", "Predation", "Acidification", "Calcification", "Bioerosion", "Damage_Wave", "Damage_Human", "Chemical", "Pollution", "Fishing", "Other", "paperID", "Narrative.Review", "Quant.Synthesis", "Syst.Map", "MetaAnalysis", "Syst.Review", "Critical.Review", "Map", "Syst.Lit.Search")
#saveRDS(systmap, here::here("data", "Systematic Map Data.rds"))

# Critical Appraisal Data
colc <- colnames(CA)
CA2 <- select(CA, colc[2:32])
#saveRDS(CA2, here::here("data", "Critical Analysis Data updated.rds"))


# Bibliometric Data
 
# colnames(bibtex) <- c("Authors", "Keywords.Author", "Keywords.Index", "Affil.Location.All", "References", "ISO", "Abstract", "ArticleNo", "CODEN", "Corres.Info", "Type", "DOI", "FundingAgency", "ISBN", "ISSN", "Journal", "Language", "CiteNo.WoS", "PartNo", "Page.Count", "Pages", "Publisher", "PubMed.ID", "Database", "Title", "URL", "Volume", "Year", "Funding.Text", "Author.University", "Author.Dept", "RefCount", "Citation.1", "Citation.2")
# #saveRDS(bibtex, here::here("data", "BibData_clean.rds"))
```

## Just Systematic-Like Papers
```{r}
systematics <- subset(systmap, Syst.Map == 1 | MetaAnalysis == 1 | Syst.Review == 1 | Syst.Lit.Search == 1)
systematic <- systematics[!str_detect(systematics$Article_Type, "Protocol"),]
#saveRDS(systematic, here::here("Rdata", "Systematic-Like Only Extracted Data.rds"))
```


# MetaTag Extraction
```{r}
# add these extracted data to new columns of my own name that's understandable
# but use all original columns for data processing

bibtex2 <- metaTagExtraction(bibtex, Field = "AU1_CO", sep = ";") # Country of affiliation for first author
bibtex2 <- metaTagExtraction(bibtex2, Field = "AU_CO", sep = ";") # Country of affiliation for co-authors
bibtex2 <- metaTagExtraction(bibtex2, Field = "CR_AU", sep = ";") # First author of each cited reference
bibtex2 <- metaTagExtraction(bibtex2, Field = "CR_SO", sep = ";") # Source of each cited reference
bibtex2 <- metaTagExtraction(bibtex2, Field = "AU1_UN", sep = ";") # Uni affiliation for first author
bibtex2 <- metaTagExtraction(bibtex2, Field = "AU_UN", sep = ";") # Uni affiliation for co-authors
bibtex2 <- metaTagExtraction(bibtex2, Field = "SR", sep = ";") # Reference
#saveRDS(bibtex2, here::here("data", "BibliometricData.rds"))
```

# Merge Bibliometric Data
```{r}
all.bib <- full_join(bib1, bib2)

#saveRDS(all.bib, here::here("Rdata", "FullBibliometricData.rds"))
```

## Just Systematic-Like Papers
```{r}
title.matches <- sapply(bib$TI, function(x) {
  max(agrep(x, CA2$Title.of.Paper, max.distance = 0.3, ignore.case = TRUE))
})

syst.bib <- bib[!is.infinite(title.matches), ]
syst.bib <- syst.bib[-c(4, 5, 12, 16), ]
# missing Bove paper since this was not found in Scopus records
#saveRDS(syst.bib, here::here("Rdata", "Systematic-Like Only Bibliometrics.rds"))
```

## Just Systematic-Like Papers that have Code
```{r}
CA_reduced <- select(CA, !ends_with("Comments"))
# Rename columns to be only until first .
colnames(CA_reduced) <- c("TITLE", "Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15")

have.code <- subset(CA_reduced, Q15 == "Yes")

code.matches <- sapply(syst.bib$TI, function(x) {
  max(agrep(x, have.code$TITLE, max.distance = 0.3, ignore.case = TRUE))
})

code.bib <- syst.bib[!is.infinite(code.matches), ]
# missing Bove paper since this was not found in Scopus records
```

## Just Systematic-Like Papers that have Data
```{r}
have.data <- subset(CA_reduced, Q14 == "Yes")

data.matches <- sapply(syst.bib$TI, function(x) {
  max(agrep(x, have.data$TITLE, max.distance = 0.3, ignore.case = TRUE))
})

data.bib <- syst.bib[!is.infinite(data.matches), ]
# missing Bove paper since this was not found in Scopus records
```


# Systematic Mapping
make little tables of (table()) of frequencies of the categorical variables (report year would be a hist or scatterplot)
## Type of review
```{r}
# plot types of reviews
# ggplot(systmap) +
#   geom_bar(aes(x = Article.Type)) +
#   theme_bw()
review.type <- systmap %>% select(c(Narrative.Review, Quant.Synthesis, MetaAnalysis, Syst.Map, Syst.Review, Critical.Review, Map))
review_type <- data.frame(Review = c("Narrative Review", "Quantitative Synthesis", "Meta-Analysis", "Systematic Map", "Systematic Review", "Critical Review", "Map"),
                          Frequency = colSums(review.type))
# count plot
ggplot(review_type, aes(x = Review, y = Frequency)) +
  geom_col() +
  geom_text(aes(label = Frequency), vjust = -0.5) +
  labs(x = "Review Type") +
  theme_bw()

# percentages
as.percent <- function(count, column){
  signif((count/sum(column))*100, digits = 3)
}
review_type$Percent <- lapply(review_type$Frequency, as.percent, column =  review_type$Frequency)
review_type$Percent <- unlist(review_type$Percent)

review_type$Review <- factor(review_type$Review, levels = review_type$Review[order(review_type$Percent, decreasing = TRUE)])

# percent plot
ggplot(review_type, aes(x = Review, y = Percent)) +
  geom_col() +
  geom_text(aes(label = Percent), vjust = -0.5) +
  ylim(0, 100) +
  labs(x = "Review Type", y = "Percent (%)") +
  theme_bw()
```

## Systematic map methodology/purpose data
```{r, height = 15, width = 15}
# plot the yes and no/unclear of each term? - group by section of form on one figure
comp <- count(systmap2$Quantitative.synthesis..Computational.Synthesis..Mathematical.Synthesis.)
percent_comp <- comp %>% mutate(Synthesis = x, percent = (round((freq/sum(freq))*100, digits = 0)))

met1 <- ggplot(percent_comp, aes(x = Synthesis, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Quantitative Synthesis", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


method <- count(systmap2$Methodological)
percent_method <- method %>% mutate(Method = x, percent = (round((freq/sum(freq))*100, digits = 0)))

met2 <- ggplot(percent_method, aes(x = Method, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Methodological", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


remote <- count(systmap2$Remote.Sensing..satellite.sensing..remote.collection..satellite.collection.)
percent_remote <- remote %>% mutate(Sensing = x, percent = (round((freq/sum(freq))*100, digits = 0)))

met3 <- ggplot(percent_remote, aes(x = Sensing, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Remote Sensing", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


mon <- count(systmap2$Monitoring..Surveying..routine.observation..annual.surveys.)
percent_mon <- mon %>% mutate(Monitor = x, percent = (round((freq/sum(freq))*100, digits = 0)))

met4 <- ggplot(percent_mon, aes(x = Monitor, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Monitoring", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


eco <- count(systmap2$Ecosystem.Management.Conservation)
percent_eco <- eco %>% mutate(Conservation = x, percent = (round((freq/sum(freq))*100, digits = 0)))

met5 <- ggplot(percent_eco, aes(x = Conservation, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Conservation", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


met1 + met2 + met3 + met4 + met5 + plot_annotation(tag_levels = "A")



# SYSTEMATIC ONLY ----
comp <- count(systematic$Quantitative)
percent_comp <- comp %>% mutate(Synthesis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method <- count(systematic$Methodological)
percent_method <- method %>% mutate(Method = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

remote <- count(systematic$Remote_Sensing)
percent_remote <- remote %>% mutate(Sensing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

mon <- count(systematic$Monitoring)
percent_mon <- mon %>% mutate(Monitor = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

eco <- count(systematic$Conservation)
percent_eco <- eco %>% mutate(Conservation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method_info <- data.frame(Term = c("Quantitative", "Methodological", "Remote Sensing", "Monitoring", "Conservation"),
                          Yes = c(percent_comp$percent[2], percent_method$percent[2], percent_remote$percent[2], percent_mon$percent[2], percent_eco$percent[2]),
                          No = c(percent_comp$percent[1], percent_method$percent[1], percent_remote$percent[1], percent_mon$percent[1], percent_eco$percent[1]))

method_info <- method_info %>% arrange(Yes)

method_long <- method_info %>% 
  pivot_longer(
    cols = !Term,
    names_to = "Presence",
    values_to = "Percent"
)
method_long$Term <- as.factor(method_long$Term)
method_long$Yes_Percent <- ifelse(method_long$Presence == "Yes",  print(method_long$Percent), print(NA))

extract_method <- ggplot(data = method_long, x = reorder(Term, -Percent), y = Percent) +
  geom_col(mapping = aes(x = reorder(Term, -Percent), y = Percent, fill = Presence), width = 0.7,
           position = "stack", color = "black") +
  scale_x_discrete(limits = rev(levels(method_long$Term))) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.75) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("gainsboro", "black"), name = "Presence:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Extracted Terms - Systematic-Like Papers Only")
```


## Systematic map health outcomes data
```{r, height = 15, width = 15}
# some kind of plot that if (col) is "yes" then print the col_name (this would be per paper - don't really like)

# plot the yes and no/unclear of each term? - group by section of form on one figure
bleach <- count(systmap2$Bleaching)
percent_bleach <- bleach %>% mutate(Bleaching = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out1 <- ggplot(percent_bleach, aes(x = Bleaching, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


dis <- count(systmap2$Disease)
percent_dis <- dis %>% mutate(Disease = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out2 <- ggplot(percent_dis, aes(x = Disease, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


bioc <- count(systmap2$Biodiversity.of.coral.reef.assemblage)
percent_bioc <- bioc %>% mutate(Coral_Biodiv = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out3 <- ggplot(percent_bioc, aes(x = Coral_Biodiv, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Coral Biodiversity", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


biom <- count(systmap2$Biodiversity.of.microbial.community)
percent_biom <- biom %>% mutate(Micro_Biodiv = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out4 <- ggplot(percent_biom, aes(x = Micro_Biodiv, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Microbial Biodiversity", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


RTA <- count(systmap2$Coral.and.or.Coral.Reef.Resilience.Tolerance.Adaptability)
percent_RTA <- RTA %>% mutate(Resilience = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out5 <- ggplot(percent_RTA, aes(x = Resilience, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


RRR <- count(systmap2$Coral.Recovery.Restoration.Rehabilitation)
percent_RRR <- RRR %>% mutate(Recovery = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out6 <- ggplot(percent_RRR, aes(x = Recovery, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


cover <- count(systmap2$Coral.Cover.Mortality)
percent_cover <- cover %>% mutate(Cover = x, percent = (round((freq/sum(freq))*100, digits = 0)))

out7 <- ggplot(percent_cover, aes(x = Cover, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Coral Cover", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()

out1 + out2 + out3 + out4 + out5 + out6 + out7 + plot_annotation(tag_levels = "A")


# SYSTEMATIC ONLY----
bleach <- count(systematic$Bleaching)
percent_bleach <- bleach %>% mutate(Bleaching = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

dis <- count(systematic$Disease)
percent_dis <- dis %>% mutate(Disease = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

bioc <- count(systematic$Biodiversity_Reef)
percent_bioc <- bioc %>% mutate(Coral_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

biom <- count(systematic$Biodiversity_Microbes)
percent_biom <- biom %>% mutate(Micro_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RTA <- count(systematic$Resilience)
percent_RTA <- RTA %>% mutate(Resilience = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RRR <- count(systematic$Recovery)
percent_RRR <- RRR %>% mutate(Recovery = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cover <- count(systematic$Mortality)
percent_cover <- cover %>% mutate(Cover = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

outcomes <- data.frame(Term = c("Bleaching", "Disease", "Reef Biodiversity", "Microbe Biodiversity", "Resilience", "Recovery", "Mortality"),
                          Yes = c(percent_bleach$percent[2], percent_dis$percent[2], percent_bioc$percent[2], percent_biom$percent[2], percent_RTA$percent[2], percent_RRR$percent[2], percent_cover$percent[2]),
                          No = c(percent_bleach$percent[1], percent_dis$percent[1], percent_bioc$percent[1], percent_biom$percent[1], percent_RTA$percent[1], percent_RRR$percent[1], percent_cover$percent[1]))

outcomes <- outcomes %>% arrange(-Yes)

outcome_long <- outcomes %>% 
  pivot_longer(
    cols = !Term,
    names_to = "Presence",
    values_to = "Percent"
)
outcome_long$Term <- as.factor(outcome_long$Term)
outcome_long$Term <- fct_reorder(outcome_long$Term, outcome_long$Percent, max, .desc = TRUE)
outcome_long$Yes_Percent <- ifelse(outcome_long$Presence == "Yes",  print(outcome_long$Percent), print(NA))

ggplot(data = outcome_long, aes(x = reorder(Term, -Percent), y = Percent)) +
  geom_col(mapping = aes(fill = Presence), width = 0.7,
           position = "stack", color = "black") +
  scale_x_discrete(limits = c("Disease", "Microbe Biodiversity", "Reef Biodiversity", "Recovery", "Bleaching", "Mortality", "Resilience")) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("gainsboro", "lightseagreen"), name = "Presence:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Extracted Terms - Systematic-Like Papers Only")
```


## Systematic map driver data
```{r, height = 15, width = 15}
# plot the yes and no/unclear of each term? - group by section of form on one figure
virus <- count(systmap2$Pathogen..virus..bacteria..viral.agent.)
percent_path <- virus %>% mutate(Pathogen = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr1 <- ggplot(percent_path, aes(x = Pathogen, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


micro <- count(systmap2$Microbiome..holobiont..symbionts.)
percent_micro <- micro %>% mutate(Microbiome = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr2 <- ggplot(percent_micro, aes(x = Microbiome, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


sym <- count(systmap2$Symbiosis)
percent_sym <- sym %>% mutate(Symbiosis = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr3 <- ggplot(percent_sym, aes(x = Symbiosis, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


cc <- count(systmap2$Climate.Change)
percent_cc <- cc %>% mutate(Climate = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr4 <- ggplot(percent_cc, aes(x = Climate, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Climate Change", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


AI <- count(systmap2$Anthropogenic.Impacts)
percent_AI <- AI %>% mutate(Anthropogenic = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr5 <- ggplot(percent_AI, aes(x = Anthropogenic, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


hab <- count(systmap2$Habitat.Range.Shift..change.in.geographic.range.)
percent_hab <- hab %>% mutate(Habitat = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr6 <- ggplot(percent_hab, aes(x = Habitat, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Habitat Shift", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


sst <- count(systmap2$Increased.Sea.Surface.Temperature)
percent_sst <- sst %>% mutate(SST = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr7 <- ggplot(percent_sst, aes(x = SST, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Increased SST", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


uv <- count(systmap2$UV.Radiation)
percent_uv <- uv %>% mutate(UV = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr8 <- ggplot(percent_uv, aes(x = UV, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "UV Radiation", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


pred <- count(systmap2$Predation..e.g...invasive.species.)
percent_pred <- pred %>% mutate(Predation = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr9 <- ggplot(percent_pred, aes(x = Predation, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


acid <- count(systmap2$Acidification..acidify..lowered.pH..more.acidic.waters.)
percent_acid <- acid %>% mutate(Acidification = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr10 <- ggplot(percent_acid, aes(x = Acidification, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


calc <- count(systmap2$Calcification..calcify..calcified.)
percent_calc <- calc %>% mutate(Calcification = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr11 <- ggplot(percent_calc, aes(x = Calcification, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


erod <- count(systmap2$Bioerosion..biodegradation...breakdown.of.organic.material.through.bacterial.means.)
percent_ero <- erod %>% mutate(Bioerosion = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr12 <- ggplot(percent_ero, aes(x = Bioerosion, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


wave <- count(systmap2$Physical.Damage...Wave.Action..e.g...storms.cyclones.hurricanes.)
percent_wave <- wave %>% mutate(Wave_Damage = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr13 <- ggplot(percent_wave, aes(x = Wave_Damage, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Physical Damage - Wave", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


hum <- count(systmap2$Physical.Damage...Human..e.g...boating..tourism..land.development..shipwreck.)
percent_hum <- hum %>% mutate(Human_Damage = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr14 <- ggplot(percent_hum, aes(x = Human_Damage, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Physical Damage - Human", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


chem <- count(systmap2$Chemical.Exposure..e.g...oils..hormones..drugs..sewage.exposure..increased.nutrients..herbicides.)
percent_chem <- chem %>% mutate(Chemical = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr15 <- ggplot(percent_chem, aes(x = Chemical, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Chemical Pollution", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


pol <- count(systmap2$Non.chemical.pollution..e.g...sedimentation..eutrophication..solids.)
percent_pol <- pol %>% mutate(Pollution = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr16 <- ggplot(percent_pol, aes(x = Pollution, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(x = "Non-Chemical Pollution", y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


fish <- count(systmap2$Fishing)
percent_fish <- fish %>% mutate(Fishing = x, percent = (round((freq/sum(freq))*100, digits = 0)))

dr17 <- ggplot(percent_fish, aes(x = Fishing, y = percent)) +
  geom_col(aes(), width = 0.5) +
  labs(y = "Percent (%)") +
  ylim(0,100)+
  geom_text(aes(label = percent), vjust = -0.5) +
  theme_bw()


dr1 + dr2 + dr3 + dr4 + dr5 + dr6 + dr7 + dr8 + dr9 + dr10 + dr11 + dr12 + dr13 + dr14 + dr15 + dr16 + dr17 + plot_annotation(tag_levels = "A")


# SYSTEMATIC ONLY ----
virus <- count(systematic$Pathogen)
percent_path <- virus %>% mutate(Pathogen = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

micro <- count(systematic$Microbiome)
percent_micro <- micro %>% mutate(Microbiome = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sym <- count(systematic$Symbiosis)
percent_sym <- sym %>% mutate(Symbiosis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cc <- count(systematic$Climate_Change)
percent_cc <- cc %>% mutate(Climate = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

AI <- count(systematic$Anthropogenic)
percent_AI <- AI %>% mutate(Anthropogenic = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hab <- count(systematic$Habitat_Range)
percent_hab <- hab %>% mutate(Habitat = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sst <- count(systematic$SST)
percent_sst <- sst %>% mutate(SST = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

uv <- count(systematic$UV)
percent_uv <- uv %>% mutate(UV = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pred <- count(systematic$Predation)
percent_pred <- pred %>% mutate(Predation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

acid <- count(systematic$Acidification)
percent_acid <- acid %>% mutate(Acidification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

calc <- count(systematic$Calcification)
percent_calc <- calc %>% mutate(Calcification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

erod <- count(systematic$Bioerosion)
percent_ero <- erod %>% mutate(Bioerosion = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

wave <- count(systematic$Damage_Wave)
percent_wave <- wave %>% mutate(Wave_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hum <- count(systematic$Damage_Human)
percent_hum <- hum %>% mutate(Human_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

chem <- count(systematic$Chemical)
percent_chem <- chem %>% mutate(Chemical = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pol <- count(systematic$Pollution)
percent_pol <- pol %>% mutate(Pollution = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

fish <- count(systematic$Fishing)
percent_fish <- fish %>% mutate(Fishing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))


drivers <- data.frame(Term = c("Pathogen", "Microbiome", "Symbiosis", "Climate Change", "Anthropogenic Impacts", "Range Shift", "Increase SST", "UV Radiation", "Predation", "Acidification", "Calcification", "Bioerosion", "Wave Damage", "Human Damage", "Chemical Exposure", "Pollution", "Fishing"),
                          Yes = c(percent_path$percent[2], percent_micro$percent[2], percent_sym$percent[2], percent_cc$percent[2], percent_AI$percent[2], percent_hab$percent[2], percent_sst$percent[2], percent_uv$percent[2], percent_pred$percent[2], percent_acid$percent[2], percent_calc$percent[2], percent_ero$percent[2], percent_wave$percent[2], percent_hum$percent[2], percent_chem$percent[2], percent_pol$percent[2], percent_fish$percent[2]),
                          No = c(percent_path$percent[1], percent_micro$percent[1], percent_sym$percent[1], percent_cc$percent[1], percent_AI$percent[1], percent_hab$percent[1], percent_sst$percent[1], percent_uv$percent[1], percent_pred$percent[1], percent_acid$percent[1], percent_calc$percent[1], percent_ero$percent[1], percent_wave$percent[1], percent_hum$percent[1], percent_chem$percent[1], percent_pol$percent[1], percent_fish$percent[1]))

drivers <- drivers %>% arrange(-Yes)
term.list <- list(drivers$Term)

drivers_long <- drivers %>% 
  pivot_longer(
    cols = !Term,
    names_to = "Presence",
    values_to = "Percent"
)
drivers_long$Term <- as.factor(drivers_long$Term)
drivers_long$Yes_Percent <- ifelse(drivers_long$Presence == "Yes",  print(drivers_long$Percent), print(NA))

ggplot(data = drivers_long, aes(x = reorder(Term, -Percent), y = Percent)) +
  geom_col(mapping = aes(fill = Presence), width = 0.7,
           position = "stack", color = "black") +
  scale_x_discrete(limits = rev(term.list[[1]])) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("gainsboro", "violetred4"), name = "Presence:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Extracted Terms - Systematic-Like Papers Only")
```

## Additional Terms Added
```{r}
# terms were added in "Other" column of systmap data
OtherData <- subset(systmap, Other != "")
OtherTerms <- unlist(strsplit(OtherData$Other, "; "))
OtherTerms <- unlist(strsplit(OtherTerms, ", "))

# cleaning to get similar terms to match
OtherTerms <- str_to_lower(OtherTerms)
OtherTerm <- sub("s$", "", OtherTerms)
OtherTerm <- sub("herbivory", "herbivore", OtherTerm)
OtherTerm <- sub("^trophic ", "", OtherTerm)
OtherTerm <- sub("^disease ", "", OtherTerm)
OtherTerm <- sub("^larval ", "", OtherTerm)
OtherTerm <- sub("^stress ", "", OtherTerm)
OtherTerm <- sub("^spatial ", "", OtherTerm)
OtherTerm <- sub(" rise$", "", OtherTerm)
OtherTerm <- sub("photosynthesi", "photosynthesis", OtherTerm)

# how many
count_terms <- count(OtherTerm)
AddedTerms <- subset(count_terms, count_terms$freq > 5)

# percentage 
count_terms <- count_terms %>% mutate(as.percent = (freq/335)*100)

NewTerm_df <- data.frame(Term = count_terms$x,
                         Yes = count_terms$as.percent,
                         No = (100 - count_terms$as.percent))

NewTerm_df <- NewTerm_df %>% arrange(-Yes)
Newterm.list <- list(NewTerm_df$Term)

NewTerm_long <- NewTerm_df %>% 
  pivot_longer(
    cols = !Term,
    names_to = "Presence",
    values_to = "Percent"
)

NewTerm_long$Term <- as.factor(NewTerm_long$Term)
NewTerm_long$Yes_Percent <- ifelse(NewTerm_long$Presence == "Yes",  print(NewTerm_long$Percent), print(NA))


# plot
ggplot(data = NewTerm_long, aes(x = reorder(Term, -Percent), y = Percent)) +
  geom_col(mapping = aes(fill = Presence), width = 0.7,
           position = "stack", color = "black") +
  scale_x_discrete(limits = rev(Newterm.list[[1]])) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.05) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("gainsboro", "deepskyblue"), name = "Presence:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Additional Terms")
```


# Critical Appraisal data
```{r}
# selecting only the columns with scores
CA_reduced <- select(CA, !ends_with("Comments"))

# Rename columns to be only until first .
colnames(CA_reduced) <- c("TITLE", "Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15")

ceesat_red <- select(CA_reduced, c("TITLE", "Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13"))
data_code <- select(CA_reduced, c("TITLE", "Q14", "Q15"))

# wrangling data into long format
CA_long <- gather(ceesat_red, question, score, Q1:Q13, factor_key = TRUE)
dc_long <- gather(data_code, question, score, Q14:Q15, factor_key = TRUE)

# calculating the % of scores within each questions 
CA_count <- CA_long %>% dplyr::count(score, by = question) 
CA_percent <- CA_count %>% group_by(by) %>% mutate(percent = (n/35)*100)
CA_percent <- CA_percent %>%
  dplyr::rename(
    question = by
  )
CA_percent$question <- as.factor(CA_percent$question)
CA_percent$question <- factor(CA_percent$question, levels(CA_percent$question)[length(CA_percent$question):1]) #reverse the order of questions
CA_percent$score <- as.factor(CA_percent$score)
CA_percent$score <- factor(CA_percent$score, levels(CA_percent$score)[c(3,1,2)]) #set the order of levels for assessment scores

Amber <- subset(CA_percent, score == "Amber")
(sum(Amber$percent)/13)
Green <- subset(CA_percent, score == "Green")
mean(Green$percent)
Red <- subset(CA_percent, score == "Red")
mean(Red$percent)

dc_count <- dc_long %>% dplyr::count(score, by = question) 
dc_percent <- dc_count %>% mutate(percent = (n/35)*100)
dc_percent <- dc_percent %>%
  dplyr::rename(
    question = by
  )
dc_percent$question <- as.factor(dc_percent$question)
dc_percent$question <- factor(dc_percent$question, levels(dc_percent$question)[length(dc_percent$question):1]) #reverse the order of questions
levels(dc_percent$question) <- c("Code", "Data")
dc_percent$score <- as.factor(dc_percent$score)


# Traffic light plot
ggplot(data = CA_percent, x = question, y = percent) +
  geom_col(mapping = aes(x = question, y = percent, fill = score), width = 0.7,
           position = "fill", color = "black") +
  coord_flip(ylim = c(0, 1)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("red", "yellow", "green"), name = "Score:") +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank()) + 
  ylab("Proportion") + xlab("CEESAT Question")

# Data & Code plot
ggplot(data = dc_percent, x = question, y = percent) +
  geom_col(mapping = aes(x = question, y = percent, fill = score), width = 0.7,
           position = "fill", color = "black") +
  coord_flip(ylim = c(0, 1)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("grey80", "navy"), name = "Score:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank()) +
  ylab("Proportion") + xlab("Code and Data Availability")
```


# Bibliometric Analysis
## Map of first author affiliations
```{r}
#https://www.riinu.me/2022/02/world-map-ggplot2/
#table(bibtex$Corres.Info)[1:10] #column with countries of first authors

A1_country <- bib %>% dplyr::group_by(AU1_CO) %>% dplyr::count(AU1_CO) %>% dplyr::filter(!is.na(AU1_CO))

#load map data
world_map <- map_data("world") %>% 
  filter(! long > 180) #remove countries with longitude >180 to make equal projection-like map without artifacts
table(world_map$region) #note that United Kingdom is UK and America is USA

# Format country names to match regions on the world map
A1_country$region <- str_to_title(A1_country$AU1_CO)
A1_country$region[A1_country$region == "Usa"] <- "USA" #Fix "Usa" to "USA" :
A1_country$region[A1_country$region == "United Kingdom"] <- "UK" #fix to "UK"
(A1_country$region) %in% world_map$region #check matching

## colour all regions on the map:
emptymap <- tibble(region = unique(world_map$region), n = rep(0,length(unique(world_map$region)))) #create table with all counts as 0
fullmap <- left_join(emptymap, A1_country, by = "region") #join with actual counts table
fullmap$count <- fullmap$n.x + fullmap$n.y # make new column for fixed counts
fullmap$count[is.na(fullmap$count)] <- 0 #change NA to 0 for regions with no counts

fullmap %>% 
  ggplot(aes(fill = count, map_id = region)) +
  geom_map(map = world_map) +
  expand_limits(x = world_map$long, y = world_map$lat) +
  coord_map("moll") +
  theme_map(line_size = 0.5) + 
  theme(legend.position="right") +
  scale_fill_gradient(low = "gainsboro", high = "red",
      guide = guide_colorbar(direction = "vertical.")) +
  guides(fill = guide_colourbar(barwidth = unit(15, units = "mm"), barheight = unit(20, units = "mm")))
```
### First Author Affiliations with/without reefs
```{r}
# how many papers are there where the first author countries have reefs?
with.reef <- subset(A1_country, region == "USA" |
                      region == "Australia" |
                      region == "Brazil" |
                      region == "China" |
                      region == "Colombia" |
                      region == "Dominican Republic" |
                      region == "France" |
                      region == "India" |
                      region == "Indonesia" |
                      region == "Israel" |
                      region == "Japan" |
                      region == "Kenya" |
                      region == "Malaysia" |
                      region == "Mauritius" |
                      region == "Mexico" |
                      region == "Philippines" |
                      region == "Saudi Arabia" |
                      region == "South Africa" |
                      region == "Sri Lanka" |
                      region == "Thailand")
sum(with.reef$n)

without.reef <- subset(A1_country, region == "Belgium" |
                      region == "Canada" |
                      region == "Estonia" |
                      region == "Germany" |
                      region == "Hong Kong" |
                      region == "Italy" |
                      region == "Monaco" |
                      region == "Netherlands" |
                      region == "New Zealand" |
                      region == "Portugal" |
                      region == "Sweden" |
                      region == "UK" |
                      region == "Spain")
sum(without.reef$n)

sum(A1_country$n)
reef.perc <- (sum(with.reef$n) / sum(A1_country$n))*100
no.reef.perc <- (sum(without.reef$n) / sum(A1_country$n))*100
```


## Collaboration Network
### Country Collaborations
```{r}
CountryNetwork <- biblioNetwork(bib, analysis = "collaboration", 
                           network = "countries", sep = ";")

collab <- as.matrix(CountryNetwork)
diag(collab) <- 0 #get rid of collaboration with same country

# getting rid of lower triangle (as this is duplication of info)
collab[lower.tri(collab)] <- 0 
# change row and col names case
colnames(collab) <- str_to_title(colnames(collab))
rownames(collab) <- str_to_title(rownames(collab))
# change to "USA"
colnames(collab)[colnames(collab) == "Usa"] <- "USA"
rownames(collab)[rownames(collab) == "Usa"] <- "USA"
# change to "UK"
colnames(collab)[colnames(collab) == "United Kingdom"] <- "UK"
rownames(collab)[rownames(collab) == "United Kingdom"] <- "UK"
# change to "NewZealand"
colnames(collab)[colnames(collab) == "New Zealand"] <- "NewZealand"
rownames(collab)[rownames(collab) == "New Zealand"] <- "NewZealand"
# change to "SaudiArabia"
colnames(collab)[colnames(collab) == "Saudi Arabia"] <- "SaudiArabia"
rownames(collab)[rownames(collab) == "Saudi Arabia"] <- "SaudiArabia"
# change to "SouthAfrica"
colnames(collab)[colnames(collab) == "South Africa"] <- "SouthAfrica"
rownames(collab)[rownames(collab) == "South Africa"] <- "SouthAfrica"
# change to "CostaRica"
colnames(collab)[colnames(collab) == "Costa Rica"] <- "CostaRica"
rownames(collab)[rownames(collab) == "Costa Rica"] <- "CostaRica"
# change to "DominicanRepublic"
colnames(collab)[colnames(collab) == "Dominican Republic"] <- "DominicanRepublic"
rownames(collab)[rownames(collab) == "Dominican Republic"] <- "DominicanRepublic"
# change to "HongKong"
colnames(collab)[colnames(collab) == "Hong Kong"] <- "HongKong"
rownames(collab)[rownames(collab) == "Hong Kong"] <- "HongKong"
# change to "SriLanka"
colnames(collab)[colnames(collab) == "Sri Lanka"] <- "SriLanka"
rownames(collab)[rownames(collab) == "Sri Lanka"] <- "SriLanka"

circos.clear()
colours.random <- c(USA = "black", Australia = "black", UK = "darkslategray", China = "black", Israel = "gainsboro", Germany = "gray", Brazil = "lightgrey", Canada = "slategray", France = "lightslategray", Mexico = "lightslategray", Netherlands = "darkgray", Sweden = "darkslategray", SaudiArabia = "grey80", SouthAfrica = "gainsboro", Colombia = "grey66", Italy = "dimgray", Kenya = "grey80", Philippines = "darkslategray", Portugal = "grey66", India = "slategray", Fiji = "darkslategray", Japan = "darkslategray", Malaysia = "slategray", Monaco = "lightgrey", NewZealand = "gainsboro", Denmark = "lightgrey", Norway = "grey66", Panama = "grey80", Singapore = "slategray", Switzerland = "gainsboro", Thailand = "black", Venezuela = "slategray", Austria = "lightgrey", Belgium = "grey66", CostaRica = "darkgray", DominicanRepublic = "darkslategray", Guam = "dimgray", Indonesia = "black", Samoa = "grey", Spain = "lightgrey", Tanzania = "lightgrey", Mauritius = "darkslategray", Greece = "lightslategray", Egypt = "darkslategray", Ireland = "gray", HongKong = "black")

# reference: https://hdr.undp.org/system/files/documents/global-report-document/hdr2021-22pdf_1.pdf (above medium HDI = developed)
colours.develop <- c(USA = "black", Australia = "black", UK = "black", China = "dimgray", Israel = "black", Germany = "black", Brazil = "dimgray", Canada = "black", France = "black", Mexico = "dimgray", Netherlands = "black", Sweden = "black", SaudiArabia = "black", SouthAfrica = "dimgray", Colombia = "dimgray", Italy = "black", Kenya = "gray", Philippines = "gray", Portugal = "black", India = "gray", Fiji = "dimgray", Japan = "black", Malaysia = "black", Monaco = "lightgrey", NewZealand = "black", Denmark = "black", Norway = "black", Panama = "black", Singapore = "black", Switzerland = "black", Thailand = "black", Venezuela = "gray", Austria = "black", Belgium = "black", CostaRica = "black", DominicanRepublic = "dimgray", Guam = "black", Indonesia = "dimgray", Samoa = "dimgray", Spain = "black", Tanzania = "gray", Mauritius = "black", Greece = "black", Egypt = "dimgray", Ireland = "black", HongKong = "black")

# reference: https://coral.org/en/coral-reefs-101/geography/, https://habitats.oceanplus.org/ (restricted to warm water corals)
colours.coral <- c(USA = "steelblue1", Australia = "blue", UK = "gray", China = "steelblue1", Israel = "steelblue1", Germany = "gray", Brazil = "steelblue1", Canada = "gray", France = "steelblue1", Mexico = "steelblue1", Netherlands = "gray", Sweden = "gray", SaudiArabia = "steelblue1", SouthAfrica = "steelblue1", Colombia = "steelblue1", Italy = "gray", Kenya = "steelblue1", Philippines = "blue", Portugal = "gray", India = "steelblue1", Fiji = "blue", Japan = "steelblue1", Malaysia = "steelblue1", Monaco = "gray", NewZealand = "gray", Denmark = "gray", Norway = "gray", Panama = "steelblue1", Singapore = "steelblue1", Switzerland = "gray", Thailand = "steelblue1", Venezuela = "steelblue1", Austria = "gray", Belgium = "gray", CostaRica = "steelblue1", DominicanRepublic = "steelblue1", Guam = "steelblue1", Indonesia = "blue", Samoa = "steelblue1", Spain = "gray", Tanzania = "steelblue1", Mauritius = "steelblue1", Greece = "steelblue1", Egypt = "steelblue1", Ireland = "gray", HongKong = "gray")


C.network <- chordDiagram(collab,
                          annotationTrack = "grid",
                          preAllocateTracks = 1,
                          grid.col = colours.coral)
C.network <- select(C.network, c(rn, cn, value1, col))
C.network <- C.network[order(-C.network$value1),]
C.network <- C.network %>% group_by(value1)

chordDiagram(C.network,
             annotationTrack = "grid",
             preAllocateTracks = 1,
             grid.col = colours.coral)

circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
  circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
  circos.axis(h = "top", labels.cex = 0.5, major.tick.length = 0.2, sector.index = sector.name, track.index = 2)
}, bg.border = NA)
C.network
```

### Author Collabs
```{r}
AuthorNetwork <- biblioNetwork(bib, analysis = "collaboration", network = "authors", sep = ";")


networkPlot(AuthorNetwork, Title = "Author Collaborations", cluster = "infomap")
networkPlot(AuthorNetwork, Title = "Author Collaborations", cluster = "infomap", label = FALSE, halo = TRUE)
networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, halo = TRUE, remove.isolates = TRUE)
networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE)
networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE, alpha = 0.75, community.repulsion = 0)


# try n or *degree* to see if you can remove satellites
networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE, degree = 5)
networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE, degree = 7, n = 325)

# how many authors in the central node?
str(AuthorNetwork)    # what does this data look like?
AuthorNetwork[1:10, 1:10]   
colSums(as.matrix(AuthorNetwork) == 0)   # number of connections total
hist(colSums(as.matrix(AuthorNetwork) == 0))
1190 - colSums(as.matrix(AuthorNetwork) == 0)   # number of connections with other authors
hist(1190 - colSums(as.matrix(AuthorNetwork) == 0))
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 5)   # 709 authors have more than 5 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 2)   # 966 authors have more than 2 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 10)  # 481 authors have more than 10 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 70)  # 97 authors have more than 70 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 80)  # 96 authors have more than 80 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 100) # 4 authors have more than 100 connections
sum(1190-colSums(as.matrix(AuthorNetwork) == 0) > 50)  # 103 authors have more than 50 connections

# ## using igraph
# collab <- graph_from_incidence_matrix(AuthorNetwork)
# collab <- simplify(collab)
# layout <- layout_with_fr(collab)
# 
# #plot(collab, vertex.label = NA, vertex.size = 5, edge.arrow.size = 0.5)
# ggraph(collab, layout = layout) + 
#   geom_edge_link() + 
#   geom_node_point(size = 1, alpha = 0.5, colour= 'coral') +
#   #geom_jitter() +
#   xlab("") + 
#   ylab("") + 
#   theme_void()
```
old code
```{r, include=FALSE}
AuthorNetwork <- biblioNetwork(bib, analysis = "collaboration", 
                           network = "authors", sep = ";")
a.collab <- as.matrix(AuthorNetwork)
diag(a.collab) <- 0
a.collab[lower.tri(a.collab)] <- 0
colnames(a.collab) <- str_to_title(colnames(a.collab))
rownames(a.collab) <- str_to_title(rownames(a.collab))
colnames(a.collab) <- str_remove_all(colnames(a.collab), "[:space:]")
rownames(a.collab) <- str_remove_all(rownames(a.collab), "[:space:]")

author.count <- colSums(a.collab)
author.count <- as.table(author.count) # current issue is that I can tell there are multiple names for one author here (e.g., Ainsworth T and Ainsworth TD)
author.count <- as.data.frame(author.count)
author.count <- author.count %>% subset(Freq != 0)
author.count <- arrange(author.count, desc(Freq))

A.network <- chordDiagram(a.collab,
                          annotationTrack = "grid",
                          preAllocateTracks = 1,
                          grid.col = palette.colors())
chordDiagram(A.network,
             annotationTrack = c("name", "grid"),
             #keep.diagonal = TRUE,
             #preAllocateTracks = 1,
             grid.col = palette.colors())



# what is going on with authors
df <- data.frame(bib$AU)
df$AU2

# checking how many single author papers
is.na(str_locate(df$bib.AU, ";")[,1])

df.au <- data.frame()
for (i in 1:nrow(df)) {
  if (!is.na(str_locate(df$bib.AU, ";")[,1][i])) {
    df.au <- rbind(df.au, df$bib.AU[i])
    print(df[i,])
  }
}

# list of unique author names - str_??
df_new <- df.au
authors <- unique(word(df_new$X.ACOSTA.A.DUEÑAS.LF.PIZARRO.V., sep = ";"))

# create dataframe col1 = author names (from above), col2 = 0
author_df <- data.frame(author = authors,
                        tally = 0)

# for loop (i 1:length(df)) if str_?? = true then +1 to matching col name
for (i in 1:nrow(author_df)) {
  for (j in 1:nrow(df_new)){
    if (str_detect(df_new$X.ACOSTA.A.DUEÑAS.LF.PIZARRO.V.[j], author_df$author[i])){
      author_df$tally[i] <- author_df$tally[i] +1
    }
    }
}

author_df$tally <- sort(author_df$tally, decreasing = TRUE)
print(author_df[1:15,])


# network function in bibliometrix instead. Don't use chordDiagram
```

## Most used Keywords word cloud
```{r}
# isolate keywords
keyword.data <- bib$DE

# create a corpus for cleaning text
corpkeys <- Corpus(VectorSource(keyword.data))

# cleaning
keywords <- corpkeys %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
keywords <- tm_map(keywords, content_transformer(tolower))
keywords <- tm_map(keywords, removeWords, stopwords("english"))

# create dataframe of words and freq
words <- TermDocumentMatrix(keywords) 
word.matrix <- as.matrix(words) 
wordcounts <- sort(rowSums(word.matrix),decreasing=TRUE) 
keyword.df <- data.frame(keywords = names(wordcounts),freq = wordcounts)

# word cloud
set.seed(1234) # for reproducibility
wordcloud(words = keyword.df$keywords, freq = keyword.df$freq, min.freq = 1, max.words=200, random.order=FALSE, rot.per=0.35, colors=brewer.pal(8, "Dark2"))
```

## Published Journals
```{r}
journals <- count(bib$SO)
colnames(journals) <- c("Journal", "Freq")
journals <- arrange(journals, -Freq)
# journals <- filter(journals, Freq > 2)
journals <- journals[1:10,]
iso <- count(bib$JI)
colnames(iso) <- c("Journal", "Freq")
# iso <- filter(iso, Freq > 2)
iso <- arrange(iso, -Freq)
iso <- iso[1:10,]

ggplot(journals) +
  geom_col(aes(x = reorder(Journal, Freq), y = Freq)) +
  #scale_x_discrete(guide = guide_axis(n.dodge = 3)) +
  labs(x = "Journal", y = "Frequency") +
  coord_flip() +
  theme_bw() #+
  #theme(axis.text.x = element_text(angle = 300, hjust = 0))

ggplot(iso) +
  geom_col(aes(x = reorder(Journal, -Freq), y = Freq)) +
  #scale_x_discrete(guide = guide_axis(n.dodge = 3)) +
  labs(x = "Journal", y = "Frequency") +
  #coord_flip() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 300, hjust = 0))
```

Are the reviews published in special or general focus journals?
```{r}
length(unique(bib$SO))
unique(bib$SO)
specialized <- bib$SO[str_which(bib$SO, paste(c("MARINE", "CORAL", "DISEASE", "COAST", "AQUA", "OCEAN", "HYDRO", "WATER"), collapse = "|"))] # why does this not return all instances? many indices are missing
count.special <- length(specialized)
general <- bib$SO[str_which(bib$SO, paste(c("MARINE", "CORAL", "DISEASE", "COAST", "AQUA", "OCEAN", "HYDRO", "WATER"), collapse = "|"), negate = TRUE)]
count.gen <- length(general)
length(bib$SO)
journal.type <- data.frame(Type = c("General", "Specialized"),
                           Count = c(count.gen, count.special))
journal.type$percent <- NA
journal.type$percent[1] <- (count.gen/sum(journal.type$Count))*100
journal.type$percent[2] <- (count.special/sum(journal.type$Count))*100
journal.type$percent <- signif(journal.type$percent, digits = 3)


ggplot(journal.type) +
  geom_col(aes(x = Type, y = Count)) +
  labs(x = "Journal Type", y = "Frequency") +
  theme_bw()

ggplot(journal.type) +
  geom_col(aes(x = Type, y = percent)) +
  labs(x = "Journal Type", y = "Percent (%)") +
  ylim(0,100) +
  geom_text(aes(x = Type, y = percent, label = percent), vjust = -0.5) +
  theme_bw()


# USING SCIMAG - code adapted from Losia
bib_ <- bib
dat_scimago <- readr::read_delim(file = here::here("data","scimagojr 2021.csv"), delim = ";")
attr(dat_scimago,'problems') <- NULL #remove problems attribute from the IF column
#clean-up Scimago journal names: 
dat_scimago$Journal_name <- str_replace_all(dat_scimago$Title, "[\r\n]" , " ") %>% str_to_upper() #remove newline symbols and change all to uppercase
dat_scimago$Journal_name <- gsub("&", "AND", dat_scimago$Journal_name, fixed = TRUE) #replace all & with AND
#check the overlap of journal names in two data frames
fit <- bib_$SO %in% dat_scimago$Journal_name
sum(fit) # found in scimago
sum(!fit) # not found in scimago # 15 papers
table(bib_$SO[!fit]) # look at these online at scimagojr.com journal search

# find journal names in scimago dataframe that don't match bib dataframe
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "BIOLOGY REPORTS")]  # this journal doesn't exist anymore so it's not in the 2021 collection of journals in scimago. add in manually
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "ENVIRONMENTAL MONITORING")] # cannot find this journal at all on scimag. Can find on research gate only: https://www.researchgate.net/journal/Journal-of-Environmental-Monitoring-1464-0333
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "JOURNAL OF MARINE SCIENCE AND TECHNOLOGY")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "JOURNAL OF TOXICOLOGY AND ENVIRONMENTAL HEALTH")] # this journal doesn't exist anymore so it's not in the 2021 collection of journals in scimago. add in manually
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "LIFE SUPPORT & BIOSPHERE SCIENCE")] # this journal doesn't exist anymore so it's not in the 2021 collection of journals in scimago. add in manually
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "MARINE ECOLOGY - PROGRESS SERIES")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "OCEANOGRAPHY AND MARINE BIOLOGY")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "ACTA ECOLOGICA SINICA")] # alternate name for Shengtai xuebao
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "ZHURNAL OBSH")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "COASTAL RESEARCH")] # book series
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "EARTH AND ENVIRONMENT")] # not in scimag
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "ENVIRONMENTAL SCIENCE")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "NATURAL RESOURCE")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "PETROCHEMICAL")] # this journal doesn't exist anymore so it's not in the 2021 collection of journals in scimago. add in manually
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "POLLUTION")]
#dat_scimago$Journal_name[str_detect(dat_scimago$Journal_name, "SCIENTIFIC WORLD")]

# clean sicmag data for matching
dat_scimago$Journal_name <- str_replace(dat_scimago$Journal_name, "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY. PART A, MOLECULAR ANDAMP; INTEGRATIVE PHYSIOLOGY", "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY - A MOLECULAR AND INTEGRATIVE PHYSIOLOGY")
dat_scimago$Journal_name <- str_replace(dat_scimago$Journal_name, "ENVIRONMENTAL SCIENCE ANDAMP; TECHNOLOGY", "ENVIRONMENTAL SCIENCE AND TECHNOLOGY")
dat_scimago$Journal_name <- str_replace(dat_scimago$Journal_name, "POLLUTION RESEARCH \\(DISCONTINUED\\)", "POLLUTION RESEARCH")
dat_scimago$Journal_name <- str_replace(dat_scimago$Journal_name, "SCIENTIFIC WORLD JOURNAL, THE", "THE SCIENTIFIC WORLD JOURNAL")

# clean bib journal names to match scimago data
bib_$SO <- str_replace(bib_$SO, "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY - A MOLECULAR AND INTEGRATIVE PHYSIOLOGY", "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY. PART A, MOLECULAR ANDAMP; INTEGRATIVE PHYSIOLOGY")
bib_$SO <- str_replace(bib_$SO, "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY. PART A, MOLECULAR ANDAMP; INTEGRATIVE PHYSIOLOGY", "COMPARATIVE BIOCHEMISTRY AND PHYSIOLOGY - A MOLECULAR AND INTEGRATIVE PHYSIOLOGY")
bib_$SO <- str_replace(bib_$SO, "MARINE ECOLOGY PROGRESS SERIES", "MARINE ECOLOGY - PROGRESS SERIES")
bib_$SO <- str_replace(bib_$SO, "OCEANOGRAPHY AND MARINE BIOLOGY: AN ANNUAL REVIEW", "OCEANOGRAPHY AND MARINE BIOLOGY")
bib_$SO <- str_replace(bib_$SO, "SHENGTAI XUEBAO/ ACTA ECOLOGICA SINICA", "SHENGTAI XUEBAO")
bib_$SO <- str_replace(bib_$SO, "SHENGTAI XUEBAO", "ACTA ECOLOGICA SINICA")
bib_$SO <- str_replace(bib_$SO, bib_$SO[str_detect(bib_$SO, "ZHURNAL OBSH")], "ZHURNAL OBSHCHEI BIOLOGII")
bib_$SO <- str_replace(bib_$SO, "NATURAL RESOURCE MODELING", "NATURAL RESOURCE MODELLING")
bib_$SO <- str_replace(bib_$SO, "SCIENCE \\(NEW YORK, N\\.Y\\.\\)", "SCIENCE")
bib_$SO <- str_replace(bib_$SO, "THESCIENTIFICWORLDJOURNAL", "THE SCIENTIFIC WORLD JOURNAL")

dat_scimago$Journal <- NA
for (i in 1:length(dat_scimago$Journal_name)) {
  ifelse(dat_scimago$Journal_name[i] == "JOURNAL OF MARINE SCIENCE AND TECHNOLOGY" & dat_scimago$Country[i] == "Taiwan",
         dat_scimago$Journal[i] <- "JOURNAL OF MARINE SCIENCE AND TECHNOLOGY (TAIWAN)", 
         dat_scimago$Journal[i] <- dat_scimago$Journal_name[i])
}
#saveRDS(dat_scimago, here::here("data", "scimag_fixed.rds"))

fit2 <- bib_$SO %in% dat_scimago$Journal
sum(fit2) # found in scimago
sum(!fit2) # not found in scimago
table(bib_$SO[!fit2])


dat_scimago <- dat_scimago %>% select(-c("Title", "Sourceid", "Type", "Issn", "Total Docs. (2021)", "Total Docs. (3years)", "Total Refs.", "Total Cites (3years)", "Citable Docs. (3years)", "Cites / Doc. (2years)",  "Ref. / Doc.", "Publisher", "Coverage"))

#merge data frames to append journal info to bib dataframe
bib_$Journal <- bib_$SO
scimag_bib <- left_join(bib_, dat_scimago, by = c("Journal"))
#str(scimag_bib)

scimag_bib$SJR <- as.numeric(gsub(",", ".", scimag_bib$SJR)) #change to numeric scores

#write_csv(scimag_bib, file = here::here("data","bib data with scimag.csv"))
#saveRDS(scimag_bib, file = here::here("data","bib data with scimag.rds"))
#table(scimag_bib$Categories) #one or more Categories per journal - need to separate
categories_list <- str_split(scimag_bib$Categories, "; ") #split column

##How many categories per journal? - interdisciplinarity
#count number of sub-elements in a list, for each journal
count_categories <- lengths(categories_list) #NA is counted as 1
count_categories[is.na(categories_list)] <- "no information" #change NA to text value and all to strings
categories_list_counts <- as.data.frame(table(count_categories)) ##count each value
names(categories_list_counts) <- c("subject_categories", "article_count") #rename variables 
levels(categories_list_counts$subject_categories) <- c("1 subject category", "2 subject categories", "3 subject categories", "4 subject categories", "5 subject categories", "6 subject categories", "7 subject categories", "8 subject categories", "no information") #rename levels

#plot
Colors_JCatL <-  c("1 subject category" = "black", "2 subject categories" = "lightcoral", "3 subject categories" = "lightcoral", "4 subject categories" = "lightcoral", "5 subject categories" = "lightcoral", "6 subject categories" = "lightcoral", "7 subject categories" = "lightcoral", "8 subject categories" = "lightcoral", "no information" = "#bdbdbd") #custom colors
cat.count <- ggplot(categories_list_counts, aes(x = subject_categories, y = article_count, fill = subject_categories)) + 
  geom_bar(stat = "identity") +
  coord_flip() +
  xlab("")  +
  scale_fill_manual(values = Colors_JCatL) +
  #geom_text(aes(label = article_count), color="white", size = 3, position = position_stack(vjust = 0.5)) + 
  scale_x_discrete(limits = rev(levels(categories_list_counts$subject_categories))) + #reverse order of the bars (factor)
  ylab("article count") +
  theme_bw() +
  theme(legend.position = "none")


##Most common journal categories? - dominating research subjects
#unlist and clean categories by removing last 5 characters (Q-values) from each sub-element
categories_unlisted <- gsub('.{5}$', '', unlist(categories_list))
count_categories_unlisted <- as.data.frame(table(categories_unlisted)) %>% slice_max(order_by = Freq, n = 10) #get top10 most common subject categories

#plot
Colors_JCatUL <-  c("Aquatic Science" = "seagreen3",  "Ecology, Evolution, Behavior and Systematics" = "seagreen3", "Oceanography" = "seagreen3", "Ecology" = "seagreen3", "Environmental Science (miscellaneous)" = "seagreen3", rep.int("darkgrey", 7)) #custom colors
cat.names <- ggplot(count_categories_unlisted, aes(x = reorder(categories_unlisted, Freq), y = Freq, fill = categories_unlisted)) + 
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(name = '', values = Colors_JCatUL) + #, breaks = c('Multidisciplinary'), labels = c("journal subject categories other than Biology/Ecology/Environment")) +
    #geom_text(aes(label = article_count), color = "white", size = 3, position = position_stack(vjust = 0.5)) +
  xlab("") +
  ylab("article count") +
  theme_bw() + 
  theme(legend.position = "none")

cat.count / cat.names + plot_annotation(tag_levels = "A")
```

## Publication Year
```{r}
ggplot(bib) +
  geom_bar(aes(x = PY)) +
  labs(x = "Publication Year", y = "Number of Articles") +
  theme_bw()

ggplot(count(bib$PY)) +
  geom_line(aes(x = x, y = freq)) +
  labs(x = "Publication Year", y = "Number of Articles") +
  theme_bw()
```

### Publication Year with Country
```{r}
ggplot(bib) +
  geom_histogram(aes(x = PY, colour = AU1_CO), fill = "gainsboro") + 
  theme_bw()

```


## Citations
```{r}
ggplot(bib) +      # Not pretty, too large a gap
  geom_bar(aes(x = TC)) +
  labs(x = "Number of Citations", y = "Number of Articles") +
  theme_bw()
min(bib$TC)
max(bib$TC)

# try bins?

ggplot(bib) +
  geom_jitter(aes(x = PY, y = TC)) +
  labs(x = "Publication Year", y = "Number of Citations") +
  theme_bw()


# average citations per Year (not publication year)
bib.dat <- biblioAnalysis(bib)

cit.df <- data.frame(Year = bib.dat$Years,
                     Cits = bib.dat$TCperYear)
cit.df$Year<- as.factor(cit.df$Year) 
citation.df <- aggregate(cit.df$Cits, list(cit.df$Year), mean)
colnames(citation.df) <- c("Year", "Citation_Num")

ggplot(citation.df) +
  geom_col(aes(x = Year, y = Citation_Num))+
  labs(x = "Publication Year", y = "Citations per Article") +
  theme_bw()

# 
# cit.df2 <- data.frame(Year = bib.dat$Years,
#                      Cits = bib.dat$TotalCitation)
# cit.df2$Year<- as.factor(cit.df$Year) 
# citation.df2 <- aggregate(cit.df2$Cits, list(cit.df2$Year), mean)
# colnames(citation.df2) <- c("Year", "Citation_Num")
# 
# ggplot(citation.df2) +
#   geom_col(aes(x = Year, y = Citation_Num))+
#   theme_bw()

plot(bib.dat, k = 10, pause = F)
```

## Altmetrics
```{r}
# rAltmetric doesn't work, Yefeng Yang made a new function

# function getAltmetrics(), the only parameter is doi
getAltmetrics <- function(doi = NULL,
                          foptions = list(),
                           ...) {
    if (!is.null(doi)) doi <- stringr::str_c("doi/", doi)
    identifiers <- purrr::compact(list(doi))
    if (!is.null(identifiers)) {
      ids <- identifiers[[1]]
    }
    base_url <- "http://api.altmetric.com/v1/"
    #request <- httr::GET(paste0(base_url, ids), httr::add_headers("user-agent" = "#rstats rAltmertic package https://github.com/ropensci/rAltmetric"))
    request <- httr::GET(paste0(base_url, ids))
    results <-
      jsonlite::fromJSON(httr::content(request, as = "text"), flatten = TRUE)
    results <- rlist::list.flatten(results)
    class(results) <- "altmetric"
    results
}

altmetric.crawler <- list(NULL)

for (n in 1:length(bib$DI)) { # change your file name
# format altmetric object
  format.Altmetric <- function(altmetric.object) {
  stats <- altmetric.object[grep("^cited", names(altmetric.object))]
  stats <- data.frame(stats, stringsAsFactors = FALSE)
  data.frame(paper_title = altmetric.object$title,
             journal = altmetric.object$journal,
             doi = altmetric.object$doi,
             #subject = altmetric.object$subjects,
             Altmetric.score = altmetric.object$score,
             stats = stats)
}
# JASON formate
  altmetric.crawler[[n]] <- suppressMessages(try(list(format.Altmetric(getAltmetrics(doi = bib$DI[n]))))) #  # change your file name
 
# create a dataframe function
  altmetric_df <- function(altmetric.object) {
  df <- data.frame(t(unlist(altmetric.object)), stringsAsFactors = FALSE)
  }
  #altmetric.crawler[[n]]  <-  try(list(altmetric_df(getAltmetrics(doi = DOIs[n]))))
  # create a function to summarize Altmetric object
  summary.altmetric <- function(x, ...) {
  if (inherits(x, "altmetric"))  {
string <- "Altmetrics on: \"%s\" with altmetric_id: %s published in %s."
vals   <- c(x$title,  x$altmetric_id, x$journal)
 if("journal" %in% names(x)) {
  cat(do.call(sprintf, as.list(c(string, vals))))
 } else {
   string <- "Altmetrics on: \"%s\" with altmetric_id: %s"
   cat(do.call(sprintf, as.list(c(string, vals))))
 }
  cat("\n")
  stats <- x[grep("^cited", names(x))]
  stats <- data.frame(stats, stringsAsFactors = FALSE)
  print(data.frame(stats = t(stats)))
  }
}
  # crawl
 # altmetric.crawler[[n]] <- try(list(summary.altmetric(getAltmetrics(doi = DOIs[n]))))
}


# save results from altmetric.crawler and retrieve lists within lists
altmetric.crawler2 <- sapply(altmetric.crawler, function(x) {x})


# retrieve stats
all.res <- data.frame(paper_title = sapply(altmetric.crawler2, function(x)  ifelse(class(x) == "data.frame",x$paper_title,NA)),
           journal = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",x$journal,NA)),
           doi = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",x$doi,NA)),
           #subject = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",x$subject,NA)),
           Altmetric.score = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",x$Altmetric.score,0)),
           policy = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_policies_count),x$stats.cited_by_policies_count,0),0)),
           patent = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_patents_count),x$stats.cited_by_patents_count,0),0)),
           feeds = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_feeds_count),x$stats.cited_by_feeds_count,0),0)),
           news = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_msm_count),x$stats.cited_by_msm_count,0),0)),
           posts = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_posts_count),x$stats.cited_by_posts_count,0),0)),
           twitter = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_tweeters_count),x$stats.cited_by_tweeters_count,0),0)),
           wiki = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_wikipedia_count),x$stats.cited_by_wikipedia_count,0),0)),
           videos = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_videos_count),x$stats.cited_by_videos_count,0),0)),
           accounts = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_accounts_count),x$stats.cited_by_accounts_count,0),0)),
           fb = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_fbwalls_count),x$stats.cited_by_fbwalls_count,0),0)),
           rh = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_rh_count),x$stats.cited_by_rh_count,0),0)),
           rdts = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_rdts_count),x$stats.cited_by_rdts_count,0),0)),
           googlep = sapply(altmetric.crawler2, function(x) ifelse(class(x) == "data.frame",ifelse(!is.null(x$stats.cited_by_gplus_count),x$stats.cited_by_gplus_count,0),0))
           )
# saveRDS(all.res, here::here("data", "Altmetric Data.rds"))
# write.csv(all.res, here::here("data", "Altmetric Data.csv"))

all.res <- subset(all.res, is.na(doi) == FALSE)

alt.means <- colMeans(select(all.res, !c(1:3)))
colMax <- function(data) sapply(data, max, na.rm = TRUE)
alt.max <- colMax(select(all.res, !c(1:3)))

# https://api.altmetric.com/docs/call_citations.html - what all but "accounts" means

AS <- ggplot(all.res, aes(x = Altmetric.score)) +
  geom_histogram() +
  labs(x = "Altmetric Score", y = "Frequency") +
  theme_bw()

pol <- ggplot(all.res, aes(x = policy)) +
  geom_histogram() +
  labs(x = "Policy Sources", y = "Frequency") +
  theme_bw()

pat <- ggplot(all.res, aes(x = patent)) +
  geom_histogram() +
  labs(x = "Patents", y = "Frequency") +
  xlim(-0.5,5) +
  theme_bw()

feed <- ggplot(all.res, aes(x = feeds)) +
  geom_histogram() +
  labs(x = "Blogs", y = "Frequency") +
  theme_bw()

msm <- ggplot(all.res, aes(x = news)) +
  geom_histogram() +
  labs(x = "News Outlets", y = "Frequency") +
  theme_bw()

post <- ggplot(all.res, aes(x = posts)) +
  geom_histogram() +
  labs(x = "Total Posts", y = "Frequency") +
  theme_bw()

twit <- ggplot(all.res, aes(x = twitter)) +
  geom_histogram() +
  labs(x = "Twitter Users", y = "Frequency") +
  theme_bw()

wikip <- ggplot(all.res, aes(x = wiki)) +
  geom_histogram() +
  labs(x = "Wikipedia Pages", y = "Frequency") +
  theme_bw()

vid <- ggplot(all.res, aes(x = videos)) +
  geom_histogram() +
  labs(x = "YouTube Channels", y = "Frequency") +
  xlim(-0.5,5) +
  theme_bw()

acc <- ggplot(all.res, aes(x = accounts)) +
  geom_histogram() +
  labs(x = "Accounts ?", y = "Frequency") +
  theme_bw()

FB <- ggplot(all.res, aes(x = fb)) +
  geom_histogram() +
  labs(x = "Facebook Accounts", y = "Frequency") +
  theme_bw()

RH <- ggplot(all.res, aes(x = rh)) +
  geom_histogram() +
  labs(x = "Research Highlight Platforms", y = "Frequency") +
  xlim(-0.5,5) +
  theme_bw()

red <- ggplot(all.res, aes(x = rdts)) +
  geom_histogram() +
  labs(x = "Reddit Users", y = "Frequency") +
  xlim(-0.5,5) +
  theme_bw()

ggl <- ggplot(all.res, aes(x = googlep)) +
  geom_histogram() +
  labs(x = "Google+ Users", y = "Frequency") +
  xlim(-0.5,5) +
  theme_bw()

(AS + post + acc + msm + feed + wikip + twit +  FB + vid + red + ggl + RH + pol + pat) + plot_annotation(tag_levels = "A")
```

## Impact Factor
```{r}
clarivate_IF <- read.csv(here::here("data", "SamanthaBurke_JCR_JournalResults_01_2023.csv"))
clarivate_IF$Journal.name <- stringr::str_to_upper(clarivate_IF$Journal.name)

bib_wIF <- as.data.frame(bib$SO)
colnames(bib_wIF) <- c("Journal.name")
matches <- left_join(bib_wIF, clarivate_IF, by = "Journal.name")
missing_journals <- subset(matches, is.na(X2021.JIF) == TRUE)
unique(missing_journals$Journal.name)

missing_IFs <- data.frame(Journal.name = unique(missing_journals$Journal.name),
                          missing_2021.JIF = c(NA, 4.653, 2.166, 4.640, 3.229, 4.328, 10.618, 3.644, 5.247, NA, 4.926, 3.752, 3.307, 2.476, 2.744, 1.022, 2.381, 3.109, 4.315, 3.172, 3.031, 3.737, 7.378, 2.134, 3.167, NA, 0.553, 4.755, NA, 2.646, 3.173, 0.37, 2.945, NA, 0.462, 9.470, 4.177, NA, 4.379, 3.061, 5.092, 2.698, 0.15, 0.746, NA, NA, 1.932, 2.179, 2.822, 3.605, 1.274, 2.915, 4.416, 1.805, NA, 4.218, 3.737, 3.392, 20.589, 1.559, 14.340, 13.757, 1.769, 14.284, 0.803, 9.711, 0.896, 8.071, 4.295, 3.906, 2.795, 5.135, 2.941, 1.254, 3.842, 6.000, 5.492, 0.306, 0.917, 3.734, 6.671, 4.181, 37.214, NA, 2.358, 4.474, 5.019, 3.530, 1.871, 7.290, 3.084, 2.200, 12.038, 1.609, 0.42, 4.498, 1.394, 1.736, 4.296, NA, 4.520, 5.349, NA, NA, 1.689, 3.298, NA, 1.110, 2.453, 63.832, 2.170, NA, 11.357, 1.219, 3.925, NA))


full_match <- left_join(matches, missing_IFs, by = "Journal.name")
full_match$IF <- coalesce(full_match$X2021.JIF,full_match$missing_2021.JIF)
journal_wIF <- full_match %>% select(Journal.name, ISSN, eISSN, IF)
colnames(journal_wIF) <- c("SO", "SN", "EI", "IF")

mean(na.omit(journal_wIF$IF))
```

## H-Index
```{r}
H_Index <- na.omit(scimag_bib$`H index`)
max(H_Index)
min(H_Index)
mean(H_Index)
```


# Figure 2 - systematic map for methodology information

## Figure 2A - review type
```{r}
review.type <- systmap %>% select(c(Narrative.Review, Quant.Synthesis, MetaAnalysis, Syst.Map, Syst.Review, Critical.Review, Syst.Lit.Search))
review_type <- data.frame(Review = c("Narrative Review", "Quantitative Synthesis", "Meta-Analysis", "Systematic Map", "Systematic Review", "Critical Review", "Systematic Literature Search"),
                          Frequency = colSums(review.type))

as.percent <- function(count, column){
  signif((count/sum(column))*100, digits = 3)
}
review_type$Percent <- lapply(review_type$Frequency, as.percent, column =  review_type$Frequency)
review_type$Percent <- unlist(review_type$Percent)
# percent plot
article.type <- ggplot(review_type, aes(x = reorder(Review, Percent), y = Percent)) +
  geom_col(fill = "black") +
  geom_text(aes(label = Percent), hjust = -0.45) +
  coord_flip(ylim = c(0, 100)) +
  labs(x = "Review Type", y = "Percent (%)") +
  theme_bw() +
  theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12))#, axis.text.x = element_text(angle = -45, vjust = 0.05, hjust = -0.01))
```

## Figure 2B - extracted methodology information
```{r}
comp <- count(systmap$Quantitative)
percent_comp <- comp %>% mutate(Synthesis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method <- count(systmap$Methodological)
percent_method <- method %>% mutate(Method = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

remote <- count(systmap$Remote_Sensing)
percent_remote <- remote %>% mutate(Sensing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

mon <- count(systmap$Monitoring)
percent_mon <- mon %>% mutate(Monitor = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

eco <- count(systmap$Conservation)
percent_eco <- eco %>% mutate(Conservation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method_info <- data.frame(Term = c("Quantitative", "Methodological", "Remote Sensing", "Monitoring", "Conservation"),
                          Yes = c(percent_comp$percent[2], percent_method$percent[2], percent_remote$percent[2], percent_mon$percent[2], percent_eco$percent[2]),
                          No = c(percent_comp$percent[1], percent_method$percent[1], percent_remote$percent[1], percent_mon$percent[1], percent_eco$percent[1]))

method_info$Term <- as.factor(method_info$Term)

extract_met <- ggplot(method_info) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "black") +
  scale_x_discrete(limits = rev(levels(method_info$Term))) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5)+ 
  ylab("Percent (%)") + xlab("Extracted Terms") +
  theme_bw() +
  theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12)) 


# method_info <- method_info %>% arrange(Yes)
# 
# method_long <- method_info %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# method_long$Term <- as.factor(method_long$Term)
# method_long$Yes_Percent <- ifelse(method_long$Presence == "Yes",  print(method_long$Percent), print(NA))
# 
# extract_method <- ggplot(data = method_long, x = Term, y = Percent) +
#   geom_col(mapping = aes(x = Term, y = Percent, fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = rev(levels(method_long$Term))) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.75) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "black"), name = "Presence:") +
#   theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("Percent (%)") + xlab("Extracted Terms")
```

## Combine plots
```{r}
(article.type / extract_met) + plot_annotation(tag_levels = "A")
```


# Figure 3 - systematic map information
## Figure 3A - drivers
```{r}
virus <- count(systmap$Pathogen)
percent_path <- virus %>% mutate(Pathogen = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

micro <- count(systmap$Microbiome)
percent_micro <- micro %>% mutate(Microbiome = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sym <- count(systmap$Symbiosis)
percent_sym <- sym %>% mutate(Symbiosis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cc <- count(systmap$Climate_Change)
percent_cc <- cc %>% mutate(Climate = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

AI <- count(systmap$Anthropogenic)
percent_AI <- AI %>% mutate(Anthropogenic = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hab <- count(systmap$Habitat_Range)
percent_hab <- hab %>% mutate(Habitat = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sst <- count(systmap$SST)
percent_sst <- sst %>% mutate(SST = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

uv <- count(systmap$UV)
percent_uv <- uv %>% mutate(UV = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pred <- count(systmap$Predation)
percent_pred <- pred %>% mutate(Predation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

acid <- count(systmap$Acidification)
percent_acid <- acid %>% mutate(Acidification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

calc <- count(systmap$Calcification)
percent_calc <- calc %>% mutate(Calcification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

erod <- count(systmap$Bioerosion)
percent_ero <- erod %>% mutate(Bioerosion = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

wave <- count(systmap$Damage_Wave)
percent_wave <- wave %>% mutate(Wave_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hum <- count(systmap$Damage_Human)
percent_hum <- hum %>% mutate(Human_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

chem <- count(systmap$Chemical)
percent_chem <- chem %>% mutate(Chemical = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pol <- count(systmap$Pollution)
percent_pol <- pol %>% mutate(Pollution = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

fish <- count(systmap$Fishing)
percent_fish <- fish %>% mutate(Fishing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))


drivers <- data.frame(Term = c("Pathogen", "Microbiome", "Symbiosis", "Climate Change", "Anthropogenic Impacts", "Range Shift", "Increase SST", "UV Radiation", "Predation", "Acidification", "Calcification", "Bioerosion", "Wave Damage", "Human Damage", "Chemical Exposure", "Pollution", "Fishing"),
                          Yes = c(percent_path$percent[2], percent_micro$percent[2], percent_sym$percent[2], percent_cc$percent[2], percent_AI$percent[2], percent_hab$percent[2], percent_sst$percent[2], percent_uv$percent[2], percent_pred$percent[2], percent_acid$percent[2], percent_calc$percent[2], percent_ero$percent[2], percent_wave$percent[2], percent_hum$percent[2], percent_chem$percent[2], percent_pol$percent[2], percent_fish$percent[2]),
                          No = c(percent_path$percent[1], percent_micro$percent[1], percent_sym$percent[1], percent_cc$percent[1], percent_AI$percent[1], percent_hab$percent[1], percent_sst$percent[1], percent_uv$percent[1], percent_pred$percent[1], percent_acid$percent[1], percent_calc$percent[1], percent_ero$percent[1], percent_wave$percent[1], percent_hum$percent[1], percent_chem$percent[1], percent_pol$percent[1], percent_fish$percent[1]))

drivers <- drivers %>% arrange(-Yes)
term.list <- list(drivers$Term)

drivers$Term <- as.factor(drivers$Term)

drivers$Term <- reorder(drivers$Term, drivers$Yes)

driver_plot <- ggplot(drivers) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "violetred4") +
  scale_x_discrete(limits = levels(drivers$Term)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Extracted Terms") +
  theme_bw()

# drivers_long <- drivers %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# drivers_long$Term <- as.factor(drivers_long$Term)
# drivers_long$Yes_Percent <- ifelse(drivers_long$Presence == "Yes",  print(drivers_long$Percent), print(NA))
# 
# driver_plot <- ggplot(data = drivers_long, aes(x = reorder(Term, -Percent), y = Percent)) +
#   geom_col(mapping = aes(fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = rev(term.list[[1]])) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "violetred4"), name = "Presence:") +
#   theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("Percent (%)") + xlab("Extracted Terms")

```

## Figure 3B - outcomes
```{r}
bleach <- count(systmap$Bleaching)
percent_bleach <- bleach %>% mutate(Bleaching = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

dis <- count(systmap$Disease)
percent_dis <- dis %>% mutate(Disease = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

bioc <- count(systmap$Biodiversity_Reef)
percent_bioc <- bioc %>% mutate(Coral_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

biom <- count(systmap$Biodiversity_Microbes)
percent_biom <- biom %>% mutate(Micro_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RTA <- count(systmap$Resilience)
percent_RTA <- RTA %>% mutate(Resilience = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RRR <- count(systmap$Recovery)
percent_RRR <- RRR %>% mutate(Recovery = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cover <- count(systmap$Mortality)
percent_cover <- cover %>% mutate(Cover = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

outcomes <- data.frame(Term = c("Bleaching", "Disease", "Reef Biodiversity", "Microbe Biodiversity", "Resilience", "Recovery", "Mortality"),
                          Yes = c(percent_bleach$percent[2], percent_dis$percent[2], percent_bioc$percent[2], percent_biom$percent[2], percent_RTA$percent[2], percent_RRR$percent[2], percent_cover$percent[2]),
                          No = c(percent_bleach$percent[1], percent_dis$percent[1], percent_bioc$percent[1], percent_biom$percent[1], percent_RTA$percent[1], percent_RRR$percent[1], percent_cover$percent[1]))

outcomes <- outcomes %>% arrange(-Yes)

outcomes$Term <- as.factor(outcomes$Term)
outcomes$Term <- reorder(outcomes$Term, outcomes$Yes)

out_plot <- ggplot(outcomes) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "lightseagreen") +
  scale_x_discrete(limits = levels(outcomes$Term)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Extracted Terms") +
  theme_bw()

# outcome_long <- outcomes %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# outcome_long$Term <- as.factor(outcome_long$Term)
# outcome_long$Term <- fct_reorder(outcome_long$Term, outcome_long$Percent, max, .desc = TRUE)
# outcome_long$Yes_Percent <- ifelse(outcome_long$Presence == "Yes",  print(outcome_long$Percent), print(NA))
# 
# out_plot <- ggplot(data = outcome_long, aes(x = Term, y = Percent)) +
#   geom_col(mapping = aes(fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = c("Microbe Biodiversity", "Disease", "Reef Biodiversity", "Mortality", "Recovery", "Bleaching", "Resilience")) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "lightseagreen"), name = "Presence:") +
#   theme(legend.position = "bottom", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("Percent (%)") + xlab("Extracted Terms")
```

## Combine plots
```{r}
( driver_plot / out_plot) + 
  plot_layout(widths = c(1, 1), heights = c(1, 0.75)) +
  plot_annotation(tag_levels = "A")
```


# Figure 4 - bibliometric map collaboration information

## Figure 4A - map of first author affiliations
```{r}
A1_country <- bib %>% dplyr::group_by(AU1_CO) %>% dplyr::count(AU1_CO) %>% dplyr::filter(!is.na(AU1_CO))

#load map data
world_map <- map_data("world") %>% 
  filter(! long > 180) #remove countries with longitude >180 to make equal projection-like map without artifacts

# Format country names to match regions on the world map
A1_country$region <- str_to_title(A1_country$AU1_CO)
A1_country$region[A1_country$region == "Usa"] <- "USA" #Fix "Usa" to "USA" :
A1_country$region[A1_country$region == "United Kingdom"] <- "UK" #fix to "UK"

## colour all regions on the map:
emptymap <- tibble(region = unique(world_map$region), n = rep(0,length(unique(world_map$region)))) #create table with all counts as 0
fullmap <- left_join(emptymap, A1_country, by = "region") #join with actual counts table
fullmap$Article.Count <- fullmap$n.x + fullmap$n.y # make new column for fixed counts
fullmap$Article.Count[is.na(fullmap$Article.Count)] <- 0 #change NA to 0 for regions with no counts

authormap <- fullmap %>% 
  ggplot(aes(fill = Article.Count, map_id = region)) +
  geom_map(map = world_map) +
  expand_limits(x = world_map$long, y = world_map$lat) +
  coord_map("moll") +
  theme_map(line_size = 10) + 
  theme(text = element_text(family = "sans"), legend.text = element_text("Number of Articles Published")) +
  scale_fill_gradient(low = "grey94", high = "darkblue",
      guide = guide_colorbar(direction = "vertical.")) +
  guides(fill = guide_colourbar(barwidth = unit(10, units = "mm"), barheight = unit(30, units = "mm")))
```

## Figure 4B - Country collaborations
```{r}
CountryNetwork <- biblioNetwork(bib, analysis = "collaboration", 
                           network = "countries", sep = ";")

collab <- as.matrix(CountryNetwork)
diag(collab) <- 0 #get rid of collaboration with same country

# getting rid of lower triangle (as this is duplication of info)
collab[lower.tri(collab)] <- 0 
# change row and col names case
colnames(collab) <- str_to_title(colnames(collab))
rownames(collab) <- str_to_title(rownames(collab))
# change to "USA"
colnames(collab)[colnames(collab) == "Usa"] <- "USA"
rownames(collab)[rownames(collab) == "Usa"] <- "USA"
# change to "UK"
colnames(collab)[colnames(collab) == "United Kingdom"] <- "UK"
rownames(collab)[rownames(collab) == "United Kingdom"] <- "UK"
# change to "NewZealand"
colnames(collab)[colnames(collab) == "New Zealand"] <- "NewZealand"
rownames(collab)[rownames(collab) == "New Zealand"] <- "NewZealand"
# change to "SaudiArabia"
colnames(collab)[colnames(collab) == "Saudi Arabia"] <- "SaudiArabia"
rownames(collab)[rownames(collab) == "Saudi Arabia"] <- "SaudiArabia"
# change to "SouthAfrica"
colnames(collab)[colnames(collab) == "South Africa"] <- "SouthAfrica"
rownames(collab)[rownames(collab) == "South Africa"] <- "SouthAfrica"
# change to "CostaRica"
colnames(collab)[colnames(collab) == "Costa Rica"] <- "CostaRica"
rownames(collab)[rownames(collab) == "Costa Rica"] <- "CostaRica"
# change to "DominicanRepublic"
colnames(collab)[colnames(collab) == "Dominican Republic"] <- "DominicanRepublic"
rownames(collab)[rownames(collab) == "Dominican Republic"] <- "DominicanRepublic"
# change to "HongKong"
colnames(collab)[colnames(collab) == "Hong Kong"] <- "HongKong"
rownames(collab)[rownames(collab) == "Hong Kong"] <- "HongKong"
# change to "SriLanka"
colnames(collab)[colnames(collab) == "Sri Lanka"] <- "SriLanka"
rownames(collab)[rownames(collab) == "Sri Lanka"] <- "SriLanka"

# reference: https://coral.org/en/coral-reefs-101/geography/, https://habitats.oceanplus.org/ (restricted to warm water corals)
# colours.coral <- c(USA = "steelblue1", Australia = "blue", UK = "gray", China = "steelblue1", Israel = "steelblue1", Germany = "gray", Brazil = "steelblue1", Canada = "gray", France = "steelblue1", Mexico = "steelblue1", Netherlands = "gray", Sweden = "gray", SaudiArabia = "steelblue1", SouthAfrica = "steelblue1", Colombia = "steelblue1", Italy = "gray", Kenya = "steelblue1", Philippines = "blue", Portugal = "gray", India = "steelblue1", Fiji = "blue", Japan = "steelblue1", Malaysia = "steelblue1", Monaco = "gray", NewZealand = "gray", Denmark = "gray", Norway = "gray", Panama = "steelblue1", Singapore = "steelblue1", Switzerland = "gray", Thailand = "steelblue1", Venezuela = "steelblue1", Austria = "gray", Belgium = "gray", CostaRica = "steelblue1", DominicanRepublic = "steelblue1", Guam = "steelblue1", Indonesia = "blue", Samoa = "steelblue1", Spain = "gray", Tanzania = "steelblue1", Mauritius = "steelblue1", Greece = "steelblue1", Egypt = "steelblue1", Ireland = "gray", HongKong = "gray")
colours.coral <- c(USA = "steelblue1", Australia = "blue", UK = "goldenrod1", China = "steelblue1", Israel = "steelblue1", Germany = "gray", Brazil = "steelblue1", Canada = "gray", France = "goldenrod1", Mexico = "steelblue1", Netherlands = "goldenrod1", Sweden = "gray", SaudiArabia = "steelblue1", SouthAfrica = "gray", Colombia = "steelblue1", Italy = "gray", Kenya = "steelblue1", Philippines = "blue", Portugal = "gray", India = "steelblue1", Fiji = "blue", Japan = "steelblue1", Malaysia = "steelblue1", Monaco = "gray", NewZealand = "gray", Denmark = "gray", Norway = "gray", Panama = "steelblue1", Singapore = "steelblue1", Switzerland = "gray", Thailand = "steelblue1", Venezuela = "steelblue1", Austria = "gray", Belgium = "gray", CostaRica = "steelblue1", DominicanRepublic = "steelblue1", Guam = "steelblue1", Indonesia = "blue", Samoa = "steelblue1", Spain = "gray", Tanzania = "steelblue1", Mauritius = "steelblue1", Greece = "gray", Egypt = "steelblue1", Ireland = "gray", HongKong = "gray")

C.network <- chordDiagram(collab,
                          annotationTrack = "grid",
                          preAllocateTracks = 1,
                          grid.col = colours.coral)
C.network <- select(C.network, c(rn, cn, value1, col))
C.network <- C.network[order(-C.network$value1),]
C.network <- C.network %>% group_by(value1)

chordDiagram(C.network,
             annotationTrack = "grid",
             preAllocateTracks = 1,
             grid.col = colours.coral)

circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
  circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
}, bg.border = NA)

circos.clear()
```

# Figure 5 - author collaboration central node
```{r}
AuthorNetwork <- biblioNetwork(bib, analysis = "collaboration", network = "authors", sep = ";")

a.net <- networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE, degree = 7, n = 325)
```


# Figure 6
## Figure 6A - publication timeline
```{r}
timeline <- ggplot(bib) +
  geom_bar(aes(x = PY)) +
  labs(x = "Year", y = "Number of Articles Published") +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 6B - Citation
```{r}
outlier <- which(bib$TC > 2000)
bib$label <- NA
bib$label[outlier[1]] <- print("Hughes et al., 2003")
bib$label[outlier[2]] <- print("Hoegh-Guldberg et al., 2007")

citation <- ggplot(bib) +
  geom_jitter(aes(x = PY, y = TC), size = 3, alpha = 0.25) +
  labs(x = "Publication Year", y = "Number of Citations") +
  geom_label(aes(x = PY, y = TC, label = label), hjust = -0.1, vjust = 0.4) +
  geom_hline(yintercept = mean(bib$TC), col = "deepskyblue", linewidth = 1.5, alpha = 0.5) +
  theme(element_text(size = 12)) +
  theme_bw()

```

## Combine plots
```{r}

(timeline / citation) + plot_annotation(tag_levels = "A")
```


# Figure 7 - Journal information
## Figure 7A - subject area
```{r}
categories_list <- str_split(scimag_bib$Categories, "; ") #split column

##How many categories per journal? - interdisciplinarity
#count number of sub-elements in a list, for each journal
count_categories <- lengths(categories_list) #NA is counted as 1
count_categories[is.na(categories_list)] <- "no information" #change NA to text value and all to strings
categories_list_counts <- as.data.frame(table(count_categories)) ##count each value
names(categories_list_counts) <- c("subject_categories", "article_count") #rename variables 
levels(categories_list_counts$subject_categories) <- c("1 Subject Category", "2 Subject Categories", "3 Subject Categories", "4 Subject Categories", "5 Subject Categories", "6 Subject Categories", "7 Subject Categories", "8 Subject Categories", "no information") #rename levels

#plot
Colors_JCatL <-  c("1 Subject Category" = "black", "2 Subject Categories" = "lightcoral", "3 Subject Categories" = "lightcoral", "4 Subject Categories" = "lightcoral", "5 Subject Categories" = "lightcoral", "6 Subject Categories" = "lightcoral", "7 Subject Categories" = "lightcoral", "8 Subject Categories" = "lightcoral", "no information" = "#bdbdbd") #custom colors
cat.count <- ggplot(categories_list_counts, aes(x = subject_categories, y = article_count, fill = subject_categories)) + 
  geom_bar(stat = "identity") +
  coord_flip() +
  xlab("Number of Subjects")  +
  scale_fill_manual(values = Colors_JCatL) +
  #geom_text(aes(label = article_count), color="white", size = 3, position = position_stack(vjust = 0.5)) + 
  scale_x_discrete(limits = rev(levels(categories_list_counts$subject_categories))) + #reverse order of the bars (factor)
  ylab("Article Count") +
  theme_bw() +
  theme(legend.position = "none")

```

## Figure 7B - common subjects
```{r}
#unlist and clean categories by removing last 5 characters (Q-values) from each sub-element
categories_unlisted <- gsub('.{5}$', '', unlist(categories_list))
count_categories_unlisted <- as.data.frame(table(categories_unlisted)) %>% slice_max(order_by = Freq, n = 10) #get top10 most common subject categories

Colors_JCatUL <-  c("Aquatic Science" = "seagreen3",  "Ecology, Evolution, Behavior and Systematics" = "seagreen3", "Oceanography" = "seagreen3", "Ecology" = "seagreen3", "Environmental Science (miscellaneous)" = "seagreen3", rep.int("darkgrey", 7)) #custom colors
cat.names <- ggplot(count_categories_unlisted, aes(x = reorder(categories_unlisted, Freq), y = Freq, fill = categories_unlisted)) + 
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(name = '', values = Colors_JCatUL) + #, breaks = c('Multidisciplinary'), labels = c("journal subject categories other than Biology/Ecology/Environment")) +
    #geom_text(aes(label = article_count), color = "white", size = 3, position = position_stack(vjust = 0.5)) +
  xlab("Category Subjects") +
  ylab("Article Count") +
  theme_bw() + 
  theme(legend.position = "none")
```

## Figure 7C - top journals
```{r}
iso <- count(bib$JI)
colnames(iso) <- c("Journal", "Freq")
# iso <- filter(iso, Freq > 2)
iso <- arrange(iso, -Freq)
iso <- iso[1:10,]

journal.name <- ggplot(iso) +
  geom_col(aes(x = reorder(Journal, Freq), y = Freq)) +
  #scale_x_discrete(guide = guide_axis(n.dodge = 3)) +
  labs(x = "Journal", y = "Article Count") +
  coord_flip() +
  theme_bw()
```

## Combine plots
```{r}
cat.count / cat.names / journal.name + plot_annotation(tag_levels = "A")
```


# Figure 8 - Altmetrics

## Figure 8A - total altmetric score
```{r}
all.res <- subset(all.res, is.na(doi) == FALSE)

AS <- ggplot(all.res, aes(x = Altmetric.score)) +
  geom_histogram() +
  labs(x = "Altmetric Score", y = "Frequency") +
  #geom_vline(xintercept = mean(all.res$Altmetric.score), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 8B - policy
```{r}
pol <- ggplot(all.res, aes(x = policy)) +
  geom_histogram() +
  labs(x = "Policy Sources", y = "Article Count") +
  #geom_vline(xintercept = mean(all.res$policy), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 8C - patent
```{r}
pat <- ggplot(all.res, aes(x = patent)) +
  geom_histogram() +
  labs(x = "Patents", y = "Article Count") +
  #geom_vline(xintercept = mean(all.res$patent), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 8D - mainstream news
```{r}
msm <- ggplot(all.res, aes(x = news)) +
  geom_histogram() +
  labs(x = "News Outlets", y = "Article Count") +
  #geom_vline(xintercept = mean(all.res$news), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 8E - Wikipedia
```{r}
wikip <- ggplot(all.res, aes(x = wiki)) +
  geom_histogram() +
  labs(x = "Wikipedia Pages", y = "Article Count") +
  #geom_vline(xintercept = mean(all.res$wiki), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Figure 8F - twitter
```{r}
twit <- ggplot(all.res, aes(x = twitter)) +
  geom_histogram() +
  labs(x = "Twitter Users", y = "Article Count") +
  #geom_vline(xintercept = mean(all.res$twitter), col = "deepskyblue", linewidth = 1.5, alpha = 0.75) +
  theme(element_text(size = 12)) +
  theme_bw()
```

## Combine plots
```{r}
((AS +  pol + pat)/ (msm + wikip + twit)) + plot_annotation(tag_levels = "A")
```


# Figure 9 - Critical Appraisal

## Figure 9A - summary plot
```{r}
# selecting only the columns with scores
CA_reduced <- select(CA, !ends_with("Comments"))

# Rename columns to be only until first .
colnames(CA_reduced) <- c("TITLE", "Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15")

ceesat_red <- select(CA_reduced, c("TITLE", "Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13"))

# wrangling data into long format
CA_long <- gather(ceesat_red, question, score, Q1:Q13, factor_key = TRUE)

# calculating the % of scores within each questions 
CA_count <- CA_long %>% dplyr::count(score, by = question) 
CA_percent <- CA_count %>% mutate(percent = ((n/35)*100))
CA_percent <- CA_percent %>%
  dplyr::rename(
    question = by
  )
CA_percent$question <- as.factor(CA_percent$question)
CA_percent$question <- factor(CA_percent$question, levels(CA_percent$question)[length(CA_percent$question):1]) #reverse the order of questions
CA_percent$score <- as.factor(CA_percent$score)
CA_percent$score <- factor(CA_percent$score, levels(CA_percent$score)[c(3,1,2)]) #set the order of levels for assessment scores

sum.plot <- ggplot(data = CA_percent, x = question, y = percent) +
  geom_col(mapping = aes(x = question, y = percent, fill = score), width = 0.7,
           position = "stack", color = "black") +
  coord_flip(ylim = c(0, 100)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("red", "gold1", "green"), name = "Score:") +
  theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percentage (%) of Systematic-like Papers") + xlab("CEESAT Question")

```

## Figure 9B - data and code plot
```{r}
data_code <- select(CA_reduced, c("TITLE", "Q14", "Q15"))
dc_long <- gather(data_code, question, score, Q14:Q15, factor_key = TRUE)

dc_count <- dc_long %>% dplyr::count(score, by = question) 
dc_percent <- dc_count %>% mutate(percent = signif((n/35)*100, digits = 3))
dc_percent <- dc_percent %>%
  dplyr::rename(
    question = by
  )

dc_percent <- subset(dc_percent, dc_percent$score == "Yes")

dc_percent$question <- as.factor(dc_percent$question)
dc_percent$question <- factor(dc_percent$question, levels(dc_percent$question)[length(dc_percent$question):1]) #reverse the order of questions
levels(dc_percent$question) <- c("Code", "Data")
#dc_percent$score <- as.factor(dc_percent$score)


dc.plot <- ggplot(dc_percent) +
  geom_col(aes(x = question, y = percent), width = 0.7, color = "black", fill = "navy") +
  scale_x_discrete(limits = levels(dc_percent$question)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = question, y = percent, label = percent), hjust = -0.5) + 
  ylab("Percentage (%) of Systematic-like Papers") + xlab("Code and Data Availability") +
  theme_bw() +
  theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12))

# dc.plot <- ggplot(data = dc_percent, x = question, y = percent) +
#   geom_col(mapping = aes(x = question, y = percent, fill = score), width = 0.7,
#            position = "stack", color = "black") +
#   coord_flip(ylim = c(0, 100)) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("grey80", "navy"), name = "Score:") +
#   theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank()) + 
#   ylab("Percentage (%) of Systematic-like Papers") + xlab("Code and Data Availability")
```

## Combine plots
```{r}
(sum.plot / dc.plot) + 
  plot_layout(widths = c(1, 1), heights = c(1, 0.5)) + 
  plot_annotation(tag_levels = "A")
```


# Supplementary Figures
## Figure S1 - extracted terms through time
```{r}
title.matches2 <- sapply(systmap$Title, function(x) {
  max(agrep(x, bib$TI, max.distance = 0.3, ignore.case = TRUE))
})

systmap$Year <- NA
for (i in 1:nrow(systmap)){
  systmap$Year[i] <- bib$PY[title.matches2[i]]
}
```

### Figure S1A - methodology
```{r}
comp.time <- subset(systmap, Quantitative == "Yes")
comp.count <- count(comp.time$Year)

method.time <- subset(systmap, Methodological == "Yes")
method.count <- count(method.time$Year)

remote.time <- subset(systmap, Remote_Sensing == "Yes")
remote.count <- count(remote.time$Year)

mon.time <- subset(systmap, Monitoring == "Yes")
mon.count <- count(mon.time$Year)

eco.time <- subset(systmap, Conservation == "Yes")
eco.count <- count(eco.time$Year)

met.year.df <- data.frame(Term = c(rep("Quantitative", nrow(comp.count)), rep("Methodological", nrow(method.count)), rep("Remote Sensing", nrow(remote.count)), rep("Monitoring", nrow(mon.count)), rep("Conservation", nrow(eco.count))),
                      Year = c(comp.count[,1], method.count[,1], remote.count[,1], mon.count[,1], eco.count[,1]),
                      Count = c(comp.count[,2], method.count[,2], remote.count[,2], mon.count[,2], eco.count[,2]))

palette.colors(palette = "Okabe-Ito")
met.time.plot <- ggplot(met.year.df) +
  geom_line(aes(x = Year, y = Count, color = Term)) +
  scale_colour_manual(values = c("black", "orange", "skyblue", "#CC79A7", "#009E73")) +
  #scale_color_viridis(discrete = TRUE) +
  ylab("Number of Articles") +
  theme(text = element_text(size = 12)) +
  theme_bw()
```

### Figure S1B - drivers
```{r}
virus.time <- subset(systmap, Pathogen == "Yes")
virus.count <- count(virus.time$Year)

micro.time <- subset(systmap, Microbiome == "Yes")
micro.count <- count(micro.time$Year)

sym.time <- subset(systmap, Symbiosis == "Yes")
sym.count <- count(sym.time$Year)

cc.time <- subset(systmap, Climate_Change == "Yes")
cc.count <- count(cc.time$Year)

AI.time <- subset(systmap, Anthropogenic == "Yes")
AI.count <- count(AI.time$Year)

hab.time <- subset(systmap, Habitat_Range == "Yes")
hab.count <- count(hab.time$Year)

sst.time <- subset(systmap, SST == "Yes")
sst.count <- count(sst.time$Year)

uv.time <- subset(systmap, UV == "Yes")
uv.count <- count(uv.time$Year)

pred.time <- subset(systmap, Predation == "Yes")
pred.count <- count(pred.time$Year)

acid.time <- subset(systmap, Acidification == "Yes")
acid.count <- count(acid.time$Year)

calc.time <- subset(systmap, Calcification == "Yes")
calc.count <- count(calc.time$Year)

erod.time <- subset(systmap, Bioerosion == "Yes")
erod.count <- count(erod.time$Year)

wave.time <- subset(systmap, Damage_Wave == "Yes")
wave.count <- count(wave.time$Year)

hum.time <- subset(systmap, Damage_Human == "Yes")
hum.count <- count(hum.time$Year)

chem.time <- subset(systmap, Chemical == "Yes")
chem.count <- count(chem.time$Year)

pol.time <- subset(systmap, Pollution == "Yes")
pol.count <- count(pol.time$Year)

fish.time <- subset(systmap, Fishing == "Yes")
fish.count <- count(fish.time$Year)

drive.year.df <- data.frame(Term = c(rep("Pathogen", nrow(virus.count)), rep("Microbiome", nrow(micro.count)), rep("Symbiosis", nrow(sym.count)), rep("Climate Change", nrow(cc.count)), rep("Anthropogenic Impacts", nrow(AI.count)), rep("Range Shift", nrow(hab.count)), rep("Increase SST", nrow(sst.count)), rep("UV Radiation", nrow(uv.count)), rep("Predation", nrow(pred.count)), rep("Acidification", nrow(acid.count)), rep("Calcification", nrow(calc.count)), rep("Bioerosion", nrow(erod.count)), rep("Wave Damage", nrow(wave.count)), rep("Human Damage", nrow(hum.count)), rep("Chemical Exposure", nrow(chem.count)), rep("Pollution", nrow(pol.count)), rep("Fishing", nrow(fish.count))),
                      Year = c(virus.count[,1], micro.count[,1], sym.count[,1], cc.count[,1], AI.count[,1], hab.count[,1], sst.count[,1], uv.count[,1], pred.count[,1], acid.count[,1], calc.count[,1], erod.count[,1], wave.count[,1], hum.count[,1], chem.count[,1], pol.count[,1], fish.count[,1]),
                      Count = c(virus.count[,2], micro.count[,2], sym.count[,2], cc.count[,2], AI.count[,2], hab.count[,2], sst.count[,2], uv.count[,2], pred.count[,2], acid.count[,2], calc.count[,2], erod.count[,2], wave.count[,2], hum.count[,2], chem.count[,2], pol.count[,2], fish.count[,2]))

drive.year.df$Term <- as.factor(drive.year.df$Term)

drive.time.plot <- ggplot(drive.year.df) +
  geom_line(aes(x = Year, y = Count, color = Term, linetype = Term)) +
  scale_colour_manual(values = rep(c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999"), 2)) +     # too many terms for this palette
  scale_linetype_manual(values = rep(c("solid", "dashed"), each = 9)) +
  ylab("Number of Articles") +
  theme(text = element_text(size = 12), legend.title = "Term") +
  theme_bw()

# this plot is so chaotic
# splitting into 3-4 categories for visual purposes
cc.drive.year.df <- data.frame(Term = as.factor(c(rep("Climate Change", nrow(cc.count)), rep("Increase SST", nrow(sst.count)), rep("UV Radiation", nrow(uv.count)), rep("Acidification", nrow(acid.count)), rep("Calcification", nrow(calc.count)), rep("Bioerosion", nrow(erod.count)))),
                               Year = c(cc.count[,1], sst.count[,1], uv.count[,1], acid.count[,1], calc.count[,1], erod.count[,1]),
                               Count = c(cc.count[,2], sst.count[,2], uv.count[,2], acid.count[,2], calc.count[,2], erod.count[,2]))
ccterm_count <- aggregate(Count ~ Year, data = cc.drive.year.df, sum)
ccterm_count$Category <- print("Climate Change")

ai.drive.year.df <- data.frame(Term = as.factor(c(rep("Anthropogenic Impacts", nrow(AI.count)), rep("Human Damage", nrow(hum.count)), rep("Chemical Exposure", nrow(chem.count)), rep("Pollution", nrow(pol.count)), rep("Fishing", nrow(fish.count)))),
                               Year = c(AI.count[,1], hum.count[,1], chem.count[,1], pol.count[,1], fish.count[,1]),
                               Count = c(AI.count[,2], hum.count[,2], chem.count[,2], pol.count[,2], fish.count[,2]))
aiterm_count <- aggregate(Count ~ Year, data = ai.drive.year.df, sum)
aiterm_count$Category <- print("Anthropogenic Impacts")

symb.drive.year.df <- data.frame(Term = as.factor(c(rep("Pathogen", nrow(virus.count)), rep("Microbiome", nrow(micro.count)), rep("Symbiosis", nrow(sym.count)))),
                                 Year = c(virus.count[,1], micro.count[,1], sym.count[,1]),
                                 Count = c(virus.count[,2], micro.count[,2], sym.count[,2]))
symbterm_count <- aggregate(Count ~ Year, data = symb.drive.year.df, sum)
symbterm_count$Category <- print("Microbe")

other.drive.year.df <- data.frame(Term = as.factor(c(rep("Range Shift", nrow(hab.count)), rep("Predation", nrow(pred.count)), rep("Wave Damage", nrow(wave.count)))),
                                  Year = c(hab.count[,1], pred.count[,1], wave.count[,1]),
                                  Count = c(hab.count[,2], pred.count[,2], wave.count[,2]))
otherterm_count <- aggregate(Count ~ Year, data = other.drive.year.df, sum)
otherterm_count$Category <- print("Other")

alldrive.df <- data.frame(Term = c(ccterm_count$Category, aiterm_count$Category, symbterm_count$Category, otherterm_count$Category),
                          Year = c(ccterm_count$Year, aiterm_count$Year, symbterm_count$Year, otherterm_count$Year),
                          Count = c(ccterm_count$Count, aiterm_count$Count, symbterm_count$Count, otherterm_count$Count))

drive.time.plot2 <- ggplot(alldrive.df) +
  geom_line(aes(x = Year, y = Count, color = Term)) +
  scale_colour_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#CC79A7")) +
  ylab("Number of Articles") +
  theme(text = element_text(size = 12), legend.title = "Term") +
  theme_bw()

```

### Figure S1C - outcome
```{r}
bleach.time <- subset(systmap, Bleaching == "Yes")
bleach.count <- count(bleach.time$Year)

dis.time <- subset(systmap, Disease == "Yes")
dis.count <- count(dis.time$Year)

bioc.time <- subset(systmap, Biodiversity_Reef == "Yes")
bioc.count <- count(bioc.time$Year)

biom.time <- subset(systmap, Biodiversity_Microbes == "Yes")
biom.count <- count(biom.time$Year)

RTA.time <- subset(systmap, Resilience == "Yes")
RTA.count <- count(RTA.time$Year)

RRR.time <- subset(systmap, Recovery == "Yes")
RRR.count <- count(RRR.time$Year)

cover.time <- subset(systmap, Mortality == "Yes")
cover.count <- count(cover.time$Year)


out.year.df <- data.frame(Term = c(rep("Bleaching", nrow(bleach.count)), rep("Disease", nrow(dis.count)), rep("Reef Biodiversity", nrow(bioc.count)), rep("Microbe Biodiversity", nrow(biom.count)), rep("Resilience", nrow(RTA.count)), rep("Recovery", nrow(RRR.count)), rep("Mortality", nrow(cover.count))),
                      Year = c(bleach.count[,1], dis.count[,1], bioc.count[,1], biom.count[,1], RTA.count[,1], RRR.count[,1], cover.count[,1]),
                      Count = c(bleach.count[,2], dis.count[,2], bioc.count[,2], biom.count[,2], RTA.count[,2], RRR.count[,2], cover.count[,2]))


out.time.plot <- ggplot(out.year.df) +
  geom_line(aes(x = Year, y = Count, color = Term)) +
  scale_colour_manual(values = c("#000000", "#E69F00", "#56B4E9", "#009E73", "#D55E00", "#CC79A7", "#999999")) +
  #scale_color_viridis(discrete = TRUE) +
  ylab("Number of Articles") +
  theme(text = element_text(size = 12)) +
  theme_bw()
```

### Combine plots
```{r}
( met.time.plot / drive.time.plot2 / out.time.plot) + 
  plot_layout(widths = c(1, 1, 1), heights = c(0.8, 1, 0.8)) +
  plot_annotation(tag_levels = "A")
```


## Figure S2 - systematic-like only terms
### Figure S2A - systematic-like only methodology
```{r}
comp.syst <- count(systematic$Quantitative)
percent_comp.syst <- comp.syst %>% mutate(Synthesis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method.syst <- count(systematic$Methodological)
percent_method.syst <- method.syst %>% mutate(Method = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

remote.syst <- count(systematic$Remote_Sensing)
percent_remote.syst <- remote.syst %>% mutate(Sensing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

mon.syst <- count(systematic$Monitoring)
percent_mon.syst <- mon.syst %>% mutate(Monitor = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

eco.syst <- count(systematic$Conservation)
percent_eco.syst <- eco.syst %>% mutate(Conservation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

method_info.syst <- data.frame(Term = c("Quantitative", "Methodological", "Remote Sensing", "Monitoring", "Conservation"),
                          Yes = c(percent_comp.syst$percent[2], percent_method.syst$percent[2], percent_remote.syst$percent[2], percent_mon.syst$percent[2], percent_eco.syst$percent[2]),
                          No = c(percent_comp.syst$percent[1], percent_method.syst$percent[1], percent_remote.syst$percent[1], percent_mon.syst$percent[1], percent_eco.syst$percent[1]))

method_info.syst <- method_info.syst %>% arrange(Yes)

for (i in 1:nrow(method_info.syst)) {
  method_info.syst$Yes[i] <- ifelse(is.na(method_info.syst$Yes[i]) == TRUE, 0, method_info.syst$Yes[i])
}

method_info.syst$Term <- reorder(method_info.syst$Term, method_info.syst$Yes)

met_plot.syst <- ggplot(method_info.syst) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "black") +
  scale_x_discrete(limits = levels(method_info.syst$Term)) +
  coord_flip(ylim = c(0, 100)) +
  theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5) + 
  ylab("") + xlab("") +
  theme_bw()

# method_long.syst <- method_info.syst %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# method_long.syst$Term <- as.factor(method_long.syst$Term)
# method_long.syst$Yes_Percent <- ifelse(method_long.syst$Presence == "Yes",  print(method_long.syst$Percent), print(NA))
# 
# met_plot.syst <- ggplot(data = method_long.syst, x = reorder(Term, -Percent), y = Percent) +
#   geom_col(mapping = aes(x = reorder(Term, -Percent), y = Percent, fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = c("Remote Sensing", "Methodological", "Monitoring", "Quantitative", "Conservation")) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.75) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "black"), name = "Presence:") +
#   theme(legend.position = "right", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("") + xlab("")
```

### Figure S2B - systematic-like only drivers
```{r}
virus.syst <- count(systematic$Pathogen)
percent_path.syst <- virus.syst %>% mutate(Pathogen = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

micro.syst <- count(systematic$Microbiome)
percent_micro.syst <- micro.syst %>% mutate(Microbiome = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sym.syst <- count(systematic$Symbiosis)
percent_sym.syst <- sym.syst %>% mutate(Symbiosis = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cc.syst <- count(systematic$Climate_Change)
percent_cc.syst <- cc.syst %>% mutate(Climate = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

AI.syst <- count(systematic$Anthropogenic)
percent_AI.syst <- AI.syst %>% mutate(Anthropogenic = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hab.syst <- count(systematic$Habitat_Range)
percent_hab.syst <- hab.syst %>% mutate(Habitat = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

sst.syst <- count(systematic$SST)
percent_sst.syst <- sst.syst %>% mutate(SST = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

uv.syst <- count(systematic$UV)
percent_uv.syst <- uv.syst %>% mutate(UV = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pred.syst <- count(systematic$Predation)
percent_pred.syst <- pred.syst %>% mutate(Predation = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

acid.syst <- count(systematic$Acidification)
percent_acid.syst <- acid.syst %>% mutate(Acidification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

calc.syst <- count(systematic$Calcification)
percent_calc.syst <- calc.syst %>% mutate(Calcification = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

erod.syst <- count(systematic$Bioerosion)
percent_ero.syst <- erod.syst %>% mutate(Bioerosion = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

wave.syst <- count(systematic$Damage_Wave)
percent_wave.syst <- wave.syst %>% mutate(Wave_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

hum.syst <- count(systematic$Damage_Human)
percent_hum.syst <- hum.syst %>% mutate(Human_Damage = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

chem.syst <- count(systematic$Chemical)
percent_chem.syst <- chem.syst %>% mutate(Chemical = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

pol.syst <- count(systematic$Pollution)
percent_pol.syst <- pol.syst %>% mutate(Pollution = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

fish.syst <- count(systematic$Fishing)
percent_fish.syst <- fish.syst %>% mutate(Fishing = x, percent = (signif((freq/sum(freq))*100, digits = 3)))


drivers.syst <- data.frame(Term = c("Pathogen", "Microbiome", "Symbiosis", "Climate Change", "Anthropogenic Impacts", "Range Shift", "Increase SST", "UV Radiation", "Predation", "Acidification", "Calcification", "Bioerosion", "Wave Damage", "Human Damage", "Chemical Exposure", "Pollution", "Fishing"),
                          Yes = c(percent_path.syst$percent[2], percent_micro.syst$percent[2], percent_sym.syst$percent[2], percent_cc.syst$percent[2], percent_AI.syst$percent[2], percent_hab.syst$percent[2], percent_sst.syst$percent[2], percent_uv.syst$percent[2], percent_pred.syst$percent[2], percent_acid.syst$percent[2], percent_calc.syst$percent[2], percent_ero.syst$percent[2], percent_wave.syst$percent[2], percent_hum.syst$percent[2], percent_chem.syst$percent[2], percent_pol.syst$percent[2], percent_fish.syst$percent[2]),
                          No = c(percent_path.syst$percent[1], percent_micro.syst$percent[1], percent_sym.syst$percent[1], percent_cc.syst$percent[1], percent_AI.syst$percent[1], percent_hab.syst$percent[1], percent_sst.syst$percent[1], percent_uv.syst$percent[1], percent_pred.syst$percent[1], percent_acid.syst$percent[1], percent_calc.syst$percent[1], percent_ero.syst$percent[1], percent_wave.syst$percent[1], percent_hum.syst$percent[1], percent_chem.syst$percent[1], percent_pol.syst$percent[1], percent_fish.syst$percent[1]))

drivers.syst <- drivers.syst %>% arrange(-Yes)
term.list.syst <- list(drivers.syst$Term)

for (i in 1:nrow(drivers.syst)) {
  drivers.syst$Yes[i] <- ifelse(is.na(drivers.syst$Yes[i]) == TRUE, 0, drivers.syst$Yes[i])
}

drivers.syst$Term <- reorder(drivers.syst$Term, drivers.syst$Yes)

driver_plot.syst <- ggplot(drivers.syst) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "violetred4") +
  scale_x_discrete(limits = levels(drivers.syst$Term)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("") + xlab("Extracted Terms - Systematic-Like Papers Only") +
  theme_bw()

# drivers_long.syst <- drivers.syst %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# drivers_long.syst$Term <- as.factor(drivers_long.syst$Term)
# drivers_long.syst$Yes_Percent <- ifelse(drivers_long.syst$Presence == "Yes",  print(drivers_long.syst$Percent), print(NA))
# 
# driver_plot.syst <- ggplot(data = drivers_long.syst, aes(x = reorder(Term, -Percent), y = Percent)) +
#   geom_col(mapping = aes(fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = rev(term.list.syst[[1]])) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "violetred4"), name = "Presence:") +
#   theme(legend.position = "right", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("") + xlab("Extracted Terms - Systematic-Like Papers Only")
```

### Figure S2C - systematic-like only outcome
```{r}
bleach.syst <- count(systematic$Bleaching)
percent_bleach.syst <- bleach.syst %>% mutate(Bleaching = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

dis.syst <- count(systematic$Disease)
percent_dis.syst <- dis.syst %>% mutate(Disease = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

bioc.syst <- count(systematic$Biodiversity_Reef)
percent_bioc.syst <- bioc.syst %>% mutate(Coral_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

biom.syst <- count(systematic$Biodiversity_Microbes)
percent_biom.syst <- biom.syst %>% mutate(Micro_Biodiv = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RTA.syst <- count(systematic$Resilience)
percent_RTA.syst <- RTA.syst %>% mutate(Resilience = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

RRR.syst <- count(systematic$Recovery)
percent_RRR.syst <- RRR.syst %>% mutate(Recovery = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

cover.syst <- count(systematic$Mortality)
percent_cover.syst <- cover.syst %>% mutate(Cover = x, percent = (signif((freq/sum(freq))*100, digits = 3)))

outcomes.syst <- data.frame(Term = c("Bleaching", "Disease", "Reef Biodiversity", "Microbe Biodiversity", "Resilience", "Recovery", "Mortality"),
                          Yes = c(percent_bleach.syst$percent[2], percent_dis.syst$percent[2], percent_bioc.syst$percent[2], percent_biom.syst$percent[2], percent_RTA.syst$percent[2], percent_RRR.syst$percent[2], percent_cover.syst$percent[2]),
                          No = c(percent_bleach.syst$percent[1], percent_dis.syst$percent[1], percent_bioc.syst$percent[1], percent_biom.syst$percent[1], percent_RTA.syst$percent[1], percent_RRR.syst$percent[1], percent_cover.syst$percent[1]))

outcomes.syst <- outcomes.syst %>% arrange(-Yes)

outcomes.syst$Term <- as.factor(outcomes.syst$Term)
outcomes.syst$Term <- reorder(outcomes.syst$Term, outcomes.syst$Yes)

out_plot.syst <- ggplot(outcomes.syst) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "lightseagreen") +
  scale_x_discrete(limits = levels(outcomes.syst$Term)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.5) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("") +
  theme_bw()

# outcome_long.syst <- outcomes.syst %>% 
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# outcome_long.syst$Term <- as.factor(outcome_long.syst$Term)
# outcome_long.syst$Term <- fct_reorder(outcome_long.syst$Term, outcome_long.syst$Percent, max, .desc = TRUE)
# outcome_long.syst$Yes_Percent <- ifelse(outcome_long.syst$Presence == "Yes",  print(outcome_long.syst$Percent), print(NA))
# 
# out_plot.syst <- ggplot(data = outcome_long.syst, aes(x = reorder(Term, -Percent), y = Percent)) +
#   geom_col(mapping = aes(fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = c("Disease", "Microbe Biodiversity", "Reef Biodiversity", "Recovery", "Bleaching", "Mortality", "Resilience")) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.5) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "lightseagreen"), name = "Presence:") +
#   theme(legend.position = "right", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("Percent (%)") + xlab("")
```

### Combine plots
```{r}
(met_plot.syst / driver_plot.syst / out_plot.syst) +
  plot_layout(widths = c(1, 1, 1), heights = c(0.5, 1, 0.5)) + 
  plot_annotation(tag_levels = "A")
```


## Figure S3 - keyword cloud
```{r}
# isolate keywords
keyword.data <- bib$DE

# create a corpus for cleaning text
corpkeys <- Corpus(VectorSource(keyword.data))

# cleaning
keywords <- corpkeys %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
keywords <- tm_map(keywords, content_transformer(tolower))
keywords <- tm_map(keywords, removeWords, stopwords("english"))

# create dataframe of words and freq
words <- TermDocumentMatrix(keywords) 
word.matrix <- as.matrix(words) 
wordcounts <- sort(rowSums(word.matrix),decreasing=TRUE) 
keyword.df <- data.frame(keywords = names(wordcounts),freq = wordcounts)

# word cloud
set.seed(1234) # for reproducibility
wordcloud(words = keyword.df$keywords, freq = keyword.df$freq, min.freq = 1, max.words=200, random.order=FALSE, rot.per=0.35, colors=brewer.pal(8, "Dark2"))
```


## Figure S4 - country collaborations coloured by HDI
```{r}
CountryNetwork <- biblioNetwork(bib, analysis = "collaboration", 
                           network = "countries", sep = ";")

collab <- as.matrix(CountryNetwork)
diag(collab) <- 0 #get rid of collaboration with same country

# getting rid of lower triangle (as this is duplication of info)
collab[lower.tri(collab)] <- 0 
# change row and col names case
colnames(collab) <- str_to_title(colnames(collab))
rownames(collab) <- str_to_title(rownames(collab))
# change to "USA"
colnames(collab)[colnames(collab) == "Usa"] <- "USA"
rownames(collab)[rownames(collab) == "Usa"] <- "USA"
# change to "UK"
colnames(collab)[colnames(collab) == "United Kingdom"] <- "UK"
rownames(collab)[rownames(collab) == "United Kingdom"] <- "UK"
# change to "NewZealand"
colnames(collab)[colnames(collab) == "New Zealand"] <- "NewZealand"
rownames(collab)[rownames(collab) == "New Zealand"] <- "NewZealand"
# change to "SaudiArabia"
colnames(collab)[colnames(collab) == "Saudi Arabia"] <- "SaudiArabia"
rownames(collab)[rownames(collab) == "Saudi Arabia"] <- "SaudiArabia"
# change to "SouthAfrica"
colnames(collab)[colnames(collab) == "South Africa"] <- "SouthAfrica"
rownames(collab)[rownames(collab) == "South Africa"] <- "SouthAfrica"
# change to "CostaRica"
colnames(collab)[colnames(collab) == "Costa Rica"] <- "CostaRica"
rownames(collab)[rownames(collab) == "Costa Rica"] <- "CostaRica"
# change to "DominicanRepublic"
colnames(collab)[colnames(collab) == "Dominican Republic"] <- "DominicanRepublic"
rownames(collab)[rownames(collab) == "Dominican Republic"] <- "DominicanRepublic"
# change to "HongKong"
colnames(collab)[colnames(collab) == "Hong Kong"] <- "HongKong"
rownames(collab)[rownames(collab) == "Hong Kong"] <- "HongKong"
# change to "SriLanka"
colnames(collab)[colnames(collab) == "Sri Lanka"] <- "SriLanka"
rownames(collab)[rownames(collab) == "Sri Lanka"] <- "SriLanka"

circos.clear()

# reference: https://hdr.undp.org/system/files/documents/global-report-document/hdr2021-22pdf_1.pdf (above medium HDI = developed)
colours.develop <- c(USA = "#440154FF", Australia = "#440154FF", UK = "#440154FF", China = "#2D708EFF", Israel = "#440154FF", Germany = "#440154FF", Brazil = "#2D708EFF", Canada = "#440154FF", France = "#440154FF", Mexico = "#2D708EFF", Netherlands = "#440154FF", Sweden = "#440154FF", SaudiArabia = "#440154FF", SouthAfrica = "#2D708EFF", Colombia = "#2D708EFF", Italy = "#440154FF", Kenya = "#73D055FF", Philippines = "#73D055FF", Portugal = "#440154FF", India = "#73D055FF", Fiji = "#2D708EFF", Japan = "#440154FF", Malaysia = "#440154FF", Monaco = "#FDE725FF", NewZealand = "#440154FF", Denmark = "#440154FF", Norway = "#440154FF", Panama = "#440154FF", Singapore = "#440154FF", Switzerland = "#440154FF", Thailand = "#440154FF", Venezuela = "#73D055FF", Austria = "#440154FF", Belgium = "#440154FF", CostaRica = "#440154FF", DominicanRepublic = "#2D708EFF", Guam = "#440154FF", Indonesia = "#2D708EFF", Samoa = "#2D708EFF", Spain = "#440154FF", Tanzania = "#73D055FF", Mauritius = "#440154FF", Greece = "#440154FF", Egypt = "#2D708EFF", Ireland = "#440154FF", HongKong = "#440154FF")


C.network <- chordDiagram(collab,
                          annotationTrack = "grid",
                          preAllocateTracks = 1,
                          grid.col = colours.develop)
C.network <- select(C.network, c(rn, cn, value1, col))
C.network <- C.network[order(-C.network$value1),]
C.network <- C.network %>% group_by(value1)

chordDiagram(C.network,
             annotationTrack = "grid",
             preAllocateTracks = 1,
             grid.col = colours.develop)

circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
  circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
}, bg.border = NA)

```


## Figure S5 - full author collaboration plot
```{r}
AuthorNetwork <- biblioNetwork(bib, analysis = "collaboration", network = "authors", sep = ";")

networkPlot(AuthorNetwork, Title = NA, cluster = "infomap", label = FALSE, remove.isolates = TRUE)
```


## Figure S6
### Figure S6A - Systematic Papers through time
```{r}
syst.timeline <- ggplot(syst.bib) +
  geom_bar(aes(x = PY)) +
  labs(x = "Year", y = "Number of Systematic-Like Articles Published") +
  xlim(2002,2022) +
  ylim(0,8) +
  theme_bw()
```

### Figure S6B - Systematic Papers that publish data through time
```{r}
have.data <- subset(CA_reduced, Q14 == "Yes")

data.matches <- sapply(syst.bib$TI, function(x) {
  max(agrep(x, have.data$TITLE, max.distance = 0.3, ignore.case = TRUE))
})

data.bib <- syst.bib[!is.infinite(data.matches), ]
# missing Bove paper since this was not found in Scopus records

data.timeline <- ggplot(data.bib) +
  geom_bar(aes(x = PY)) +
  labs(x = "Year", y = "Number of Systematic-Like Articles Published with Data") +
  xlim(2002,2022) +
  ylim(0,8) +
  theme_bw()
```

### Figure S6C - Systematic Papers that publish code through time
```{r}
have.code <- subset(CA_reduced, Q15 == "Yes")

code.matches <- sapply(syst.bib$TI, function(x) {
  max(agrep(x, have.code$TITLE, max.distance = 0.3, ignore.case = TRUE))
})

code.bib <- syst.bib[!is.infinite(code.matches), ]
# missing Bove paper since this was not found in Scopus records

code.timeline <- ggplot(code.bib) +
  geom_bar(aes(x = PY)) +
  labs(x = "Year", y = "Number of Systematic-Like Articles Published with Code") +
  xlim(2002,2022) +
  ylim(0,8) +
  theme_bw()
```

### Combine Plots
```{r}
(syst.timeline / data.timeline / code.timeline) + plot_annotation(tag_levels = "A")
```


## Figure S7 - Additional Terms
```{r}
# terms were added in "Other" column of systmap data
OtherData <- subset(systmap, Other != "")
OtherTerms <- unlist(strsplit(OtherData$Other, "; "))
OtherTerms <- unlist(strsplit(OtherTerms, ", "))

# cleaning to get similar terms to match
OtherTerms <- str_to_lower(OtherTerms)
OtherTerm <- sub("s$", "", OtherTerms)
OtherTerm <- sub("herbivory", "herbivore", OtherTerm)
OtherTerm <- sub("^trophic ", "", OtherTerm)
OtherTerm <- sub("^disease ", "", OtherTerm)
OtherTerm <- sub("^larval ", "", OtherTerm)
OtherTerm <- sub("^stress ", "", OtherTerm)
OtherTerm <- sub("^spatial ", "", OtherTerm)
OtherTerm <- sub(" rise$", "", OtherTerm)
OtherTerm <- sub("photosynthesi", "photosynthesis", OtherTerm)

# how many
count_terms <- count(OtherTerm)
AddedTerms <- subset(count_terms, count_terms$freq > 5)

# percentage 
count_terms <- count_terms %>% mutate(as.percent = signif((freq/335)*100, digits = 3))

# ggplot(count_terms) +
#   geom_bar(aes())

NewTerm_df <- data.frame(Term = count_terms$x,
                         Yes = count_terms$as.percent,
                         No = (100 - count_terms$as.percent))

NewTerm_df <- NewTerm_df %>% arrange(-Yes)
# Newterm.list <- list(NewTerm_df$Term)
# 
# NewTerm_long <- NewTerm_df %>%
#   pivot_longer(
#     cols = !Term,
#     names_to = "Presence",
#     values_to = "Percent"
# )
# 
# NewTerm_long$Term <- as.factor(NewTerm_long$Term)
# NewTerm_long$Yes_Percent <- ifelse(NewTerm_long$Presence == "Yes",  print(NewTerm_long$Percent), print(NA))


NewTerm_df$Term <- as.factor(NewTerm_df$Term)
NewTerm_df$Term <- reorder(NewTerm_df$Term, NewTerm_df$Yes)

figS7B <- ggplot(NewTerm_df) +
  geom_col(aes(x = Term, y = Yes), width = 0.7, color = "black", fill = "deepskyblue") +
  scale_x_discrete(limits = levels(NewTerm_df$Term)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Term, y = Yes, label = Yes), hjust = -0.23) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Additional Terms") +
  theme_bw()


# # plot
# ggplot(data = NewTerm_long, aes(x = reorder(Term, -Percent), y = Percent)) +
#   geom_col(mapping = aes(fill = Presence), width = 0.7,
#            position = "stack", color = "black") +
#   scale_x_discrete(limits = rev(Newterm.list[[1]])) +
#   coord_flip(ylim = c(0, 100)) +
#   geom_text(aes(x = Term, y = Percent, label = Yes_Percent), hjust = -0.05) +
#   guides(fill = guide_legend(reverse = TRUE)) +
#   scale_fill_manual(values = c("gainsboro", "deepskyblue"), name = "Presence:") +
#   theme(legend.position = "bottom", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
#   ylab("Percent (%)") + xlab("Additional Terms")
```

```{r}
# separate systmap$Other into different cols using stringr
OtherData <- subset(systmap, Other != "")
# define max number of terms in Other
max_terms <- 3
# Split Other into different columns
split_terms <- str_split_fixed(OtherData$Other, "; ", max_terms)
# Create column names for the split items
column_names <- paste0("term", 1:max_terms)
# Assign the split items to the dataset
OtherData[column_names] <- split_terms

biochem_response <- c('bioaccumulation', 'chlorophyll', 'functional gene', 'genetic', 'genetic response', 'immune response', 'metabolism', 'photobiology', 'photosynthesis', 'pigmentation', 'productivity', 'respiration', 'signaling molecules')
OtherData$biochem <- NA

biophys_drivers <- c('accretion', 'depth', 'salinity', 'sea level', 'water flow', 'wind')
OtherData$biophys <- NA

reef_process <- c('aging', 'algae', 'competition', 'connectivity', 'growth', 'herbivory', 'hybridization', 'interaction', 'phase shift', 'recruitment', 'reproduction', 'succession', 'transmission')
OtherData$reefp <- NA

other <- c('aquaculture', 'captivity', 'evolution', 'feedback loop', 'life history', 'marine protected area', 'physiology', 'science communication', 'soft coral')
OtherData$other_term <- NA

term_df <- data.frame(biochem = biochem_response,
                      biophys = c(biophys_drivers, rep("", times = 7)),
                      reefp = reef_process,
                      other_term = c(other, rep("", times = 4))
)

# match OtherData$Other to term_df
term.matcher <- function(data, term)
  sapply(data, function(x) {
  max(agrep(x, term, max.distance = 0.25, ignore.case = TRUE))
})

for (i in 1:nrow(OtherData)) {
  term.matches1 <- term.matcher(OtherData$term1[i], term_df$biochem)
  OtherData$biochem[i] <- ifelse(!is.infinite(term.matches1), 1, 0)
  term.matches2 <- term.matcher(!is.na(OtherData$term2[i]), term_df$biochem)
  OtherData$biochem[i] <- ifelse(!is.infinite(term.matches2), 1, ifelse(OtherData$biochem[i] == 1, 1, 0))
  term.matches3 <- term.matcher(!is.na(OtherData$term3[i]), term_df$biochem)
  OtherData$biochem[i] <- ifelse(!is.infinite(term.matches3), 1, ifelse(OtherData$biochem[i] == 1, 1, 0))
  term.matches1 <- term.matcher(OtherData$term1[i], term_df$biophys)
  OtherData$biophys[i] <- ifelse(!is.infinite(term.matches1), 1, 0)
  term.matches2 <- term.matcher(!is.na(OtherData$term2[i]), term_df$biophys)
  OtherData$biophys[i] <- ifelse(!is.infinite(term.matches2), 1, ifelse(OtherData$biophys[i] == 1, 1, 0))
  term.matches3 <- term.matcher(!is.na(OtherData$term3[i]), term_df$biophys)
  OtherData$biophys[i] <- ifelse(!is.infinite(term.matches3), 1, ifelse(OtherData$biophys[i] == 1, 1, 0))
  term.matches1 <- term.matcher(OtherData$term1[i], term_df$reefp)
  OtherData$reefp[i] <- ifelse(!is.infinite(term.matches1), 1, 0)
  term.matches2 <- term.matcher(!is.na(OtherData$term2[i]), term_df$reefp)
  OtherData$reefp[i] <- ifelse(!is.infinite(term.matches2), 1, ifelse(OtherData$reefp[i] == 1, 1, 0))
  term.matches3 <- term.matcher(!is.na(OtherData$term3[i]), term_df$reefp)
  OtherData$reefp[i] <- ifelse(!is.infinite(term.matches3), 1, ifelse(OtherData$reefp[i] == 1, 1, 0))
  term.matches1 <- term.matcher(OtherData$term1[i], term_df$other_term)
  OtherData$other_term[i] <- ifelse(!is.infinite(term.matches1), 1, 0)
  term.matches2 <- term.matcher(!is.na(OtherData$term2[i]), term_df$other_term)
  OtherData$other_term[i] <- ifelse(!is.infinite(term.matches2), 1, ifelse(OtherData$other_term[i] == 1, 1, 0))
  term.matches3 <- term.matcher(!is.na(OtherData$term3[i]), term_df$other_term)
  OtherData$other_term[i] <- ifelse(!is.infinite(term.matches3), 1, ifelse(OtherData$other_term[i] == 1, 1, 0))
}

count_groups <- data.frame(Category = c("Biochemical Response", "Biophysical Drivers", "Reef Processes", "Other"),
                           Count = c(sum(OtherData$biochem), sum(OtherData$biophys), sum(OtherData$reefp), sum(OtherData$other_term)))

count_groups <- count_groups %>% mutate(as.percent = signif((Count/335)*100, digits = 3))

Group_df <- data.frame(Category = count_groups$Category,
                         Yes = count_groups$as.percent,
                         No = (100 - count_groups$as.percent))

Group_df <- Group_df %>% arrange(-Yes)


Group_df$Category <- as.factor(Group_df$Category)
Group_df$Category <- reorder(Group_df$Category, Group_df$Yes)

figS7A <- ggplot(Group_df) +
  geom_col(aes(x = Category, y = Yes), width = 0.7, color = "black", fill = "deepskyblue") +
  scale_x_discrete(limits = levels(Group_df$Category)) +
  coord_flip(ylim = c(0, 100)) +
  geom_text(aes(x = Category, y = Yes, label = Yes), hjust = -0.23) +
  theme(legend.position = "none", panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(), axis.title = element_text(size = 14), axis.text = element_text(size = 12)) + 
  ylab("Percent (%)") + xlab("Categories") +
  theme_bw()

```

```{r}
(figS7A / figS7B)  + 
  plot_layout(widths = c(1, 1), heights = c(0.25, 1))+
  plot_annotation(tag_levels = "A")
```


