# Session 5: Visualising with ggplot2 (and maps)
# Source: https://pierrebeaucoral.github.io/r-for-beginners/
# Open your r-course project in RStudio, then run this file line by line.

Sys.setenv(LANGUAGE = "en")   # R messages in English

library(tidyverse)
library(scales)
library(ggrepel)
library(sf)
library(rnaturalearth)
library(here)

source(here("R", "theme_cours.R"))   # theme_cours(), pal_cours, pal_regions

wdi <- read_csv(here("data", "wdi_panel.csv"), show_col_types = FALSE)

# One row per country for one year: the cross-section used below
wdi_2022 <- wdi |>
  filter(year == 2022)

wdi_2022 |>
  drop_na(life_exp) |>
  slice_min(life_exp, n = 3) |>
  select(country, iso3c, life_exp)

# CAF life expectancy is unreliable in this WDI release (14.7 in 2009, 18.8 in 2022): exclude every year
wdi <- wdi |> mutate(life_exp = if_else(iso3c == "CAF", NA_real_, life_exp))

wdi_2022 <- filter(wdi, year == 2022)

wdi_2022_ok <- wdi_2022 |>
  drop_na(gdp_pc_ppp, life_exp)        # a point needs both coordinates

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp))

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point()

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point() +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale()))

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.8) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_colour_manual(values = pal_regions, guide = "none")

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.8) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_colour_manual(values = pal_regions) +
  guides(colour = guide_legend(nrow = 4, title.position = "top")) +
  labs(
    title = str_wrap("Richer countries live longer, but gains flatten at the top", 60),
    x = "GDP per capita, PPP (constant 2021 international $, log scale)",
    y = "Life expectancy at birth (years)",
    colour = "Region",
    caption = "Source: World Bank, World Development Indicators (WDI), 2022."
  ) +
  theme_cours(base_size = 11)

focus <- c("China", "India", "Nigeria", "United States")

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point(aes(size = population, colour = region), alpha = 0.7) +
  geom_text_repel(
    data = filter(wdi_2022_ok, country %in% focus),
    aes(label = country),
    min.segment.length = 0, box.padding = 1, point.padding = 0.5, seed = 1, size = 4.2,
    bg.color = "white", bg.r = 0.15
  ) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  scale_size_area(max_size = 14, breaks = c(1e7, 1e8, 1e9),
                  labels = c("10 M", "100 M", "1 B")) +
  scale_colour_manual(values = pal_regions) +
  guides(colour = guide_legend(nrow = 3, title.position = "top", order = 1,
                               override.aes = list(size = 4)),
         size = guide_legend(title.position = "top", order = 2)) +
  theme(legend.position = "bottom", legend.box = "vertical",
        legend.box.just = "left", legend.spacing.y = unit(2, "pt")) +
  labs(
    title = str_wrap("Income and life expectancy rise together, but similar incomes can hide a 20-year gap", 60),
    x = "GDP per capita, PPP (constant 2021 international $, log scale)",
    y = "Life expectancy at birth (years)",
    colour = "Region", size = "Population",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. One point per country; countries with missing data are dropped.", 95)
  )

wdi_2022 |>
  drop_na(elec_access, life_exp, income) |>
  ggplot(aes(x = elec_access, y = life_exp, colour = income)) +
  geom_point(alpha = 0.8, size = 2) +
  scale_x_continuous(labels = label_number(suffix = "%")) +
  scale_colour_manual(values = c("Low income" = "#D55E00", "Lower middle income" = "#E69F00",
                                 "Upper middle income" = "#56B4E9", "High income" = "#0072B2")) +
  guides(colour = guide_legend(nrow = 2)) +
  labs(
    title = str_wrap("Countries without universal electricity have markedly lower life expectancy", 60),
    x = "Population with access to electricity (%)",
    y = "Life expectancy at birth (years)",
    colour = str_wrap("Income group (current World Bank classification)", 25),
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022.", 95)
  )

compare <- c("India", "China", "Nigeria", "Brazil", "United States", "Germany", "Japan", "Bangladesh")
highlight <- "India"

