avian-mamm-visualization.Rmd

---
title: 'AVIAN-MAMM: A Global Open Access Dataset of Reported Avian Flu Events in Mammals
  - Visualization'
author: 'Francesco Branda'
date: ''
output:
  pdf_document: default
---

```{r setup, include=FALSE}
# Prepare the R environment: install packages

knitr::opts_chunk$set(echo = TRUE, fig.align   = 'center')

# Define requested packages
my_packages <- c("dplyr", "ggplot2", "viridis", "ggalluvial","ggfittext","randomcoloR", "reshape2", "rworldmap", "stringr", "webr" , "knitr", "kableExtra","ggfittext" , "visNetwork", "scales", "treemap")  

# Extract not installed packages
not_installed <- my_packages[!(my_packages %in% installed.packages()[ , "Package"])]    

# Install not installed packages
if(length(not_installed)) install.packages(not_installed) 

devtools::install_github("sctyner/geomnet")

# Load libraries
library(dplyr)
library(ggplot2)
library(viridis)
library(ggalluvial)
library(ggfittext)
library(randomcoloR)
library(reshape2)
library(rworldmap)
library(stringr)
require(webr)
library(knitr)
library(kableExtra)
library(ggfittext)
library(visNetwork)
library(geomnet)
library(scales)
library(treemap)
```

# Import the data
```{r import data}
avian_flu_df <- read.csv("latest-outbreaks-mammals.csv", encoding="UTF-8") 
```

## How many avian flu events in animals have been published ? 
The figure includes follow-up history reports of outbreaks. Date of publication is used.
```{r number_reports, fig.height=5, fig.width=5, echo=FALSE, message=FALSE, warning=FALSE}
# Cumulative number of avian flu events. Date of publication is used.
# Note: this does not represent number of outbreaks but number of events, i.e. follow-up history reports are included

nb_report_day <- avian_flu_df %>%
  mutate (report_date = as.Date (report_date, format = "%Y-%m-%d")) %>%
  group_by(report_date) %>%
  dplyr::summarize(n()) %>% 
  arrange(report_date) %>% 
  dplyr::rename(number_reports = "n()") %>% 
  mutate(cum_sum=cumsum(number_reports))

# Make  the curve of the cumulative number of events
plot_day_cum <- ggplot(nb_report_day , aes(x=report_date, y=cum_sum, group = 1)) + 
  geom_line()+
  geom_point()+ # remove this line if you want to remove the dots on the plot line (or add "#" in front of this line)
  scale_x_date(limits= as.Date(c(min(nb_report_day$report_date)-7, max(nb_report_day$report_date)+7)), date_breaks = "3 week", date_labels =  "%Y-%m-%d", expand = c(0, 0))+
  xlab("Date of publication") +
  ylab("Cumulative number of events")+
  theme_bw()+
  theme(axis.text.x = element_text(angle=90,vjust = 0.5, hjust=.50))+
  theme (axis.title.x = element_text(vjust= -1),
         axis.title.y = element_text(vjust= 3),
         text = element_text(size = 10))
plot_day_cum
```
\newpage

## Where did avian flu events in animals occur? 

### For each country, how many species have been reported as infected with or exposed to avian flu?

