# Séance 3 : Manipuler les données (dplyr, tidyr, jointures)
# Source : https://pierrebeaucoral.github.io/r-for-beginners/fr/
# Ouvrez votre projet r-course dans RStudio, puis exécutez ce fichier ligne par ligne.

Sys.setenv(LANGUAGE = "en")   # messages de R en anglais
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)
# L'espérance de vie de la RCA n'est pas fiable dans cette version des WDI (14,7 en 2009, 18,8 en 2022) : on l'exclut pour toutes les années
wdi <- wdi |> mutate(life_exp = if_else(iso3c == "CAF", NA_real_, life_exp))

glimpse(wdi)

# Sans le pipe : on lit de l'intérieur vers l'extérieur
head(arrange(filter(wdi, year == 2020), desc(life_exp)), 3)

# Avec le pipe : on lit de haut en bas
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)

# Aides : sélectionner des colonnes par motif ou par 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() = du plus grand au plus petit
  select(country, gdp_pc_ppp) |>
  head(5)

wdi |>
  filter(year == 2020) |>
  mutate(gdp_bn = gdp / 1e9,                      # PIB en milliards de US$
         pop_m  = population / 1e6) |>            # population en 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 = "La plupart des pays ont un accès quasi universel à l'électricité",
       x = "Nombre de pays, 2020", y = NULL,
       caption = "Source : Banque mondiale, 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()                                    # départ : tous les pays-années
wdi |> filter(year == 2020) |> nrow()            # une année : une ligne par pays
wdi |> filter(year == 2020, region == "Sub-Saharan africa") |> nrow()   # 0 : le problème est ici

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

# Bogue 1 (erreur dans select) : la colonne s'appelle life_exp, pas life_expectancy
names(wdi)

# Bogue 2 (silencieux) : aucune ligne ne survit au filtre
wdi |> filter(year == 2019, income == "Low Income") |> nrow()
wdi |> distinct(income)   # le libellé est "Low income", avec un i minuscule

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("L'Afrique subsaharienne accuse 17 ans de retard sur l'Amérique du Nord", 60),
       x = "Espérance de vie moyenne, 2020 (ans)", y = NULL,
       caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2020. Moyenne non pondérée entre pays.", 95))

mean(c(70, 75, NA))
mean(c(70, 75, NA), na.rm = TRUE)   # ignorer les NA

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 = "Conservés", `Lost (NA)` = "Perdus (NA)")) +
  labs(title = "Supprimer les NA retire des pays de façon inégale selon les régions",
       x = "Nombre de pays, 2020", y = NULL, fill = NULL,
       caption = "Source : Banque mondiale, World Development Indicators (WDI), 2020. PIB par habitant (PPA).") +
  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)   # vérifier les libellés exacts avant de les utiliser

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("L'espérance de vie est la plus faible en Afrique subsaharienne", 60),
       x = "Espérance de vie médiane, 2020 (ans)",
       y = NULL,
       caption = str_wrap("Source : Banque mondiale, WDI. Médiane entre les pays de chaque région.", 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 est revenue sous forme de texte

life_back

# Vérification : on retrouve exactement la table de départ
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")    # toutes les lignes de data1
inner_join(data1, data2, by = "country")   # seulement les lignes présentes dans les deux
full_join(data1, data2, by = "country")    # toutes les lignes des deux
anti_join(data1, data2, by = "country")    # lignes de data1 SANS correspondance dans 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)

# première ligne de chaque pays
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 disponible", `Gini missing (NA)` = "Gini manquant (NA)")) +
  labs(title = "La jointure sur les noms a perdu les États-Unis en silence",
       x = "Nombre d'années (2000-2020)", y = NULL, fill = NULL,
       caption = "Source : Banque mondiale (WDI) et World Inequality Database (WID).") +
  theme(legend.position = "top")

# 1. anti_join : quelles clés de la table de gauche n'ont pas de correspondance à droite ?
gdp3 |>
  anti_join(wid, by = c("country" = "country_wid", "year")) |>
  distinct(country)

# 2. compter les NA créés par la jointure
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)

# Et la vérification par anti_join est maintenant vide
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("Les inégalités ont augmenté dans les trois pays ; en Chine, surtout avant 2010", 60),
       x = "Année", y = "Coefficient de Gini (WID, 0-1)",
       colour = NULL,
       caption = str_wrap("Source : Banque mondiale (WDI, PIB) jointe à la World Inequality Database (WID, Gini) par iso3c et année.", 95))

# # Chaque (iso3c, year) doit apparaître au plus une fois dans wid, sinon : erreur
# 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 milliards / 50 millions = 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 lignes, comme dans la table de gauche
gdp_pc_3

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

growth_pct(100, 110)   # vérification à la main : 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 revient sous forme de texte

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)   # identique : left_join garde toutes les lignes de `small`
joined         # XXX a des NA

small |> anti_join(elec, by = "iso3c")   # le code sans correspondance

# 1. Lire : toutes les colonnes comme du texte, pour que R ne devine rien dans votre dos
repair_raw <- read_csv(here("data", "raw", "repair_me.csv"),
                       col_types = cols(.default = col_character()))

# 2. Inspecter
read_lines(here("data", "raw", "repair_me.csv"), n_max = 3)   # le fichier tel qu'il est stocké
glimpse(repair_raw)
repair_raw |> distinct(country)
repair_raw |> count(country, year) |> filter(n > 1)
sort(unique(repair_raw$life_exp))

# 3. Réparer
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, puis un nombre
  ) |>                                                # iso_num reste du texte : "004"
  distinct()                                          # supprime la seconde ligne Ghana 2019

# 4. Vérifier la réparation
glimpse(repaired)
nrow(repaired)                                      # 12 pays-années
repaired |> count(country, year) |> filter(n > 1)   # vide : plus aucun doublon

# 5. Espérance de vie moyenne par année, avec 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

# Les cinq mêmes pays les deux années : on retire le Bénin, manquant en 2015
life_same <- repaired |>
  filter(country != "Benin") |>
  summarise(n_countries = n(), mean_life = mean(life_exp), .by = year)

life_same