Exercise 1 (20 minutes)

You practice: multiple OLS — adding a control that changes the story. Back to the carbon divide: do regions differ in emissions beyond their wealth? When you're done, assist some of the others …

1. Do the Setup if needed, then copy in this starter code — the same carbon-divide pipeline as Lecture 4 (World Bank CO2, GDP, and region), with Palau dropped and Sub-Saharan Africa set as the reference region:

pacman::p_load(tidyverse, wbstats, estimatr, modelsummary)

co2_raw <- tryCatch(
  wb_data("EN.GHG.CO2.PC.CE.AR5", start_date = 2000, end_date = 2023),
  error = function(e) readRDS("data/wb_co2_raw.rds"))
gdp_raw <- tryCatch(
  wb_data("NY.GDP.PCAP.PP.KD", start_date = 2000, end_date = 2023),
  error = function(e) readRDS("data/wb_gdp_raw.rds"))
countries <- tryCatch(
  wb_countries(), error = function(e) readRDS("data/wb_countries_raw.rds"))

Dat_co2 <- co2_raw %>%
  rename(co2 = EN.GHG.CO2.PC.CE.AR5) %>%
  select(iso3c, country, year = date, co2) %>%
  drop_na(co2) %>% group_by(country) %>% filter(year == max(year)) %>% ungroup()
Dat_gdp <- gdp_raw %>%
  rename(gdp = NY.GDP.PCAP.PP.KD) %>%
  select(iso3c, year = date, gdp) %>% drop_na(gdp) %>%
  group_by(iso3c) %>% filter(year == max(year)) %>% ungroup() %>% select(iso3c, gdp)
Dat_meta <- countries %>%
  filter(region != "Aggregates") %>% select(iso3c, region, income_level)

Dat <- Dat_co2 %>%
  inner_join(Dat_gdp,  by = "iso3c") %>%
  inner_join(Dat_meta, by = "iso3c") %>%
  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: Dat <- readRDS(url("https://willowy-quokka-e604ee.netlify.app/data/Dat_L4.rds")) |> dplyr::filter(country != "Palau")
• Still stuck? Ask in the t-R-ouble forum on Absalon.

2. Write it yourself. Fit two models and put them side by side: a bivariate one, co2 ~ region, and an adjusted one that adds wealth, co2 ~ region + gdp_k.

Hint

Fit both with lm_robust(), then pass them as a named list() to modelsummary().

R solution -> don't peek too early ;) !

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.911***	0.165
	(0.218)	(0.256)
regionEast Asia & Pacific	3.315***	0.799
	(0.943)	(0.856)
regionEurope & Central Asia	4.679***	-0.141
	(0.471)	(1.471)
regionLatin America & Caribbean	2.065***	-0.171
	(0.574)	(0.820)
regionMiddle East, North Africa, Afghanistan & Pakistan	8.728**	5.803***
	(2.618)	(1.610)
regionNorth America	9.711**	1.001
	(3.325)	(5.405)
regionSouth Asia	0.703	-0.008
	(0.499)	(0.352)
gdp_k		0.121***
		(0.035)
Num.Obs.	190	190
R2	0.241	0.465

3. Compare the region coefficients across the two models. What happens when GDP enters? They barely changeThey shrink sharply toward zeroThey all grow largerThey change sign

4. 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.

5. Write one sentence (as a # comment) interpreting that change, then compare:

Model answer

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.

6. So how is GDP best understood in the link between region and emissions? A confounder we must control to be fairA mediator: richer regions emit more largely because they are rich — wealth is the channelNeither; it is irrelevant here

---

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?

R solution

lm_robust(co2 ~ region + gdp_k, data = Dat) %>%
  tidy() %>%
  filter(str_detect(term, "region")) %>%
  arrange(desc(estimate)) %>%
  select(term, estimate)
#                                                      term estimate
# 1 regionMiddle East, North Africa, Afghanistan & Pakistan   5.8033
# 2                                     regionNorth America   1.0012
# 3                               regionEast Asia & Pacific   0.7993
# 4                                        regionSouth Asia  -0.0083
# 5                             regionEurope & Central Asia  -0.1414
# 6                         regionLatin America & Caribbean  -0.1714

The Middle East & North Africa region stays about 5.8 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.