```{r map_species, fig.height=10, fig.width=21, fig.fullwidth=TRUE, echo=FALSE, message=FALSE, warning=FALSE}
# Get the world map
world_map <- map_data("world")

# Check that countries names are similar between the two datasets (world map and AVIAN-MAMM data)  
countries_world <- sort(unique(world_map$region))
countries_avian_flu <- sort(unique(avian_flu_df$country_name))

intersect <- intersect(countries_world,countries_avian_flu)
#length(intersect) == length(countries_avian_flu) # if FALSE some country names have to be correct

# Find values to correct
#setdiff(countries_avian_flu, countries_world)

# Format the data to plot
# Correct countries names to match the world map data
avian_flu_to_plot_host <- avian_flu_df %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "United States of America", "USA")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "United Kingdom", "UK")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "Republic of Korea", "South Korea")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "Russian Federation", "Russia")) %>%
  dplyr::group_by(country_name) %>% 
  dplyr::mutate(unique_species = n_distinct(species)) %>% 
  select(country_name,unique_species) %>%
  unique()


# Add species number to world map
avian_flu.to.map <- left_join( world_map, avian_flu_to_plot_host, by = c("region"="country_name"))

# Plot the map of the number of species with SARS-CoV-2 reported cases per countries
avian_flu.to.map %>%
  filter (region != "Antarctica") %>%
  ggplot(aes(long, lat, group = group))+
  geom_polygon(aes(fill = unique_species), color = "black") + 
  scale_fill_viridis_c(name = "Number of unique (sub)species", option = 'D', trans = "log2",  breaks = trans_breaks("log2", function(x) 2^x),na.value="grey85") +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(), axis.title = element_blank(),
        plot.margin = unit(0 * c(-1.5, -1.5, -1.5, -1.5), "lines"))+ 
  xlab ("Longitude") +
  ylab ("Latitude") +
  theme_bw()+ 
  theme(legend.position="bottom") +
  guides(fill = guide_colourbar(barwidth = 20, barheight = 2)) +
  theme(text = element_text(size = 20)) 
```
&nbsp;
&nbsp;
```{r table_countries_spec, echo=FALSE, message=FALSE, warning=FALSE}
# Print the table of the number of affected species per country
colnames(avian_flu_to_plot_host) <- c("Country", "Number of species")

kable(avian_flu_to_plot_host  , booktabs = T,
             caption = "Number of unique mammal species with reported natural avian flu infection or exposure per country.")
```
\newpage

### For each country, how many outbreaks of avian flu have been reported in mammals?
```{r map_outbreaks, fig.height=10, fig.width=21, fig.fullwidth=TRUE, echo=FALSE, message=FALSE, warning=FALSE}
# Format the data to be plotted (keep outbreaks only, aggregate per country)
avian_flu_to_plot_outbreaks <- avian_flu_df %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "United States of America", "USA")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "United Kingdom", "UK")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "Republic of Korea", "South Korea")) %>%
  dplyr::mutate(country_name = replace(country_name, country_name == "Russian Federation", "Russia")) %>%
  group_by(country_name) %>% 
  tally()

# Add outbreak number to world map
avian_flu.to.map_2 <- left_join( world_map, avian_flu_to_plot_outbreaks, by = c("region"="country_name"))

avian_flu.to.map_2 %>%
  filter (region != "Antarctica") %>% # crop the map
  ggplot(aes(long, lat, group = group))+
  geom_polygon(aes(fill = n), color = "black") + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(), axis.title = element_blank(),
        plot.margin = unit(0 * c(-1.5, -1.5, -1.5, -1.5), "lines"))+ 
  xlab ("Longitude") +
  ylab ("Latitude") +
  theme_bw()+ 
  theme(legend.position="bottom") +
  scale_fill_viridis_c(name = "Number of reported outbreaks", option = 'D', trans = "log2",  breaks = trans_breaks("log2", function(x) 2^x),na.value="grey85") +  
  guides(fill = guide_colourbar(barwidth = 20, barheight = 2)) +
  theme(text = element_text(size = 20)) 

# Export Figure 2
# Font size needs to be adapted for export in eps or tiff format
Fig2 <- avian_flu.to.map_2 %>%
  filter (region != "Antarctica") %>% # crop the map
  ggplot(aes(long, lat, group = group))+
  geom_polygon(aes(fill = n), color = "black") + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(), axis.title = element_blank(),
        plot.margin = unit(0 * c(-1.5, -1.5, -1.5, -1.5), "lines"))+ 
  xlab ("Longitude") +
  ylab ("Latitude") +
  theme_bw()+ 
  theme(legend.position="bottom") +
  scale_fill_viridis_c(name = "Number of reported outbreaks", option = 'D', trans = "log2", breaks = trans_breaks("log2", function(x) 2^x),na.value="grey85") +  
  guides(fill = guide_colourbar(barwidth = 8, barheight = 1)) +
  theme(text = element_text(size = 12)) 

ggsave("Fig2.eps", Fig2 , width = 18, height = 12, units="cm", dpi = 600) 
ggsave("Fig2.tiff", Fig2 , width = 18, height = 12, units="cm", dpi = 600) 
```  
  