life_lines <- wdi |>
  filter(country %in% compare) |>
  drop_na(life_exp)

ggplot(life_lines, aes(x = year, y = life_exp, group = country)) +
  geom_line(data = filter(life_lines, country != highlight),
            colour = pal_cours[["gris"]], linewidth = 0.8) +
  geom_line(data = filter(life_lines, country == highlight),
            colour = pal_cours[["accent"]], linewidth = 1.6) +
  geom_text_repel(data = filter(life_lines, year == max(year)),
                  aes(label = country, colour = country == highlight),
                  hjust = 0, direction = "y", xlim = c(2023.4, NA), size = 4,
                  segment.color = NA, show.legend = FALSE, seed = 1) +
  scale_colour_manual(values = c(`TRUE` = pal_cours[["accent"]], `FALSE` = "grey40")) +
  scale_x_continuous(breaks = seq(2000, 2020, 5), expand = expansion(mult = c(0.02, 0.02))) +
  coord_cartesian(xlim = c(2000, 2027), clip = "off") +
  labs(
    title = str_wrap("India gained nine years of life expectancy since 2000; of these eight countries, only Bangladesh gained more", 60),
    x = "Year",
    y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000-2023. India in red; seven comparison countries in grey.", 95)
  )

life_region <- wdi |>
  drop_na(life_exp) |>
  summarise(life_med = median(life_exp), .by = c(region, year))

ggplot(life_region, aes(x = year, y = life_med, colour = region)) +
  geom_line(linewidth = 1.1) +
  geom_text_repel(
    data = filter(life_region, year == max(year)),
    aes(label = str_wrap(region, 20)),
    hjust = 0, direction = "y", xlim = c(2023.4, NA), size = 3.4, lineheight = 0.85,
    segment.color = NA, show.legend = FALSE
  ) +
  scale_colour_manual(values = pal_regions, guide = "none") +
  scale_x_continuous(breaks = seq(2000, 2020, 5), expand = expansion(mult = c(0.02, 0.05))) +
  coord_cartesian(xlim = c(2000, 2028), clip = "off") +
  labs(
    title = str_wrap("Sub-Saharan Africa gained 12 years since 2000 but still trails every other region", 60),
    x = "Year",
    y = "Median life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000-2023. Median across countries of each region (countries are not weighted by population).", 95)
  )

gdp_region <- wdi_2022_ok |>
  summarise(gdp_med = median(gdp_pc_ppp), n_countries = n(), .by = region)

gdp_region |>
  mutate(region = fct_reorder(region, gdp_med)) |>
  ggplot(aes(x = gdp_med, y = region)) +
  geom_col(fill = pal_cours[["principal"]], width = 0.7) +
  geom_text(aes(label = label_dollar(accuracy = 0.1, scale_cut = cut_short_scale())(gdp_med)), hjust = -0.1, size = 4) +
  scale_x_continuous(labels = label_dollar(scale_cut = cut_short_scale()), expand = expansion(mult = c(0, 0.15))) +
  scale_y_discrete(labels = label_wrap(28)) +
  labs(
    title = str_wrap("A typical North American country is about seven times richer than one in South Asia", 60),
    x = "Median GDP per capita, PPP (constant 2021 international $)",
    y = NULL,
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. Median across countries; North America has only 3 countries.", 95)
  )

niveaux <- c("Low income", "Lower middle income", "Upper middle income", "High income")

