# Session 3: Data manipulation (dplyr, tidyr, joins)
# 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)   # dplyr, tidyr, forcats, ggplot2, readr...
library(here)
source(here("R", "theme_cours.R"))

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

glimpse(wdi)

# Without the pipe: read from the inside out
head(arrange(filter(wdi, year == 2020), desc(life_exp)), 3)

# With the pipe: read from top to bottom
wdi |>
  filter(year == 2020) |>
  arrange(desc(life_exp)) |>
  head(3)

wdi |>
  filter(year == 2020, region == "South Asia") |>
  select(country, life_exp, gdp_pc_ppp)

wdi |>
  filter(year == 2020, iso3c %in% c("FRA", "DEU", "ESP")) |>
  select(country, life_exp)

wdi |> select(country, year, life_exp) |> head(3)

# Helpers: select columns by pattern or by type
wdi |> select(country, starts_with("gdp")) |> head(3)
wdi |> select(contains("_")) |> head(3)
wdi |> select(where(is.numeric)) |> names()

wdi |>
  filter(year == 2020) |>
  arrange(desc(gdp_pc_ppp)) |>   # desc() = largest first
  select(country, gdp_pc_ppp) |>
  head(5)

wdi |>
  filter(year == 2020) |>
  mutate(gdp_bn = gdp / 1e9,                      # GDP in billions of US$
         pop_m  = population / 1e6) |>            # population in millions
  select(country, gdp_bn, pop_m) |>
  arrange(desc(gdp_bn)) |>
  head(3)

wdi |> rename(life_expectancy = life_exp) |> names()

wdi_2020 <- wdi |>
  filter(year == 2020) |>
  mutate(
    high_life = if_else(life_exp >= 75, "75+ years", "under 75"),
    elec_level = case_when(
      elec_access >= 99 ~ "universal (99%+)",
      elec_access >= 80 ~ "high (80-99%)",
      elec_access >= 50 ~ "medium (50-80%)",
      elec_access <  50 ~ "low (<50%)",
      .default = "no data"
    )
  )

wdi_2020 |> count(elec_level)

wdi_2020 |>
  count(elec_level) |>
  ggplot(aes(x = n, y = fct_reorder(elec_level, n))) +
  geom_col(fill = pal_cours[["principal"]]) +
  labs(title = "Most countries have near-universal electricity access",
       x = "Number of countries, 2020", y = NULL,
       caption = "Source: World Bank, World Development Indicators (WDI), 2020.")

try({
# wdi |>
#   filter(year = 2020, region == "Sub-Saharan africa") |>
#   select(country, life_exp, gdp_pc) |>
#   arrange(desc(life_exp))
})

try({
# wdi |>
#   filter(year == 2020, region == "Sub-Saharan africa") |>
#   select(country, life_exp, gdp_pc) |>
#   arrange(desc(life_exp))
})

wdi |>
  filter(year == 2020, region == "Sub-Saharan africa") |>
  select(country, life_exp, gdp_pc_ppp) |>
  arrange(desc(life_exp))

wdi |> nrow()                                    # start: all country-years
wdi |> filter(year == 2020) |> nrow()            # one year: one row per country
wdi |> filter(year == 2020, region == "Sub-Saharan africa") |> nrow()   # 0: the problem is here

wdi |> distinct(region)

wdi |>
  filter(year == 2020, region == "Sub-Saharan Africa") |>
  select(country, life_exp, gdp_pc_ppp) |>
  arrange(desc(life_exp)) |>
  glimpse()

# wdi |>
#   filter(year == 2019, income == "Low Income") |>
#   mutate(gdp_pc_k = gdp_pc_ppp / 1000) |>
#   select(country, gdp_pc_k, life_expectancy) |>
#   arrange(desc(gdp_pc_k))

# Bug 1 (error in select): the column is life_exp, not life_expectancy
names(wdi)

# Bug 2 (silent): no row survives the filter
wdi |> filter(year == 2019, income == "Low Income") |> nrow()
wdi |> distinct(income)   # the label is "Low income", lower-case i

low_2019 <- wdi |>
  filter(year == 2019, income == "Low income") |>
  mutate(gdp_pc_k = gdp_pc_ppp / 1000) |>
  select(country, gdp_pc_k, life_exp) |>
  arrange(desc(gdp_pc_k))

