# Séance 4 : Statistiques descriptives et régression
# 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")
library(tidyverse)
library(here)
library(scales)
library(gt)
library(fixest)
library(modelsummary)

source(here::here("R", "theme_cours.R"))
set.seed(2024)   # le seul tirage aléatoire de la page : nuages de points simulés

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

wdi |>
  filter(year %in% c(2000, 2010, 2015, 2019, 2021, 2023)) |>
  group_by(year) |>
  summarise(
    countries       = n(),
    `life exp.`     = sum(!is.na(life_exp)),
    `GDP pc`        = sum(!is.na(gdp_pc_ppp)),
    `both`          = sum(!is.na(life_exp) & !is.na(gdp_pc_ppp)),
    `+ Gini`        = sum(!is.na(life_exp) & !is.na(gdp_pc_ppp) & !is.na(gini)),
    `electricity`   = sum(!is.na(elec_access))
  ) |>
  gt() |>
  tab_header(title = "Nombre de pays avec données non manquantes") |>
  tab_source_note("Source : Banque mondiale, World Development Indicators (WDI). Pays uniquement (agrégats exclus).")

wdi |> filter(iso3c == "CAF") |> select(year, life_exp) |> print(n = Inf)

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

cs <- wdi |>
  filter(year == 2019, !is.na(life_exp), !is.na(gdp_pc_ppp)) |>
  mutate(
    income = factor(income, levels = c("Low income", "Lower middle income",
                                       "Upper middle income", "High income")),
    log_gdp_pc = log(gdp_pc_ppp)
  )

n_cs <- nrow(cs)
n_cs

mean(cs$life_exp)
median(cs$life_exp)
sd(cs$life_exp)
quantile(cs$life_exp, probs = c(0.1, 0.25, 0.5, 0.75, 0.9))
IQR(cs$life_exp)

summary(cs$life_exp)
summary(cs$gdp_pc_ppp)

mm <- tibble(
  statistic = c("Mean", "Median"),
  value     = c(mean(cs$gdp_pc_ppp), median(cs$gdp_pc_ppp))
)

ggplot(cs, aes(gdp_pc_ppp)) +
  geom_histogram(bins = 30, fill = "grey80", colour = "white") +
  geom_vline(data = mm, aes(xintercept = value, colour = statistic), linewidth = 1.1) +
  scale_colour_manual(values = c(Mean = pal_cours[["accent"]], Median = pal_cours[["principal"]]),
                      labels = c(Mean = "Moyenne", Median = "Médiane")) +
  scale_x_continuous(labels = label_dollar()) +
  labs(
    title = "Quelques pays riches tirent la moyenne au-dessus de la médiane",
    x = "PIB par habitant (en $ PPA constants de 2021)", y = "Nombre de pays", colour = NULL,
    caption = paste0("Source : Banque mondiale, World Development Indicators (WDI), 2019. N = ", n_cs, " pays.")
  ) +
  theme(legend.position = "top")

cs |>
  select(gdp_pc_ppp, log_gdp_pc) |>
  pivot_longer(everything(), names_to = "variable", values_to = "value") |>
  # factor() : fixer l'ordre des panneaux et donner à chaque variable un libellé lisible
  mutate(variable = factor(variable,
                           levels = c("gdp_pc_ppp", "log_gdp_pc"),
                           labels = c("GDP per capita (2021 PPP $)", "Log of GDP per capita"))) |>
  ggplot(aes(value)) +
  geom_histogram(bins = 25, fill = pal_cours[["principal"]], colour = "white") +
  facet_wrap(~ variable, scales = "free_x") +
  scale_x_continuous(labels = label_number(big.mark = ",")) +
  labs(
    title = "La plupart des pays sont pauvres ; quelques-uns sont très riches",
    x = NULL, y = "Nombre de pays",
    caption = paste0("Source : Banque mondiale, World Development Indicators (WDI), 2019. N = ", n_cs, " pays.")
  )

by_inc <- cs |>
  group_by(income) |>
  summarise(
    N          = n(),
    mean_le    = mean(life_exp),
    median_le  = median(life_exp),
    sd_le      = sd(life_exp),
    iqr_le     = IQR(life_exp)
  )