wdi_2022 |>
  mutate(income = factor(income, levels = niveaux, ordered = TRUE)) |>
  ggplot(aes(x = income)) +
  geom_bar(fill = pal_cours[["principal"]], width = 0.7) +
  geom_text(aes(label = after_stat(count)), stat = "count", vjust = -0.4, size = 4.5) +
  scale_x_discrete(labels = label_wrap(14)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(
    title = str_wrap("High-income countries are the largest group, but 25 countries are still low-income", 60),
    x = "Income group (current World Bank classification)",
    y = "Number of countries",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. All economies in the WDI panel.", 95)
  )

wdi_2022_ok |>
  slice_max(gdp_pc_ppp, n = 10) |>
  mutate(country = fct_reorder(country, gdp_pc_ppp)) |>
  ggplot(aes(x = gdp_pc_ppp, y = country, fill = region)) +
  geom_col(width = 0.7) +
  scale_x_continuous(labels = label_dollar(scale_cut = cut_short_scale()), expand = expansion(mult = c(0, 0.05))) +
  scale_fill_manual(values = pal_regions) +
  guides(fill = guide_legend(nrow = 4)) +
  labs(
    title = str_wrap("Nine of the ten richest countries are small economies or financial centres; the United States is the exception", 60),
    x = "GDP per capita, PPP (constant 2021 international $)",
    y = NULL, fill = "Region",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022.", 95)
  )

ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp)) +
  geom_histogram(aes(y = after_stat(density)), bins = 25,
                 fill = pal_cours[["principal"]], alpha = 0.6, colour = "white") +
  geom_density(colour = pal_cours[["accent"]], linewidth = 1) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = str_wrap("World incomes span more than a 100-fold range, with no single typical country", 60),
    x = "GDP per capita, PPP (constant 2021 international $, log scale)",
    y = "Density",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. One observation per country; red line = kernel density estimate.", 95)
  )

wdi_2022_ok |>
  mutate(income = factor(income, levels = niveaux, ordered = TRUE)) |>
  ggplot(aes(x = income, y = gdp_pc_ppp)) +
  geom_boxplot(outliers = FALSE, fill = "grey92", colour = "grey30", width = 0.55) +
  geom_jitter(width = 0.15, alpha = 0.5, colour = pal_cours[["principal"]], size = 1.6) +
  scale_x_discrete(labels = label_wrap(14)) +
  scale_y_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = str_wrap("Median income rises steeply from one income group to the next, but neighbouring groups overlap", 60),
    x = "Income group (current World Bank classification)",
    y = str_wrap("GDP per capita, PPP (constant 2021 international $, log scale)", 40),
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. Box = quartiles and median; dots = individual countries.", 95)
  )

wdi |>
  filter(year %in% c(2000, 2022)) |>
  drop_na(life_exp) |>
  ggplot(aes(x = life_exp)) +
  geom_histogram(binwidth = 2.5, fill = pal_cours[["principal"]], colour = "white") +
  facet_wrap(~ year) +
  labs(
    title = str_wrap("The low-life-expectancy tail shrank between 2000 and 2022", 60),
    x = "Life expectancy at birth (years)",
    y = "Number of countries",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000 and 2022. One observation per country and year.", 95)
  )

world <- ne_countries(scale = "medium", returnclass = "sf") |>
  filter(name != "Antarctica")          # no data there, and it eats a fifth of the map

dim(world)
world |>
  st_drop_geometry() |>
  select(name, iso_a3, iso_a3_eh, adm0_a3) |>
  slice(1:5)

# WDI countries with no match in the map, for each candidate code
unmatched <- function(code) {
  wdi_2022 |>
    anti_join(st_drop_geometry(world), by = c("iso3c" = code)) |>
    pull(country)
}

# How many WDI economies stay unmatched with each code?
length(unmatched("iso_a3"))
length(unmatched("iso_a3_eh"))
length(unmatched("adm0_a3"))

unmatched("iso_a3_eh")

world |>
  st_drop_geometry() |>
  anti_join(wdi_2022, by = c("iso_a3_eh" = "iso3c")) |>
  select(name, iso_a3_eh) |>
  head(8)

world_gdp <- world |>
  select(iso_a3_eh, geometry) |>
  left_join(select(wdi_2022, iso3c, country, gdp_pc_ppp, elec_access),
            by = c("iso_a3_eh" = "iso3c"))

sum(is.na(world_gdp$gdp_pc_ppp))   # shapes drawn in grey (no data)

cuts   <- c(0, 2500, 5000, 10000, 20000, 40000, Inf)
labs_g <- c("< $2,500", "$2,500-5,000", "$5,000-10,000", "$10,000-20,000",
            "$20,000-40,000", "> $40,000")

