# Session 6: Advanced visualisation and Quarto reports
# 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")   # messages in English

library(tidyverse)
library(here)
library(scales)
library(patchwork)
library(plotly)
library(crosstalk)

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

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

# 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))

# Cross-section used in most examples: latest year with good coverage
wdi_2022 <- wdi |>
  filter(year == 2022) |>
  drop_na(gdp_pc_ppp, life_exp, region)

ggplot(wdi_2022, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point(aes(colour = region), alpha = 0.7, show.legend = FALSE) +
  scale_colour_manual(values = pal_regions) +
  scale_x_log10(breaks = c(1e3, 1e4, 1e5),
                labels = label_dollar(scale_cut = cut_short_scale())) +
  facet_wrap(~ region, ncol = 3, labeller = label_wrap_gen(20)) +
  labs(
    title = "Richer countries live longer, in every region",
    x = "GDP per capita, PPP (constant 2021 int. $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2022. One point = one country."
  )

trend_region <- wdi |>
  drop_na(life_exp, region) |>
  summarise(life_exp = mean(life_exp), .by = c(region, year))

ggplot(trend_region, aes(x = year, y = life_exp)) +
  geom_line(aes(colour = region), linewidth = 1, show.legend = FALSE) +
  scale_colour_manual(values = pal_regions) +
  scale_x_continuous(breaks = c(2000, 2010, 2020)) +
  facet_wrap(~ region, ncol = 3, labeller = label_wrap_gen(20)) +
  labs(
    title = "Life expectancy rose in every region, despite a COVID-19 dip",
    x = NULL, y = "Life expectancy at birth (years)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2000-2023. Unweighted mean of countries in each region."
  )

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

wdi |>
  filter(year %in% c(2000, 2022), income %in% income_levels) |>
  drop_na(gdp_pc_ppp, life_exp) |>
  mutate(income = factor(income, levels = income_levels)) |>
  ggplot(aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point(colour = pal_cours[["principal"]], alpha = 0.6) +
  scale_x_log10(breaks = c(1e3, 1e4, 1e5),
                labels = label_dollar(scale_cut = cut_short_scale())) +
  facet_grid(income ~ year, labeller = label_wrap_gen(10)) +
  labs(
    title = str_wrap("Within income groups, the income-health link looks weaker", 60),
    x = "GDP per capita, PPP (constant 2021 int. $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000 and 2022. One point = one country. Rows: income group (current World Bank classification).", 95)
  )

trend_region |>
  filter(region %in% c("Sub-Saharan Africa", "North America",
                       "South Asia", "Europe & Central Asia")) |>
  ggplot(aes(x = year, y = life_exp)) +
  geom_line(linewidth = 1, colour = pal_cours[["principal"]]) +
  facet_wrap(~ region, scales = "free_y", nrow = 1) +
  labs(
    title = str_wrap("Free y-axes: North America's 3-year gain looks as large as Africa's 11-year gain", 60),
    subtitle = str_wrap("Careful: the panels no longer share a scale", 75),
    x = NULL, y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000-2023. Unweighted mean of countries.", 95)
  )

trend_region |>
  filter(region %in% c("Sub-Saharan Africa", "North America",
                       "South Asia", "Europe & Central Asia")) |>
  ggplot(aes(x = year, y = life_exp)) +
  geom_line(linewidth = 1, colour = pal_cours[["principal"]]) +
  facet_wrap(~ region, nrow = 1) +
  labs(
    title = "Fixed y-axes: the regions differ in level, and Africa and South Asia gain most",
    x = NULL, y = "Life expectancy at birth (years)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2000-2023. Unweighted mean of countries."
  )

elec <- wdi |>
  filter(income %in% income_levels) |>
  drop_na(elec_access) |>
  mutate(income = factor(income, levels = income_levels))

elec_mean <- elec |>
  summarise(elec_access = mean(elec_access), .by = c(income, year))

ggplot(elec, aes(x = year, y = elec_access)) +
  geom_line(aes(group = iso3c), colour = pal_cours[["gris"]],
            linewidth = 0.3, alpha = 0.6) +
  geom_line(data = elec_mean, colour = pal_cours[["principal"]], linewidth = 1.2) +
  scale_x_continuous(breaks = c(2000, 2010, 2020)) +
  facet_wrap(~ income, labeller = label_wrap_gen(20)) +
  labs(
    title = "Electricity access is near-universal in rich countries, not in poor ones",
    x = NULL, y = "Population with access to electricity (%)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2000-2023. Grey: countries; blue: unweighted group mean."
  )

# Before: defaults only
ggplot(wdi_2022, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point() +
  theme_gray()

# After: message title, units, emphasis, source
wdi_2022 |>
  mutate(focus = region == "Sub-Saharan Africa") |>
  ggplot(aes(x = gdp_pc_ppp, y = life_exp, colour = focus)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = c("FALSE" = pal_cours[["gris"]],
                                 "TRUE"  = pal_cours[["accent"]]),
                      guide = "none") +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = str_wrap("Sub-Saharan Africa (red) sits at the low end of income and health", 45),
    x = "GDP per capita, PPP (2021 int. $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2022.", 55)
  )