by_inc

ggplot(cs, aes(income, life_exp)) +
  geom_boxplot(fill = "grey90", colour = pal_cours[["principal"]], outlier.shape = NA) +
  geom_jitter(width = 0.15, alpha = 0.5, colour = pal_cours[["principal"]]) +
  labs(
    title = str_wrap("Les pays plus riches vivent plus longtemps ; la dispersion est proche d'un groupe à l'autre", 60),
    x = NULL, y = "Espérance de vie à la naissance (ans)",
    caption = str_wrap(paste0("Source : Banque mondiale, World Development Indicators (WDI), 2019. N = ", n_cs, " pays."), 95)
  )

cs |>
  group_by(region) |>
  summarise(
    N               = n(),
    gdp_mean        = mean(gdp_pc_ppp),
    gdp_median      = median(gdp_pc_ppp),
    le_mean         = mean(life_exp),
    le_sd           = sd(life_exp)
  ) |>
  arrange(desc(le_mean)) |>
  gt() |>
  cols_label(
    region     = "Région",
    N          = "N",
    gdp_mean   = "Moyenne",
    gdp_median = "Médiane",
    le_mean    = "Moyenne",
    le_sd      = "Écart-type"
  ) |>
  tab_spanner("PIB par habitant ($ PPA 2021)", c(gdp_mean, gdp_median)) |>
  tab_spanner("Espérance de vie (ans)", c(le_mean, le_sd)) |>
  fmt_currency(c(gdp_mean, gdp_median), decimals = 0) |>
  fmt_number(c(le_mean, le_sd), decimals = 1) |>
  tab_header(title = "Revenu et espérance de vie par région, 2019") |>
  tab_source_note("Source : Banque mondiale, World Development Indicators (WDI). Une observation par pays ; N = nombre de pays. Les régions sont les catégories de la Banque mondiale.")

elec <- cs |> filter(!is.na(elec_access))
nrow(cs) - nrow(elec)   # pays écartés

elec |>
  group_by(income) |>
  summarise(N = n(), mean = mean(elec_access), median = median(elec_access)) |>
  gt() |>
  fmt_number(c(mean, median), decimals = 1) |>
  tab_header(title = "Accès à l'électricité (% de la population), 2019") |>
  tab_source_note("Source : Banque mondiale, World Development Indicators (WDI). N = nombre de pays avec données.")

low <- elec |> filter(income == "Low income")
c(mean = mean(low$elec_access), median = median(low$elec_access))

n_sim <- 200
x <- rnorm(n_sim)
make_y <- function(r, x) r * x + sqrt(1 - r^2) * rnorm(length(x))

x_u <- runif(n_sim, -2, 2)
sim <- bind_rows(
  tibble(x = x,   y = make_y(0.9, x),  panel = "Strong positive"),
  tibble(x = x,   y = make_y(0, x),    panel = "No relationship"),
  tibble(x = x,   y = make_y(-0.9, x), panel = "Strong negative"),
  tibble(x = x_u, y = x_u^2,           panel = "Strong but non-linear")
) |>
  group_by(panel) |>
  mutate(r = cor(x, y)) |>
  ungroup() |>
  # fct_inorder() : garder les panneaux dans leur ordre d'apparition, pas alphabétique
  mutate(panel = fct_inorder(paste0(panel, " (r = ", sprintf("%.2f", r), ")")))

ggplot(sim, aes(x, y)) +
  geom_point(alpha = 0.6, colour = pal_cours[["principal"]]) +
  facet_wrap(~ panel, scales = "free") +
  labs(
    title = str_wrap("r mesure l'association linéaire, pas l'existence d'un lien entre deux variables", 60),
    x = "X", y = "Y",
    caption = str_wrap("Source : données simulées.", 95)
  )

r_lin <- cor(cs$gdp_pc_ppp, cs$life_exp, use = "complete.obs")
r_log <- cor(cs$log_gdp_pc, cs$life_exp, use = "complete.obs")
c(levels = r_lin, logs = r_log)

ct <- cor.test(cs$log_gdp_pc, cs$life_exp)
ct