```{r table_countries_outb, echo=FALSE, message=FALSE, warning=FALSE}

avian_flu_to_plot_outbreaks_2 <- avian_flu_to_plot_outbreaks %>%
  arrange(-n)

colnames(avian_flu_to_plot_outbreaks_2) <- c("Country", "Number of outbreaks")

kable(avian_flu_to_plot_outbreaks_2, booktabs = T,
             caption = "Number of reported avian flu outbreaks per country.")
```

\newpage

## How many avian flu cases (infections and/or exposures) have been reported per animal species? 

```{r species_families, echo=FALSE, message=FALSE, warning=FALSE}
table_species <- avian_flu_df %>%
  dplyr::mutate(number_cases = as.numeric(as.character(total_cases))) %>%
  dplyr::group_by(species) %>% 
  dplyr::summarise(number_cases_tot= sum(number_cases)) %>% 
  dplyr::arrange(-number_cases_tot)
```

The dataset reports `r sum(table_species$number_cases_tot)` cases in total **(we counted one case for any missing value on the number of cases)**.

```{r species_table, echo=FALSE, message=FALSE, warning=FALSE}
# Table of the data
colnames(table_species) <- c("Species", "Number_cases")
kable(table_species, format = "latex", booktabs = T,
             caption = "Number of reported avian flu cases (infections and/or exposures) per animal species.") %>%
          kable_styling(latex_options = "scale_down")
```

\newpage

```{r treemap, echo=FALSE, message=FALSE, warning=FALSE}
treemap(table_species,
            index=c("Species","Number_cases"),
            vSize="Number_cases",
            type="index",
        title="",                            
        fontsize.labels=c(25,16), 
    fontcolor.labels=c("white","yellow"),
    fontface.labels=c(2,1),  
    bg.labels=c("transparent"), 
    align.labels=list(
        c("center", "center"), 
        c("right", "bottom")),  
    overlap.labels=0.5,                      
    inflate.labels=F) 
```

\newpage

## How many avian flu cases (infections or exposures) have been reported per animal host and country?
```{r species_countries_scientific, fig.height=14, fig.width=18, echo=FALSE, message=FALSE, warning=FALSE}
# Use a heat map to visualize the data
avian_flu_df %>%
  dplyr::mutate(number_cases = as.numeric(as.character(total_cases))) %>%
  dplyr::group_by(country_name, species) %>% 
  dplyr::summarize (number_cases_tot = sum(number_cases)) %>% 
  ggplot(aes(country_name, species, fill= number_cases_tot)) + 
  geom_tile(color = "white") +
  coord_fixed() + 
  xlab("") + 
  ylab ("Animal host")+
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.5)) +
  theme (axis.title.x = element_text(vjust= -1), 
         axis.title.y = element_text(vjust= 3))+
  scale_fill_viridis_c(option = 'D', trans = "log2",  breaks = trans_breaks("log2", function(x) 2^x), name = "Number of cases", na.value = "firebrick4") + # NA values (was NS in the original dataset) are colored in firebrick
  guides(fill = guide_colourbar(barwidth = 2, barheight = 10)) +
  theme(text = element_text(size = 24)) 
```

\newpage

## How many avian flu outbreaks have been reported per animal host and country?

The following figure displays animal host names at species level.<br>  
```{r outbreaks, fig.height=14, fig.width=18,echo=FALSE, message=FALSE, warning=FALSE}
avian_flu_df %>%
  dplyr::group_by(country_name, species) %>% 
  tally() %>% # Count the number of outbreaks
  ggplot(aes(country_name, species, fill= n)) +
  geom_tile(color = "white") +
  coord_fixed() + 
  xlab("") + 
  ylab ("Animal host (species level)")+
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.5), 
        axis.title.x = element_text(vjust= -1),
        axis.title.y = element_text(vjust= 3))+
   scale_fill_viridis_c(option = 'D', trans = "log2", breaks = trans_breaks("log2", function(x) 2^x), name = "Number of outbreaks", na.value = "firebrick4")+
 guides(fill = guide_colourbar(barwidth = 2, barheight = 10)) +
  theme(text = element_text(size = 24)) 
```

\newpage