nrow(low_2019)   # N
low_2019

wdi |>
  filter(year == 2020) |>
  group_by(region) |>
  summarise(
    n_countries = n(),
    median_life = median(life_exp, na.rm = TRUE),
    median_gdp_pc = median(gdp_pc_ppp)
  )

life_by_region <- wdi |>
  filter(year == 2020) |>
  group_by(region) |>
  summarise(mean_life = mean(life_exp, na.rm = TRUE), n_countries = sum(!is.na(life_exp)))

life_by_region

life_by_region |>
  ggplot(aes(x = mean_life, y = fct_reorder(region, mean_life))) +
  geom_col(fill = pal_cours[["principal"]]) +
  geom_text(aes(label = paste0("n = ", n_countries)), hjust = 1.1, colour = "white") +
  labs(title = str_wrap("Sub-Saharan Africa trails North America by 17 years", 60),
       x = "Mean life expectancy at birth, 2020 (years)", y = NULL,
       caption = str_wrap("Source: World Bank, World Development Indicators (WDI), 2020. Unweighted mean across countries.", 95))

mean(c(70, 75, NA))
mean(c(70, 75, NA), na.rm = TRUE)   # ignore the NAs

wdi |> summarise(across(everything(), \(x) sum(is.na(x))))

wdi |>
  filter(year == 2020) |>
  summarise(
    n_countries   = n(),
    n_with_gdp_pc = sum(!is.na(gdp_pc_ppp)),
    n_with_gini   = sum(!is.na(gini)),
    median_gdp_pc = median(gdp_pc_ppp, na.rm = TRUE),
    median_gini   = median(gini, na.rm = TRUE)
  )

gdp_n_region <- wdi |>
  filter(year == 2020) |>
  group_by(region) |>
  summarise(
    n_countries   = n(),
    n_gdp_pc      = sum(!is.na(gdp_pc_ppp)),
    median_gdp_pc = median(gdp_pc_ppp, na.rm = TRUE),
    .groups = "drop"
  ) |>
  arrange(desc(median_gdp_pc))

gdp_n_region

gdp_n_region |>
  mutate(lost = n_countries - n_gdp_pc) |>
  select(region, `Kept` = n_gdp_pc, `Lost (NA)` = lost) |>
  pivot_longer(-region, names_to = "status", values_to = "countries") |>
  mutate(status = factor(status, levels = c("Lost (NA)", "Kept"))) |>
  ggplot(aes(x = countries, y = fct_reorder(region, countries, sum), fill = status)) +
  geom_col() +
  scale_fill_manual(values = c(Kept = pal_cours[["principal"]], `Lost (NA)` = pal_cours[["accent"]]),
                    labels = c(Kept = "Kept", `Lost (NA)` = "Lost (NA)")) +
  labs(title = "Dropping NAs removes countries unevenly across regions",
       x = "Number of countries, 2020", y = NULL, fill = NULL,
       caption = "Source: World Bank, World Development Indicators (WDI), 2020. GDP per capita (PPP).") +
  theme(legend.position = "top")

wdi |>
  filter(year == 2020) |>
  summarise(mean_life = mean(life_exp, na.rm = TRUE), n = n(), .by = income)

wdi |> filter(year == 2020) |> count(income, sort = TRUE)

unique(wdi$income)   # check the exact labels before using them

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

wdi <- wdi |>
  mutate(income = factor(income, levels = income_levels))

wdi |> filter(year == 2020) |> count(income)

wdi |>
  filter(year == 2020) |>
  summarise(median_life = median(life_exp, na.rm = TRUE), .by = region) |>
  ggplot(aes(x = median_life, y = fct_reorder(region, median_life))) +
  geom_col(fill = pal_cours[["principal"]]) +
  labs(title = str_wrap("Life expectancy is lowest in Sub-Saharan Africa", 60),
       x = "Median life expectancy at birth, 2020 (years)",
       y = NULL,
       caption = str_wrap("Source: World Bank, WDI. Median across countries of each region.", 95))

life_long <- wdi |>
  filter(country %in% c("Algeria", "Brazil", "Ghana", "India", "Norway"),
         year %in% c(2000, 2010, 2020)) |>
  select(country, year, life_exp)