cor(cs$elec_access, cs$life_exp, use = "complete.obs")
cor.test(cs$elec_access, cs$life_exp)

cs |>
  filter(!is.na(elec_access)) |>
  ggplot(aes(elec_access, life_exp)) +
  geom_point(alpha = 0.6, colour = pal_cours[["principal"]]) +
  labs(
    title = str_wrap("L'espérance de vie augmente avec l'accès à l'électricité, mais le lien s'infléchit", 60),
    x = "Accès à l'électricité (% de la population)",
    y = "Espérance de vie à la naissance (ans)",
    caption = str_wrap("Source : Banque mondiale, World Development Indicators (WDI), 2019.", 95)
  )

m1 <- feols(life_exp ~ log_gdp_pc, data = cs, vcov = "hetero")
summary(m1)

b1 <- coef(m1)[["log_gdp_pc"]]

cs_fit <- cs |>
  mutate(fitted_le = fitted(m1)) |>
  arrange(gdp_pc_ppp)

# dix pays répartis le long de la distribution des revenus
shown <- cs_fit |> slice(round(seq(5, nrow(cs_fit) - 4, length.out = 10)))

ggplot(cs_fit, aes(gdp_pc_ppp, life_exp)) +
  geom_point(colour = "grey70", alpha = 0.7) +
  geom_line(aes(y = fitted_le), colour = pal_cours[["principal"]], linewidth = 1.1) +
  geom_segment(data = shown, aes(xend = gdp_pc_ppp, yend = fitted_le),
               colour = pal_cours[["accent"]], linewidth = 0.9) +
  geom_point(data = shown, colour = pal_cours[["accent"]], size = 2.5) +
  scale_x_log10(labels = label_dollar()) +
  labs(
    title = "Un résidu est l'écart entre un pays et la droite",
    x = "PIB par habitant ($ PPA 2021, échelle log)",
    y = "Espérance de vie à la naissance (ans)",
    caption = paste0("Source : Banque mondiale, World Development Indicators (WDI), 2019. N = ", n_cs, " pays.")
  )

ggplot(cs, aes(gdp_pc_ppp, life_exp)) +
  geom_point(alpha = 0.6, colour = pal_cours[["principal"]]) +
  geom_smooth(method = "lm", formula = y ~ x, colour = pal_cours[["accent"]],
              fill = pal_cours[["accent"]], alpha = 0.15) +
  scale_x_log10(labels = label_dollar()) +
  labs(
    title = "Les pays plus riches ont une espérance de vie plus longue",
    x = "PIB par habitant ($ PPA 2021, échelle log)",
    y = "Espérance de vie à la naissance (ans)",
    caption = paste0("Source : Banque mondiale, World Development Indicators (WDI), 2019. N = ", n_cs, " pays.")
  )

m_el <- feols(life_exp ~ elec_access, data = cs, vcov = "hetero")
summary(m_el)
r2(m_el, "r2")

n_gini <- cs |> filter(!is.na(gini)) |> nrow()
tibble(sample = c("All countries", "With Gini"), N = c(n_cs, n_gini))

cs_g <- cs |>
  filter(!is.na(gini)) |>
  mutate(region = relevel(factor(region), ref = "Sub-Saharan Africa"))
cs <- cs |>
  mutate(region = relevel(factor(region), ref = "Sub-Saharan Africa"))

m2 <- feols(life_exp ~ log_gdp_pc,               data = cs_g, vcov = "hetero")
m3 <- feols(life_exp ~ log_gdp_pc + gini,        data = cs_g, vcov = "hetero")
m4 <- feols(life_exp ~ log_gdp_pc + gini + region, data = cs_g, vcov = "hetero")

nice_names <- function(x) {
  x |>
    str_replace("^log_gdp_pc$", "Log du PIB par habitant") |>
    str_replace("^gini$", "Indice de Gini (0-100)") |>
    str_replace("^region", "Région : ") |>
    str_replace("^\\(Intercept\\)$", "Constante")
}

