pacman::p_load(tidyverse, wbstats, estimatr, modelsummary)
wb_raw <- wb_data(c(co2 = "EN.GHG.CO2.PC.CE.AR5",
gdp = "NY.GDP.PCAP.PP.KD"),
start_date = 2000, end_date = 2026)
wb_meta <- wb_countries()
Dat <- wb_raw %>%
drop_na(co2, gdp) %>%
group_by(country) %>%
filter(date == max(date)) %>%
ungroup() %>%
inner_join(wb_meta %>% select(iso3c, region, income_level), by = "iso3c") %>%
select(iso3c, country, date, co2, gdp, region, income_level) %>%
filter(country != "Palau") %>% # the artefact from Lecture 4
mutate(gdp_k = gdp / 1000,
region = fct_relevel(region, "Sub-Saharan Africa"))
Stuck?
- Did
pacman::p_load(...)run without red errors? Fix that first.- Odd results? Session → Restart R, then rerun from the top.
- World Bank API down? Load the cached, ready-made data instead: download
Dat_L4.rdsinto your project folder, thenDat <- readRDS("Dat_L4.rds") %>% filter(country != "Palau")- Still stuck? Ask in the t-R-ouble forum on Absalon.
co2 ~ region,
and an adjusted one that adds wealth,
co2 ~ region + gdp_k.Fit both with lm_robust(), then pass them as a named
list() to modelsummary().
ols_biv <- lm_robust(co2 ~ region, data = Dat)
ols_adj <- lm_robust(co2 ~ region + gdp_k, data = Dat)
modelsummary(list("CO2 (bivariate)" = ols_biv, "CO2 (+ GDP)" = ols_adj),
stars = TRUE, gof_map = c("nobs", "r.squared"))
| CO2 (bivariate) | CO2 (+ GDP) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 0.922*** | 0.208 |
| (0.224) | (0.248) | |
| regionEast Asia & Pacific | 3.320*** | 0.866 |
| (0.926) | (0.820) | |
| regionEurope & Central Asia | 4.557*** | -0.064 |
| (0.463) | (1.379) | |
| regionLatin America & Caribbean | 2.109*** | -0.118 |
| (0.573) | (0.808) | |
| regionMiddle East, North Africa, Afghanistan & Pakistan | 8.617*** | 5.895*** |
| (2.571) | (1.657) | |
| regionNorth America | 9.667** | 1.360 |
| (3.236) | (5.194) | |
| regionSouth Asia | 0.756 | 0.039 |
| (0.517) | (0.368) | |
| gdp_k | 0.114*** | |
| (0.033) | ||
| Num.Obs. | 190 | 190 |
| R2 | 0.241 | 0.455 |
Compare the region coefficients across the two models. What happens when GDP enters?
Read off the North America coefficient in the adjusted model: compared with Sub-Saharan Africa, North America emits about tonnes more per person once wealth is held constant (two digits; point or comma is fine) — down from 9.7 t in the bivariate model.
Write one sentence (as a # comment)
interpreting that change, then compare:
Most of North America's raw emissions gap over Sub-Saharan Africa disappears once we compare countries of similar wealth — the regional gap largely reflects the fact that richer regions are rich, not something regional beyond wealth.
Bonus (for the fast): one region stays high even after controlling for GDP — it emits more than its wealth alone would predict. Which is it, and why?
lm_robust(co2 ~ region + gdp_k, data = Dat) %>%
tidy() %>%
filter(term != "(Intercept)" & term != "gdp_k") %>% # Keep only the region terms
select(term, estimate)
# term estimate
# 1 regionEast Asia & Pacific 0.866
# 2 regionEurope & Central Asia -0.064
# 3 regionLatin America & Caribbean -0.118
# 4 regionMiddle East, North Africa, Afghanistan & Pakistan 5.895
# 5 regionNorth America 1.360
# 6 regionSouth Asia 0.039
The Middle East & North Africa region stays about 5.9 t above the reference even at equal GDP — the oil-and-gas economies emit far more than their income alone implies. Controlling for wealth reveals which regions are genuinely emission-intensive.