---
title: "lab7_slides_supp"
author: "Brian Leung"
date: "`r Sys.Date()`"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

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

## Prerequisite

```{r}
# download packages 
packages <- c("tidyverse", "tidygraph", "ggraph", "reshape2", "cluster", "circlize")
not_installed <- setdiff(packages, rownames(installed.packages()))
if (length(not_installed)) install.packages(not_installed)

# install developmental version of ggsankey
devtools::install_github("davidsjoberg/ggsankey")

# load packages
library(tidyverse)
library(tidygraph)
library(ggraph)
library(reshape2)
library(cluster)
library(circlize)
library(ggsankey)
```

## Medici
```{r}
# load data
medici <- read.table("data/medici.txt")

# coerce into a tibble graph
medici_graph <- as_tbl_graph(medici, directed = FALSE)

# preliminary visualization
ggraph(medici_graph) +
  geom_node_point() +
  geom_edge_link() +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_graph()

# add network measures
medici_graph <-
  medici_graph %>%
  mutate(
    degree = centrality_degree(), # Calculate degree centrality
    community = group_leading_eigen() # Implement community-detection algorithm
  )

# visualization
ggraph(medici_graph) +
  geom_edge_link(alpha = 0.5) +
  geom_node_point(aes(size = degree, color = factor(community)), show.legend = FALSE) +
  scale_color_brewer(palette = "Set1") +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_graph()

# save output
width = 10
ggsave("figures/medici.png", width = width, height = width/2)
```

## Global migration flow data
```{r}
# load data
migrat2010 <- read_csv("data/migrat2010.csv")

# preliminary visualization
migrat2010_graph <- as_tbl_graph(migrat2010)

# unidirectional graph
ggraph(migrat2010_graph) +
  geom_edge_link(alpha = 0.5) +
  geom_node_point() +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_graph()

# bidirectional graph
ggraph(migrat2010_graph) +
  geom_edge_parallel(
    start_cap = circle(1.25, 'mm'),
    end_cap = circle(1.25, 'mm'),
    arrow = arrow(length = unit(2, 'mm')),
    sep = unit(1.25, 'mm'),
    alpha = 0.5
  ) +
  geom_node_point() +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_graph()

# heatmap
migrat2010 %>%
  ggplot(aes(y = origRegion, x = destRegion, fill = flow)) +
  geom_tile(color = "white", size = 0.2) +
  coord_equal() +
  theme(panel.background = element_blank())

# Three wrangling tasks:
# 1. Make NA values explicit
# 2. Turn flow into a categorical variable
# 3. Cluster analysis and sorting

quantile(migrat2010$flow, na.rm = TRUE)
breaks <- c(1000, 5000, 10000, 50000, 100000, Inf)
labels <- c("1000-5000", "5000-10000", "10000-50000", "50000-100000", ">100000")

migrat2010 <- 
  migrat2010 %>%
  expand(origRegion, destRegion) %>%
  left_join(migrat2010, by = c("origRegion", "destRegion")) %>%
  mutate(flowCat = cut(flow, breaks = breaks, labels = labels))

# cluster analysis
migrat2010_matrix <-
  migrat2010 %>%
  reshape2::acast(origRegion ~ destRegion, value.var = "flow", fill = 0)

migrat2010_hclust <- 
  migrat2010_matrix %>%
  dist() %>%
  hclust(method = "ward.D")

plot(migrat2010_hclust)

# check the order
countryOrder <- migrat2010_hclust$order
names(countryOrder) <- migrat2010_hclust$labels
print(countryOrder)

# Sort the countries using the order produced by cluster analysis
countryLevels <- row.names(migrat2010_matrix)[countryOrder]
 
# Re-level `origRegion` and `destRegion` according to the level
migrat2010 <- 
  migrat2010 %>%
  mutate(
    origRegion = factor(origRegion, levels = rev(countryLevels)),
    destRegion = factor(destRegion, levels = countryLevels)
  ) %>%
  mutate(
    Asia = str_detect(origRegion, "Asia")
  )

# visualize in heatmap 
migrat2010 %>%
  ggplot(aes(y = origRegion, x = destRegion, fill = flowCat)) +
  geom_tile(color = "white", size = 0.2) +
  # Scale fill values with "Blues" palette and "grey90" for NAs
  scale_fill_brewer(palette = "Blues",
                    na.value = "grey90", 
                    breaks = c(rev(labels), NA)) +
  # Put x-axis labels on top
  scale_x_discrete(position = "top") +
  #coord_equal() +
  theme(
    panel.background = element_blank(),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    # Rotate and align x-axis labels
    axis.text.x.top = element_text(angle = 90, hjust = 0),
    legend.key.height = grid::unit(0.8, "cm"),
    legend.key.width = grid::unit(0.2, "cm")
  ) +
  guides(fill = guide_legend(title = "Migration flow")) +
  labs(y = "Origin", x = "Destination") + 
  facet_grid(Asia ~ ., scales = "free", space = "free")

width = 8
ggsave("figures/migrat_heatmap.png", width = width, height = width)
```