modelsummary(
  output = "gt",   # un tableau gt : même aspect que les autres tableaux, sans script externe
  fmt = \(x) nb(x, 3),   # séparateur décimal de la langue de la page
  list("(1) Tous" = m1, "(2) Avec Gini" = m2, "(3) + Gini" = m3, "(4) + Région" = m4),
  coef_rename = nice_names,
  stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
  gof_map = c("nobs", "r.squared"),
  title = "Espérance de vie, coupe transversale 2019 (MCO)",
  notes = c("Variable dépendante : espérance de vie à la naissance (ans). Erreurs-types robustes à l'hétéroscédasticité entre parenthèses.",
            "Catégorie de référence des régions : Afrique subsaharienne.",
            "Source : Banque mondiale, World Development Indicators (WDI), 2019. Les colonnes (2) à (4) n'utilisent que les pays disposant d'un indice de Gini.")
)

gdp3    <- coef(m3)[["log_gdp_pc"]]
se_gdp3 <- se(m3)[["log_gdp_pc"]]
b <- round(c(coef(m1)[["log_gdp_pc"]], coef(m2)[["log_gdp_pc"]],
             coef(m3)[["log_gdp_pc"]], coef(m4)[["log_gdp_pc"]]), 2)
bg  <- coef(m3)[["gini"]]
bg4 <- coef(m4)[["gini"]]
se3 <- se(m3)[["gini"]]
p3  <- pvalue(m3)[["gini"]]
p4  <- pvalue(m4)[["gini"]]

modelplot(
  list("(1) Tous" = m1, "(2) Avec Gini" = m2, "(3) + Gini" = m3, "(4) + Région" = m4),
  coef_rename = nice_names,
  coef_omit = "Intercept|region"
) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey40") +
  scale_colour_manual(values = unname(c(pal_cours[["gris"]], pal_cours[["principal"]],
                                        pal_regions[["Sub-Saharan Africa"]], pal_cours[["accent"]]))) +
  labs(
    title = "Le gradient de revenu est stable ; le coefficient de Gini est proche de zéro",
    x = "Coefficient (ans d'espérance de vie par unité), IC à 95 %", y = NULL, colour = NULL,
    caption = "Source : Banque mondiale, World Development Indicators (WDI), 2019. Colonnes (2)-(4) : pays ayant un indice de Gini."
  ) +
  theme(legend.position = "bottom")

cs_el  <- cs |> filter(!is.na(elec_access))   # même échantillon pour les deux modèles
m1_el  <- feols(life_exp ~ log_gdp_pc,               data = cs_el, vcov = "hetero")
m_ctrl <- feols(life_exp ~ log_gdp_pc + elec_access, data = cs_el, vcov = "hetero")
etable(m1_el, m_ctrl)

cs_t <- cs |>
  mutate(group = if_else(region == "Sub-Saharan Africa", "Sub-Saharan Africa", "Rest of the world"))

cs_t |>
  group_by(group) |>
  summarise(N = n(), mean = mean(life_exp), sd = sd(life_exp))

tt <- t.test(life_exp ~ group, data = cs_t)   # Welch par défaut
tt

gap <- unname(tt$estimate[1] - tt$estimate[2])

# (a) revenu élevé contre le reste
cs_h <- cs |> mutate(group = if_else(income == "High income", "High income", "Other"))
count(cs_h, group)
t.test(life_exp ~ group, data = cs_h)

# (b) deux modèles dans un seul tableau
ma <- feols(life_exp ~ log_gdp_pc,              data = cs, vcov = "hetero")
mb <- feols(life_exp ~ log_gdp_pc + elec_access, data = cs, vcov = "hetero")
modelsummary(
  output = "gt",   # un tableau gt : même aspect que les autres tableaux, sans script externe
  fmt = \(x) nb(x, 3),   # séparateur décimal de la langue de la page
  list("(1)" = ma, "(2)" = mb),
  coef_rename = c("log_gdp_pc" = "Log du PIB par habitant",
                  "elec_access" = "Accès à l'électricité (%)",
                  "(Intercept)" = "Constante"),
  stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
  gof_map = c("nobs", "r.squared"),
  notes = c("Variable dépendante : espérance de vie (ans). Erreurs-types robustes entre parenthèses.",
            "* p < 0,1, ** p < 0,05, *** p < 0,01. Source : Banque mondiale, World Development Indicators (WDI), 2019.")
)