world_gdp |>
  mutate(gdp_class = cut(gdp_pc_ppp, breaks = cuts, labels = labs_g, right = FALSE)) |>
  ggplot() +
  geom_sf(aes(fill = fct_na_value_to_level(gdp_class, "No data")),   # NA becomes a level
          colour = "white", linewidth = 0.1) +
  scale_fill_manual(values = c(hcl.colors(6, "YlGnBu", rev = TRUE), "grey85"),  # one colour per level, in order
                    drop = FALSE) +
  coord_sf(crs = "+proj=robin") +
  guides(fill = guide_legend(nrow = 2, title.position = "top")) +
  labs(
    x = NULL, y = NULL,
    title = str_wrap("Income per head is highest in the North and lowest in Sub-Saharan Africa", 60),
    fill = "GDP per capita, PPP (constant 2021 international $)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. Light grey: no data.", 95)
  ) +
  theme(axis.text = element_blank(), panel.grid.major = element_line(colour = "grey94"))

bottom3 <- world_gdp |>
  drop_na(elec_access) |>
  slice_min(elec_access, n = 3, with_ties = FALSE) |>
  mutate(label = paste0(str_remove(country, ",.*"), " ", round(elec_access), "%")) |>
  st_point_on_surface()                       # one point inside each polygon

bottom3$label                                  # check which countries are annotated

world_gdp |>
  mutate(elec_class = cut(elec_access, breaks = c(0, 50, 80, 99, 100.01),
                          labels = c("< 50%", "50-80%", "80-99%", "99-100%"),
                          right = FALSE)) |>
  ggplot() +
  geom_sf(aes(fill = fct_na_value_to_level(elec_class, "No data")),
          colour = "white", linewidth = 0.1) +
  geom_label_repel(data = bottom3, aes(label = label, geometry = geometry),
                   stat = "sf_coordinates", size = 3.4, fill = alpha("white", 0.85),
                   min.segment.length = Inf, box.padding = 0.2, seed = 1) +  # no leader lines
  scale_fill_manual(values = c("#08306b", "#2171b5", "#9ecae1", "#eff3ff", "grey85"),  # in level order
                    drop = FALSE) +
  coord_sf(crs = "+proj=robin") +
  guides(fill = guide_legend(nrow = 2, title.position = "top")) +
  labs(
    x = NULL, y = NULL,
    title = str_wrap("20 of the 22 countries where less than half the population has electricity are in Sub-Saharan Africa", 60),
    fill = "Population with access to electricity",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. Labels (placed near the countries): the three with the lowest access. Light grey: no data.", 95)
  ) +
  theme(axis.text = element_blank(), panel.grid.major = element_line(colour = "grey94"))

# dir.create(here("output"), showWarnings = FALSE)   # output/ is created if missing
# 
# p_gapminder <- ggplot(wdi_2022_ok, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
#   geom_point() +
#   scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale()))
# 
# ggsave(here("output", "gapminder.png"), plot = p_gapminder,
#        width = 8, height = 5, units = "in", dpi = 300)    # for Word / slides
# ggsave(here("output", "gapminder.pdf"), plot = p_gapminder,
#        width = 8, height = 5)                             # vector, for LaTeX

p_life_region <- wdi_2022 |>
  drop_na(life_exp) |>
  mutate(region = fct_reorder(region, life_exp, .fun = median)) |>
  ggplot(aes(x = life_exp, y = region)) +
  geom_boxplot(fill = "grey92", colour = "grey30", width = 0.55) +
  scale_y_discrete(labels = label_wrap(28)) +
  labs(
    title = str_wrap("Sub-Saharan Africa has the lowest median life expectancy, about 9 years below the next region", 60),
    x = "Life expectancy at birth (years)",
    y = NULL,
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022. One observation per country.", 95)
  )

p_life_region

# dir.create(here("output"), showWarnings = FALSE)
# ggsave(here("output", "life_exp_region.png"), plot = p_life_region,
#        width = 8, height = 5, dpi = 300)