p_base <- ggplot(wdi_2022, aes(x = gdp_pc_ppp, y = life_exp)) +
  geom_point(colour = pal_cours[["principal"]], alpha = 0.6) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = "Life expectancy rises with income",
    x = "GDP per capita, PPP (constant 2021 int. $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2022."
  )

p_base +
  theme(
    plot.title = element_text(colour = pal_cours[["accent"]], size = 18),
    panel.grid.minor = element_blank(),
    axis.title.x = element_text(margin = margin(t = 10))
  )

theme_cours_dark <- function(base_size = 14) {
  theme_cours(base_size = base_size) +
    theme(
      plot.background  = element_rect(fill = "#1b2733", colour = NA),
      panel.background = element_rect(fill = "#1b2733", colour = NA),
      text             = element_text(colour = "grey90"),
      plot.title       = element_text(colour = "white"),
      plot.subtitle    = element_text(colour = "grey80"),
      plot.caption     = element_text(colour = "grey70"),
      axis.text        = element_text(colour = "grey80"),
      axis.title       = element_text(colour = "grey85"),
      panel.grid.major = element_line(colour = "grey30")
    )
}

p_base + geom_point(colour = "#56B4E9", alpha = 0.7) + theme_cours_dark()

# What is the current default theme? (a theme is a list of settings: print one of them)
theme_get()$legend.position

theme_cours_print <- function(base_size = 12) {
  theme_cours(base_size = base_size) +
    theme(
      legend.position    = "top",
      panel.grid.major.x = element_blank(),
      panel.grid.minor   = element_blank()
    )
}

ggplot(wdi_2022, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = pal_regions, name = NULL) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  guides(colour = guide_legend(nrow = 3)) +
  labs(
    title = "Life expectancy rises with income in every region",
    x = "GDP per capita, PPP (constant 2021 int. $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = "Source: World Bank, World Development Indicators (WDI), 2022."
  ) +
  theme_cours_print() +
  theme(legend.text = element_text(size = 9))

life_world <- wdi |>
  drop_na(life_exp, population) |>
  summarise(life_exp = weighted.mean(life_exp, w = population), .by = year)

le_2019_2021 <- life_world |>
  filter(year %in% c(2019, 2021, 2023)) |>
  arrange(year) |>
  pull(life_exp)

ggplot(life_world, aes(x = year, y = life_exp)) +
  annotate("rect", xmin = 2019.5, xmax = 2021.5, ymin = -Inf, ymax = Inf,
           fill = pal_cours[["accent"]], alpha = 0.08) +
  geom_line(linewidth = 1.1, colour = pal_cours[["principal"]]) +
  geom_vline(xintercept = 2020, linetype = "dashed", colour = pal_cours[["accent"]]) +
  annotate("text", x = 2019.3, y = min(life_world$life_exp) + 0.5,
           label = "COVID-19\n(2020-21)", hjust = 1,
           colour = pal_cours[["accent"]], fontface = "bold") +
  labs(
    title = str_wrap("Global life expectancy fell after 2019 and recovered by 2023", 60),
    x = NULL, y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2000-2023. Population-weighted mean of countries.", 95)
  )

p_scatter <- ggplot(wdi_2022, aes(x = gdp_pc_ppp, y = life_exp, colour = region)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = pal_regions, name = NULL) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(x = "GDP per capita, PPP (2021 int. $, log)", y = "Life exp. (years)")

p_box <- wdi_2022 |>
  ggplot(aes(x = fct_reorder(region, life_exp), y = life_exp, fill = region)) +
  geom_boxplot(show.legend = FALSE) +
  scale_fill_manual(values = pal_regions) +
  coord_flip() +
  labs(x = NULL, y = "Life expectancy (years)")

p_hist <- ggplot(wdi_2022, aes(x = life_exp)) +
  geom_histogram(binwidth = 2, fill = pal_cours[["principal"]], colour = "white") +
  labs(x = "Life expectancy (years)", y = "Count")