## Chord diagram
```{r}
# Create vectors of countries to be aggregated
Europe <- c("Southern Europe", "Western Europe", "Northern Europe")
CAEE <- c("Eastern Europe", "Central Asia")
Africa <- c("Eastern Africa", "Middle Africa", "Northern Africa",
"Southern Africa", "Western Africa")
LACarib <- c("Central America", "South America", "Caribbean")
SAsia <- c("South-Eastern Asia", "Southern Asia")

# Use mutate() and across() to recode
# `origRegion` and `destRegion` simultaneously
migrat2010 <-
  migrat2010 %>%
  mutate(across(c(origRegion, destRegion),
                ~ case_when(.x %in% Europe ~ "Europe",
                            .x %in% CAEE ~ "Central Asia \n& Eastern Europe",
                            .x %in% Africa ~ "Africa",
                            .x %in% LACarib ~ "Latin America \n& Caribbean",
                            .x %in% SAsia ~ "Southern Asia",
                            TRUE ~ as.character(.x)))) %>%
  mutate(across(c(origRegion, destRegion),
                ~ factor(.x, levels = c("Southern Asia",
                                        "Eastern Asia",
                                        "Western Asia",
                                        "Central Asia \n& Eastern Europe",
                                        "Europe",
                                        "Africa",
                                        "Latin America \n& Caribbean",
                                        "Northern America",
                                        "Oceania"))))

# Collapse (sum) flow values according by newly aggregated regions
migrat2010 <- 
  migrat2010 %>%
  group_by(origRegion, destRegion) %>%
  summarize(flowTotal = sum(flow, na.rm = TRUE)) %>%
  ungroup()

# Basic chord diagram
chordDiagram(migrat2010)

# Setting parameters
circos.clear()
circos.par(
  start.degree = 90, # Start at 12 o'clock
  gap.degree = 4, # Increase gaps between sectors
  track.margin = c(-0.1, 0.1), # Narrow the track margin
  points.overflow.warning = FALSE # Subdue warning messages
)
par(mar = rep(0, 4)) # no margins in the plot

# Get nice colors
colors <- RColorBrewer::brewer.pal(9, "Paired")

# More advanced settings in `chordDiagram()`
chordDiagram(
  migrat2010,
  grid.col = colors, # Set colors
  directional = 1, # Indicate chords are directional
  direction.type = c("arrows", "diffHeight"), # Directionality is illustrated by arrows and height differences
  diffHeight = -0.04, # Set height difference
  link.arr.type = "big.arrow", # Use big arrows
  link.sort = TRUE, # Sort the chords and plot the smallest chords first
  link.largest.ontop = TRUE,
)

# Save the output
dev.copy2pdf(file = "figures/migrat_chord.pdf", height = 8, width = 8)
```


## Sankey diagram
```{r}
migrat2010_sankey <- 
  migrat2010 %>%
  mutate(origRegion = str_remove(origRegion, "\n"),
         destRegion = str_remove(destRegion, "\n")) %>%
  make_long(origRegion, destRegion, value = flowTotal) 
  

ggplot(migrat2010_sankey, aes(x = x, next_x = next_x, 
                              node = node, next_node = next_node, 
                              value = value,
                              fill = node)) +
  geom_sankey(flow.alpha = 0.6, width = 0.01) +
  geom_sankey_text(data = filter(migrat2010_sankey, x == "origRegion"),
                   aes(label = paste0(node, "   "), color = node), hjust = 1) +
  geom_sankey_text(data = filter(migrat2010_sankey, x == "destRegion"),
                   aes(label = paste0("   ", node), color = node), hjust = 0) +
  scale_fill_manual(values = colors) +
  scale_color_manual(values = colors) +
  theme_sankey() +
  labs(x = NULL) +
  theme(legend.position = "none") +
  scale_x_discrete(labels = c("Origin", "Destination"))

ggsave("figures/migrat_sankey.png", width = 9, height = 9)
``` 