life_long

life_wide <- life_long |>
  pivot_wider(names_from = year, values_from = life_exp, names_prefix = "y")

life_wide

life_back <- life_wide |>
  pivot_longer(cols = starts_with("y"),
               names_to = "year", names_prefix = "y",
               values_to = "life_exp") |>
  mutate(year = as.integer(year))   # year came back as text

life_back

# Check: we recover exactly the table we started from
all.equal(as.data.frame(life_back), as.data.frame(life_long))

data1 <- tibble(country = c("US", "CN", "IN"), value1 = 1:3)
data2 <- tibble(country = c("US", "CN", "BR"), value2 = 4:6)

data1
data2

left_join(data1, data2, by = "country")    # all rows of data1
inner_join(data1, data2, by = "country")   # only rows found in both
full_join(data1, data2, by = "country")    # all rows of both
anti_join(data1, data2, by = "country")    # rows of data1 with NO match in data2

gdp3 <- read_csv(here("data", "wdi_3countries.csv"), show_col_types = FALSE)
wid  <- read_csv(here("data", "wid_gini.csv"), show_col_types = FALSE)

# first row of each country
gdp3 |> arrange(country, year) |> distinct(country, .keep_all = TRUE)
wid  |> arrange(country_wid, year) |> distinct(country_wid, .keep_all = TRUE)

naive <- gdp3 |>
  left_join(wid, by = c("country" = "country_wid", "year"))

nrow(gdp3)
nrow(naive)

naive |>
  summarise(n_rows = n(), n_gini = sum(!is.na(gini)), .by = country)

naive |>
  summarise(n_gini = sum(!is.na(gini)), n_missing = sum(is.na(gini)), .by = country) |>
  pivot_longer(c(n_gini, n_missing), names_to = "status", values_to = "years") |>
  mutate(status = recode(status, n_gini = "Gini found", n_missing = "Gini missing (NA)")) |>
  ggplot(aes(x = years, y = country, fill = status)) +
  geom_col() +
  scale_fill_manual(values = c(`Gini found` = pal_cours[["principal"]],
                               `Gini missing (NA)` = pal_cours[["accent"]]),
                    labels = c(`Gini found` = "Gini found", `Gini missing (NA)` = "Gini missing (NA)")) +
  labs(title = "The join on names silently lost the United States",
       x = "Number of years (2000-2020)", y = NULL, fill = NULL,
       caption = "Source: World Bank WDI and World Inequality Database (WID).") +
  theme(legend.position = "top")

# 1. anti_join: which keys of the left table have no match on the right?
gdp3 |>
  anti_join(wid, by = c("country" = "country_wid", "year")) |>
  distinct(country)

# 2. count the NAs created by the join
sum(is.na(naive$gini))

gdp_gini <- gdp3 |>
  left_join(select(wid, iso3c, year, gini), by = c("iso3c", "year"))

nrow(gdp3)
nrow(gdp_gini)
sum(is.na(gdp_gini$gini))

gdp_gini |>
  summarise(n_rows = n(), n_gini = sum(!is.na(gini)), .by = country)

# And the anti_join check is now empty
gdp3 |> anti_join(wid, by = c("iso3c", "year")) |> nrow()

gdp_gini |>
  ggplot(aes(x = year, y = gini, colour = country)) +
  geom_line(linewidth = 1) +
  scale_colour_manual(values = c(China = pal_cours[["accent"]],
                                 India = "#E69F00",
                                 `United States` = pal_cours[["principal"]])) +
  labs(title = str_wrap("Inequality rose in all three countries; in China, mostly before 2010", 60),
       x = "Year", y = "Gini coefficient (WID, 0-1)",
       colour = NULL,
       caption = str_wrap("Source: World Bank WDI (GDP) joined with World Inequality Database (WID, Gini) on iso3c and year.", 95))

# # Each (iso3c, year) must appear at most once in wid, otherwise: error
# gdp3 |>
#   left_join(wid, by = c("iso3c", "year"), relationship = "many-to-one")

gdp_per_capita <- function(gdp, population) {
  gdp / population
}

gdp_per_capita(gdp = 2e12, population = 50e6)   # 2,000 billion / 50 million = 40,000