((p_scatter | p_box) / p_hist) +
  plot_layout(guides = "collect", heights = c(2, 1)) +
  plot_annotation(
    title = "Life expectancy in 2022: rising with income, with large regional gaps",
    caption = "Source: World Bank, World Development Indicators (WDI), 2022. One observation = one country.",
    tag_levels = "A"
  ) &
  theme(legend.position = "bottom", legend.text = element_text(size = 9))

# Only panel B keeps a legend: both panels use the same colours
pA <- ggplot(trend_region, aes(x = year, y = life_exp, colour = region)) +
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_colour_manual(values = pal_regions) +
  labs(x = NULL, y = "Life expectancy (years)")

pB <- wdi |>
  filter(year == 2021) |>
  drop_na(elec_access, gdp_pc_ppp, region) |>
  ggplot(aes(x = gdp_pc_ppp, y = elec_access, colour = region)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = pal_regions, name = NULL,
                      labels = label_wrap(24)) +
  scale_x_log10(breaks = c(1e3, 1e4, 1e5),
                labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(x = "GDP per capita, PPP (2021 int. $, log)",
       y = "Access to electricity (% of population)")

(pA | pB) +
  plot_layout(guides = "collect") +
  plot_annotation(
    title = "Health improved everywhere; electricity access depends on income",
    caption = str_wrap("Source: World Bank, World Development Indicators (WDI). A: unweighted regional means, 2000-2023; B: countries, 2021.", 100),
    tag_levels = "A"
  ) &
  theme(legend.position = "right", legend.text = element_text(size = 9))

# dir.create(here("output"), showWarnings = FALSE)
# 
# fig <- p_scatter + labs(title = "Life expectancy rises with income",
#                         caption = "Source: World Bank, World Development Indicators (WDI), 2022.")
# 
# # Raster for Word / the web: PNG at 300 dpi
# ggsave(here("output", "life_exp_income.png"), plot = fig,
#        width = 16, height = 10, units = "cm", dpi = 300)
# 
# # Vector for LaTeX / print: PDF (scales without blur)
# ggsave(here("output", "life_exp_income.pdf"), plot = fig,
#        width = 16, height = 10, units = "cm")

wdi_hover <- wdi_2022 |>
  mutate(
    tip = paste0(
      "<b>", country, "</b><br>",
      "GDP per capita: ", label_dollar()(round(gdp_pc_ppp)), "<br>",
      "Life expectancy: ", round(life_exp, 1), " years"
    )
  )

p_int <- ggplot(wdi_hover, aes(x = gdp_pc_ppp, y = life_exp,
                               colour = region, label = tip)) +
  geom_point(alpha = 0.7) +
  scale_colour_manual(values = pal_regions, name = NULL) +
  scale_x_log10(labels = label_dollar(scale_cut = cut_short_scale())) +
  labs(
    title = "Hover over a point to identify the country",
    x = "GDP per capita, PPP (2021 int. $, log scale)",
    y = "Life expectancy at birth (years)"
  )

ggplotly(p_int, tooltip = "label") |>
  layout(height = 520, legend = list(orientation = "h", x = 0, y = -0.25),
         margin = list(b = 140, t = 60), font = list(family = "IBM Plex Sans"),
         paper_bgcolor = "#fcfaf6", plot_bgcolor = "#fcfaf6")

shared_wdi <- SharedData$new(wdi_hover, key = ~iso3c)

p_left <- plot_ly(shared_wdi, x = ~gdp_pc_ppp, y = ~life_exp,
                  color = ~region, colors = unname(pal_regions),
                  text = ~tip, hoverinfo = "text",
                  type = "scatter", mode = "markers") |>
  layout(xaxis = list(title = "GDP per capita, PPP (2021 int. $)", type = "log"),
         yaxis = list(title = "Life expectancy (years)"),
         showlegend = FALSE)

p_right <- plot_ly(shared_wdi, x = ~elec_access, y = ~life_exp,
                   color = ~region, colors = unname(pal_regions),
                   text = ~tip, hoverinfo = "text",
                   type = "scatter", mode = "markers") |>
  layout(xaxis = list(title = "Access to electricity (% of population)"),
         yaxis = list(title = "Life expectancy (years)"),
         showlegend = FALSE)

bscols(
  widths = c(12, 6, 6),
  filter_select("region", "Region", shared_wdi, ~region),
  p_left, p_right
)

p_facet_int <- ggplot(wdi_hover, aes(x = gdp_pc_ppp, y = life_exp,
                                     colour = region, label = tip)) +
  geom_point(alpha = 0.7, show.legend = FALSE) +
  scale_colour_manual(values = pal_regions) +
  scale_x_log10(breaks = c(1e3, 1e4, 1e5),
                labels = label_dollar(scale_cut = cut_short_scale())) +
  facet_wrap(~ region, ncol = 3, labeller = label_wrap_gen(20)) +
  labs(x = "GDP per capita, PPP (2021 int. $, log scale)",
       y = "Life expectancy (years)")

ggplotly(p_facet_int, tooltip = "label")

# library(shiny)
# library(tidyverse)
# library(here)
# 
# wdi <- read_csv(here("data", "wdi_panel.csv"), show_col_types = FALSE)
# 
# num_vars <- c(
#   "GDP per capita (PPP)"         = "gdp_pc_ppp",
#   "Life expectancy (years)"      = "life_exp",
#   "Gini index"                   = "gini",
#   "Access to electricity (%)"    = "elec_access"
# )
# 
# ui <- fluidPage(
#   titlePanel("WDI explorer"),
#   sidebarLayout(
#     sidebarPanel(
#       selectInput("xvar", "X-axis", choices = num_vars, selected = "gdp_pc_ppp"),
#       selectInput("yvar", "Y-axis", choices = num_vars, selected = "life_exp"),
#       selectInput("region", "Region",
#                   choices = c("All", sort(unique(wdi$region)))),
#       sliderInput("year", "Year", min = 2000, max = 2023, value = 2022,
#                   step = 1, sep = "")
#     ),
#     mainPanel(plotOutput("scatter"))
#   )
# )
# 
# server <- function(input, output) {
# 
#   # Reactive: recomputed only when region or year changes
#   wdi_sel <- reactive({
#     wdi |>
#       filter(year == input$year,
#              input$region == "All" | region == input$region)
#   })
# 
#   output$scatter <- renderPlot({
#     # .data[[ ]] reads a column whose NAME is stored in a string
#     # (this replaces the deprecated aes_string())
#     ggplot(wdi_sel(), aes(x = .data[[input$xvar]], y = .data[[input$yvar]])) +
#       geom_point(colour = "#1f4e79", alpha = 0.7, na.rm = TRUE) +
#       labs(
#         title = paste(names(num_vars)[num_vars == input$yvar], "vs",
#                       names(num_vars)[num_vars == input$xvar]),
#         x = names(num_vars)[num_vars == input$xvar],
#         y = names(num_vars)[num_vars == input$yvar],
#         caption = "Source: World Bank, WDI."
#       ) +
#       theme_minimal(base_size = 14)
#   })
# }
# 
# shinyApp(ui = ui, server = server)

elec_ssa <- wdi |>
  filter(region == "Sub-Saharan Africa", year %in% c(2000, 2022)) |>
  select(country, year, elec_access) |>
  pivot_wider(names_from = year, values_from = elec_access,
              names_prefix = "y") |>
  drop_na(y2000, y2022)

lowest12 <- elec_ssa |>
  slice_min(y2022, n = 12, with_ties = FALSE) |>
  mutate(country = fct_reorder(country, y2022),
         gain = y2022 - y2000)

mean_gain <- round(mean(lowest12$gain))

fig_note <- ggplot(lowest12, aes(y = country)) +
  geom_segment(aes(x = y2000, xend = y2022, yend = country),
               colour = pal_cours[["gris"]], linewidth = 1.2) +
  geom_point(aes(x = y2000, colour = "2000"), size = 3) +
  geom_point(aes(x = y2022, colour = "2022"), size = 3) +
  scale_colour_manual(
    values = c("2000" = pal_cours[["gris"]], "2022" = pal_cours[["accent"]]),
    name = NULL
  ) +
  scale_x_continuous(limits = c(0, 100), labels = label_percent(scale = 1)) +
  labs(
    title = str_wrap(paste0("Where access is lowest, coverage still grew by ",
                            mean_gain, " points on average since 2000"), 60),
    subtitle = str_wrap("Population with access to electricity, the 12 Sub-Saharan countries with the lowest access in 2022", 75),
    x = "Population with access to electricity (%)", y = NULL,
    caption = str_wrap(paste("Source: World Bank, World Development Indicators (WDI), 2000 and 2022.",
                             "Countries with data in both years; no population weighting."), 90)
  )

fig_note

# dir.create(here("output"), showWarnings = FALSE)
# ggsave(here("output", "note_ministre_electricity.png"), plot = fig_note,
#        width = 16, height = 11, units = "cm", dpi = 300)