## How many animals have died directly or indirectly from avian flu?
The figure below includes deaths that are directly related to an infection as well as deaths related to culling when implemented as control measure. <br>  

In the figure below, the colloquial name of the animal host is displayed.<br>  

```{r deaths, fig.height=14, fig.width=18, echo=FALSE, message=FALSE, warning=FALSE}
 avian_flu_df %>%
  dplyr::mutate(number_deaths = as.numeric(as.character(total_deaths))) %>%
  dplyr::group_by(country_name, species) %>% 
  dplyr::summarise (number_deaths_tot= sum(number_deaths)) %>% 
  ggplot(aes(country_name, species, fill= number_deaths_tot)) + 
  geom_tile(color = "white")  +
  coord_fixed()+ 
  xlab("") + 
  ylab ("Animal host")+
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.5), 
        axis.title.x = element_text(vjust= -1),
        axis.title.y = element_text(vjust= 3))+
   scale_fill_viridis_c(option = 'D', trans = "sqrt", breaks = trans_breaks("sqrt", function(x) x^2),  name = "Number of deaths", na.value = "firebrick4") + # NA values (was NS in the data) are colored in firebrick
  guides(fill = guide_colourbar(barwidth = 2, barheight = 10)) +
  theme(text = element_text(size = 24)) 
```



\newpage

## What is the case fatality rate of avian flu in the different animal species, per country? 
The case fatality rate (CFR) per species and country is obtained by dividing the total number of reported deaths in one species by the total number of reported cases for this species in the country. The CFR depends strongly on testing and does not give information on the infection fatality rate (IFR, number of deaths divided by the total number of infected individuals) or mortality rate (MR, number of deaths divided by the total at-risk population).<br>  

```{r fatality, fig.height=14, fig.width=18,echo=FALSE, message=FALSE, warning=FALSE}
# To visualize the case fatality rate (proportion of individuals who died from avian flu among all individuals diagnosed with the disease), we plot (number of reported deaths / number of reported cases) on the heat map
avian_flu_df %>%
  dplyr::mutate(number_cases = as.numeric(as.character(total_cases))) %>%
  dplyr::mutate(number_deaths = as.numeric(as.character(total_deaths))) %>%
  dplyr::group_by(country_name, species) %>% 
  dplyr::summarise (fatality_rate= sum(number_deaths, na.rm = T)/sum(number_cases, na.rm = T)) %>% 
  ggplot(aes(country_name, species, fill= fatality_rate)) +
  geom_tile(color = "white")  +
  coord_fixed()+ 
  xlab("") + 
  ylab ("Animal host")+
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.5), 
        axis.title.x = element_text(vjust= -1),
        axis.title.y = element_text(vjust= 3))+
 scale_fill_viridis_c(option = 'D', trans = "sqrt",  breaks = trans_breaks("sqrt", function(x) x ^ 2), name = "Case fatality rate", na.value = "firebrick4") +
 guides(fill = guide_colourbar(barwidth = 2, barheight = 10)) +
  theme(text = element_text(size = 24)) 

# Export Fig. 3
Fig3 <- avian_flu_df %>%
  dplyr::mutate(number_cases = as.numeric(as.character(total_cases))) %>%
  dplyr::mutate(number_deaths = as.numeric(as.character(total_deaths))) %>%
  dplyr::group_by(country_name, species) %>% 
  dplyr::summarise (fatality_rate= sum(number_deaths, na.rm = T)/sum(number_cases, na.rm = T)) %>% 
  ggplot(aes(country_name, species, fill= fatality_rate)) +
  geom_tile(color = "white")  +
  coord_fixed()+ 
  xlab("") + 
  ylab ("Animal host") +
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.3), 
        axis.title.x = element_text(vjust= -1),
        axis.title.y = element_text(vjust= 3),
        text = element_text(size = 9)) +
  scale_fill_viridis_c(option = 'D', trans = "sqrt",  breaks = trans_breaks("sqrt", function(x) x ^ 2), name = "Case fatality rate", na.value = "firebrick4") +
  guides(fill = guide_colourbar(barwidth = 1, barheight = 6)) 

ggsave("Fig3.eps", Fig3 , width = 18, height = 14, units="cm", dpi = 600) 
ggsave("Fig3.tiff", Fig3 , width = 18, height = 14, units="cm", dpi = 600) 
```

