# Session 4: Descriptive statistics and regression
# 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")
library(tidyverse)
library(here)
library(scales)
library(gt)
library(fixest)
library(modelsummary)

source(here::here("R", "theme_cours.R"))
set.seed(2024)   # the only random draw of the page: simulated scatterplots

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 = "Number of countries with non-missing data") |>
  tab_source_note("Source: World Bank, World Development Indicators (WDI). Countries only (aggregates removed).")

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

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

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 = "Mean", Median = "Median")) +
  scale_x_continuous(labels = label_dollar()) +
  labs(
    title = "A few rich countries drag the mean above the median",
    x = "GDP per capita (2021 PPP $)", y = "Number of countries", colour = NULL,
    caption = paste0("Source: World Bank, World Development Indicators (WDI), 2019. N = ", n_cs, " countries.")
  ) +
  theme(legend.position = "top")

cs |>
  select(gdp_pc_ppp, log_gdp_pc) |>
  pivot_longer(everything(), names_to = "variable", values_to = "value") |>
  # factor(): fix the panel order and give each variable a readable label
  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 = "Most countries are poor; a few are very rich",
    x = NULL, y = "Number of countries",
    caption = paste0("Source: World Bank, World Development Indicators (WDI), 2019. N = ", n_cs, " countries.")
  )

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("Richer countries live longer; the spread is similar across groups", 60),
    x = NULL, y = "Life expectancy at birth (years)",
    caption = str_wrap(paste0("Source: World Bank, World Development Indicators (WDI), 2019. N = ", n_cs, " countries."), 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     = "Region",
    N          = "N",
    gdp_mean   = "Mean",
    gdp_median = "Median",
    le_mean    = "Mean",
    le_sd      = "Std. dev."
  ) |>
  tab_spanner("GDP per capita (2021 PPP $)", c(gdp_mean, gdp_median)) |>
  tab_spanner("Life expectancy (years)", 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 = "Income and life expectancy by region, 2019") |>
  tab_source_note("Source: World Bank, World Development Indicators (WDI). One observation per country; N = number of countries. Regions are World Bank categories.")

elec <- cs |> filter(!is.na(elec_access))
nrow(cs) - nrow(elec)   # countries dropped

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 = "Access to electricity (% of population), 2019") |>
  tab_source_note("Source: World Bank, World Development Indicators (WDI). N = number of countries with data.")

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(): keep the panels in their order of appearance, not alphabetical
  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 measures linear association, not whether two variables are related", 60),
    x = "X", y = "Y",
    caption = str_wrap("Source: simulated data.", 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("Life expectancy rises with electricity access, but the link bends", 60),
    x = "Access to electricity (% of population)",
    y = "Life expectancy at birth (years)",
    caption = str_wrap("Source: World Bank, 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)

# ten countries spread along the income distribution
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 = "A residual is the gap between a country and the line",
    x = "GDP per capita (2021 PPP $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = paste0("Source: World Bank, World Development Indicators (WDI), 2019. N = ", n_cs, " countries.")
  )

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 = "Richer countries have longer life expectancy",
    x = "GDP per capita (2021 PPP $, log scale)",
    y = "Life expectancy at birth (years)",
    caption = paste0("Source: World Bank, World Development Indicators (WDI), 2019. N = ", n_cs, " countries.")
  )

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 GDP per capita") |>
    str_replace("^gini$", "Gini index (0-100)") |>
    str_replace("^region", "Region: ") |>
    str_replace("^\\(Intercept\\)$", "Constant")
}

modelsummary(
  output = "gt",   # a gt table: same look as the other tables, no external script
  fmt = \(x) nb(x, 3),   # decimal mark of the page language
  list("(1) All" = m1, "(2) With Gini" = m2, "(3) + Gini" = m3, "(4) + Region" = m4),
  coef_rename = nice_names,
  stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
  gof_map = c("nobs", "r.squared"),
  title = "Life expectancy, 2019 cross-section (OLS)",
  notes = c("Dependent variable: life expectancy at birth (years). Heteroskedasticity-robust standard errors in parentheses.",
            "Region reference category: Sub-Saharan Africa.",
            "Source: World Bank, World Development Indicators (WDI), 2019. Columns (2)-(4) use only countries with a Gini index.")
)

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) All" = m1, "(2) With Gini" = m2, "(3) + Gini" = m3, "(4) + Region" = 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 = "The income gradient is stable; the Gini coefficient is near zero",
    x = "Coefficient (years of life expectancy per unit), 95% CI", y = NULL, colour = NULL,
    caption = "Source: World Bank, World Development Indicators (WDI), 2019. Columns (2)-(4): countries with a Gini index."
  ) +
  theme(legend.position = "bottom")

cs_el  <- cs |> filter(!is.na(elec_access))   # same sample for both models
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 by default
tt

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

# (a) high-income vs the rest
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) two models in one table
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",   # a gt table: same look as the other tables, no external script
  fmt = \(x) nb(x, 3),   # decimal mark of the page language
  list("(1)" = ma, "(2)" = mb),
  coef_rename = c("log_gdp_pc" = "Log GDP per capita",
                  "elec_access" = "Electricity access (%)",
                  "(Intercept)" = "Constant"),
  stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
  gof_map = c("nobs", "r.squared"),
  notes = c("Dependent variable: life expectancy (years). Robust standard errors in parentheses.",
            "* p < 0.1, ** p < 0.05, *** p < 0.01. Source: World Bank, World Development Indicators (WDI), 2019.")
)