gdp_pc_3 <- gdp3 |>
  filter(year == 2020) |>
  left_join(select(wdi, iso3c, year, population), by = c("iso3c", "year")) |>
  mutate(gdp_pc = gdp_per_capita(gdp, population))

nrow(gdp_pc_3)   # 3 rows, as in the left table
gdp_pc_3

growth_pct <- function(before, after) {
  100 * (after - before) / before
}

growth_pct(100, 110)   # check by hand: 10

growth_3 <- wdi |>
  filter(iso3c %in% c("CHN", "IND", "USA"), year %in% c(2000, 2020)) |>
  select(country, year, gdp_pc_ppp) |>
  pivot_wider(names_from = year, values_from = gdp_pc_ppp, names_prefix = "y") |>
  mutate(growth = growth_pct(y2000, y2020))

growth_3

ssa_2020 <- wdi |>
  filter(region == "Sub-Saharan Africa", year == 2020) |>
  select(country, income, life_exp, elec_access) |>
  arrange(life_exp)

nrow(ssa_2020)   # N
head(ssa_2020, 10)

wdi |>
  filter(year == 2020) |>
  mutate(gdp_bn = gdp / 1e9) |>
  arrange(desc(gdp_bn)) |>
  select(country, income, gdp_bn) |>
  head(5)

wdi |>
  filter(year %in% c(2000, 2010, 2020)) |>
  summarise(
    n_countries   = sum(!is.na(gdp_pc_ppp)),
    median_gdp_pc = median(gdp_pc_ppp, na.rm = TRUE),
    .by = c(region, year)
  ) |>
  arrange(year, desc(median_gdp_pc))

life_region <- wdi |>
  filter(year %in% c(2000, 2010, 2020)) |>
  summarise(median_life = median(life_exp, na.rm = TRUE), .by = c(region, year)) |>
  pivot_wider(names_from = year, values_from = median_life, names_prefix = "y")

life_region

life_region |>
  pivot_longer(starts_with("y"), names_to = "year", names_prefix = "y",
               values_to = "median_life") |>
  mutate(year = as.integer(year))   # year comes back as text

elec <- wdi |>
  filter(year == 2020) |>
  select(iso3c, year, elec_access)

small <- tibble(iso3c = c("FRA", "GHA", "XXX"))

joined <- small |> left_join(elec, by = "iso3c")

nrow(small)
nrow(joined)   # same: left_join keeps all rows of `small`
joined         # XXX has NA

small |> anti_join(elec, by = "iso3c")   # the code that failed to match

# 1. Read: every column as text, so that R guesses nothing behind your back
repair_raw <- read_csv(here("data", "raw", "repair_me.csv"),
                       col_types = cols(.default = col_character()))

# 2. Inspect
read_lines(here("data", "raw", "repair_me.csv"), n_max = 3)   # the file as stored
glimpse(repair_raw)
repair_raw |> distinct(country)
repair_raw |> count(country, year) |> filter(n > 1)
sort(unique(repair_raw$life_exp))

# 3. Repair
repaired <- repair_raw |>
  mutate(
    country    = str_to_title(str_trim(country)),     # " kenya", "KENYA " -> "Kenya"
    year       = as.integer(year),
    gdp_pc_ppp = parse_number(gdp_pc_ppp),            # "2,967.7" -> 2967.7
    life_exp   = as.numeric(na_if(life_exp, "n/a"))   # "n/a" -> NA, then a number
  ) |>                                                # iso_num stays text: "004"
  distinct()                                          # drops the second Ghana 2019 row

# 4. Check the repair
glimpse(repaired)
nrow(repaired)                                      # 12 country-years
repaired |> count(country, year) |> filter(n > 1)   # empty: no duplicate left

# 5. Mean life expectancy by year, with N
life_by_year <- repaired |>
  summarise(
    n_countries = n(),
    n_life      = sum(!is.na(life_exp)),
    mean_life   = mean(life_exp, na.rm = TRUE),
    .by = year
  )

life_by_year

# Same five countries in both years: drop Benin, missing in 2015
life_same <- repaired |>
  filter(country != "Benin") |>
  summarise(n_countries = n(), mean_life = mean(life_exp), .by = year)

life_same