\newpage


## Which laboratory methods are used to detect infection with or exposure to avian flu in animals? 
```{r tests, fig.height=14, fig.width=18, echo=FALSE, message=FALSE, warning=FALSE}
# Type of tests used to detect avian flu infection in animals
test_df <- avian_flu_df %>%
  select (species, test_name_1, test_name_2, test_name_3) 

# Convert to long format
test_df_melt <- melt (test_df, id.vars=c("species"))
# remove NA value as well as records where test is not specified ("NS")
test_df_melt <- test_df_melt %>%
  filter(!is.na(value)) %>%
  filter (value != "NA")

# Heat map of the type of test used per species
test_df_melt %>%
  ggplot(aes(species, value)) +
  geom_tile(color = "white", fill = "black")  +
  coord_fixed()+ 
  xlab("") + 
  ylab ("")+
  theme(axis.text.x = element_text(angle=90, hjust=1), 
        axis.title.x = element_text(vjust= -1),
        axis.title.y = element_text(vjust= 3),
        text = element_text(size = 24)) 
```
\newpage


## What are the implemented control measures and possible outcomes after detection of avian flu infection or exposure in the different mammal species?
The computed network is saved as an .html file (network.html) in the working directory. <br>  
The network shows, for each mammal host (blue nodes), which control measures (yellow nodes) have been implemented in response to avian flu outbreaks (e.g. treatment, culling, exclusion zone) and what was the outcome (red nodes) for the mammals (e.g. death). <br>  
The size of each node in the network corresponds to the number of times the value appears in the dataset. <br>  
```{r network, echo=FALSE, message=FALSE, warning=FALSE}
# Nodes
# We create nodes weighted with the number of time the animal host, control measure or outcome occur in the dataset
nodes_col1 <- avian_flu_df %>%
 group_by(species) %>%
  tally()
colnames(nodes_col1) <- c("label", "value")

nodes_col2 <- avian_flu_df %>%
  group_by(control_measures) %>%
  tally()%>% 
  mutate(control_measures = replace(control_measures, is.na(control_measures), "not applicable"))
colnames(nodes_col2) <- c("label", "value")

nodes_col3 <- avian_flu_df %>%
  group_by(outcome) %>%
  tally()

colnames(nodes_col3) <- c("label", "value")

nodes_col4 <- rbind(nodes_col1, nodes_col2, nodes_col3)

vertices.id <- c(0:(nrow(nodes_col4)-1))

nodes <- cbind.data.frame (vertices.id, nodes_col4)
nodes <- nodes %>%
  mutate(group = case_when (label %in% c(nodes_col1$label) ~ "Host",
                            label %in%  c(nodes_col2$label)~"Control measure",
                            label %in%  c(nodes_col3$label) ~ "Outcome"))

colnames(nodes) <- c("id", "label", "value","group")

#id has to be the same like from and to columns in edges
nodes$id <- nodes$label

# Edges
edges_df <-  avian_flu_df %>%
  group_by(species, control_measures, outcome) %>%
  tally() %>% 
  mutate(control_measures = replace(control_measures, is.na(control_measures), "not applicable"))
  
edges_df1 <- edges_df %>%
  select(species, control_measures, n) 
colnames(edges_df1) <- c("from", "to", "width")

edges_df2 <- edges_df[,-1]
colnames(edges_df2) <- c("from", "to", "width")

edges_df3 <- rbind(edges_df1, edges_df2)
edges_df3 <- as.data.frame(edges_df3)

# Do not weight the edges (only nodes will be)
# Do not run this command if you want to weight the vertices
edges_df3 <- edges_df3 %>%
  select (-width) 
  
# Give font size to the node labels
nodes <- nodes %>% mutate(font.size = 20)

network <- visNetwork(nodes, edges_df3) %>%
  visPhysics(solver = "barnesHut", enabled = FALSE) %>%
  visLegend(enabled = TRUE, useGroups = TRUE) # to add legend

# Save as html file (the html file is saved in your working directory)
visSave(network, "network.html", selfcontained = TRUE, background = "white")
```