Thesis Code

Install & Load the Required Packages

## install packages
install.packages("tidyverse", repos="http://cran.us.r-project.org") ## data wrangling
install.packages("maps", repos="http://cran.us.r-project.org") ## spatial database
install.packages("mapproj", repos="http://cran.us.r-project.org") ## coordinate projection
install.packages("scico", repos="http://cran.us.r-project.org") ## color-blind accessible palettes
install.packages("prettydoc", repos="http://cran.us.r-project.org") ## rmd themes

## load packages
library(tidyverse)
library(maps)
library(mapproj)
library(scico)
library(prettydoc)

Import PubMed Data

## import .csv as tbl
thesisdata <- read.csv("/Users/stockten/Documents/Thesis/PubMed Data/thesis.csv", na.strings = c("", "NA")) %>% ## replaces blank cell w/ NA
  as_tibble()

## create tbl w/ select spatial data
spatialdata <- thesisdata %>% 
  select(Event.ID, Location, City, Country, Continent, Latitude, Longitude, Length.of.Event) %>% 
  drop_na(Event.ID) ## only include unique event reports

Create Publication Bias Bar Graphs & Choropleths

## create tbl for days observed by country
dayscountry <- spatialdata %>% 
  group_by(Country) %>%
  summarize(totaldays = sum(`Length.of.Event`), .groups = "drop") %>%
  arrange(desc(totaldays)) ## order the tbl based on # of reported events

## create tbl for lethal event number reported by country
freqcountry <- spatialdata %>% 
  count(Country) %>%
  arrange(desc(n)) ## order the tbl based on # of reported events

## create bar plot for days observed by country
daysbycountry_bar <- ggplot(dayscountry, aes(x = reorder(Country, desc(totaldays)), y = totaldays)) +
  geom_col() +
  labs(text = element_text(family = "Times New Roman", size = 8)) +
  xlab("Country") +
  ylab("# of Days Observed")
daysbycountry_bar

## create bar plot for reported events by country
countbycountry_bar <- ggplot(freqcountry, aes(x = reorder(Country, desc(n)), y = n)) +
  geom_col()+
  labs(text = element_text(family = "Times New Roman", size = 8)) +
  xlab("Country") +
  ylab("# of Lethal Heatwaves Reported")
countbycountry_bar

## generate map of world; join w/ days observed & events reported by country
worldmap <- map_data("world") %>%
  as_tibble() %>% 
  left_join(freqcountry, by = c("region" = "Country")) %>%
  left_join(dayscountry, by = c("region" = "Country"))

## create choropleth of days observed by country
daysbycountry_map <- ggplot(worldmap, aes(x = long, y = lat, group = region)) +
  geom_map(aes(map_id = region), map = worldmap) + 
  coord_map("aitoff") + ## apply map projection 
  geom_polygon(aes(group = group, fill = totaldays)) + 
  scico::scale_fill_scico(palette = "grayC", na.value = "white") +
  labs(fill = "# of Heatwave Days Observed") +
  theme(text = element_text(family = "Times New Roman", size = 12),
        legend.text = element_text(size = 8, angle = 45, vjust = 0.5),
        legend.position = "bottom",
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank()) 
daysbycountry_map

## create choropleth of reported events by country
countbycountry_map <- ggplot(worldmap, aes(x = long, y = lat, group = region)) +
  geom_map(aes(map_id = region), map = worldmap) + 
  coord_map("aitoff") + ## apply map projection 
  geom_polygon(aes(group = group, fill = n)) + 
  scico::scale_fill_scico(palette = "grayC", na.value = "white") +
  labs(fill = "# of Lethal Heatwaves Reported on PubMed from 1975-2011") +
  theme(text = element_text(family = "Times New Roman", size = 12),
        legend.text = element_text(size = 8, angle = 45, vjust = 0.5),
        legend.position = "bottom",
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank()) 
countbycountry_map