Exercise 2 (20 minutes)

(Keep working in the R script from Exercise 1 — you already have Dat, ols_biv, and ols_adj. If not, the starter code and the "Stuck?" box are in Exercise 1.)

Write it yourself: a before/after coefficient plot. Tidy both models, keep only the region terms, stack them, and plot the two sets of region coefficients side by side with position_dodge() — so you can see each gap shrink when GDP enters.

tidy() each model, filter(str_detect(term, "region")), then bind_rows(bivariate, adjusted, .id = "model"). Map color = model and add position = position_dodge(width = 0.5) to geom_pointrange().

tidy_reg <- function(m) {
  m %>% tidy() %>%
    filter(str_detect(term, "region")) %>%
    mutate(term = str_remove(term, "region"))
}

plotdata <- bind_rows(tidy_reg(ols_biv), tidy_reg(ols_adj), .id = "model") %>%
  mutate(model = if_else(model == "1", "Bivariate", "Adjusted for GDP"))

ggplot(plotdata, aes(x = reorder(term, estimate), y = estimate, color = model)) +
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.5)) +
  coord_flip() +
  labs(x = "", y = "Difference in CO2 vs. Sub-Saharan Africa (t per person)",
       color = "") +
  theme_minimal() +
  theme(legend.position = "bottom")

Which statement matches your plot?
Write it yourself: a prediction plot. From ols_adj, predict CO2 across a GDP range (say 1–130) holding region at the reference, Sub-Saharan Africa, and plot the fitted line with its confidence ribbon.

Build tibble(gdp_k = 1:130, region = "Sub-Saharan Africa"), then predict(ols_adj, newdata = ..., interval = "confidence")$fit, as_tibble() it and bind_cols() back. Plot with geom_ribbon() + geom_line().

fict <- tibble(gdp_k = 1:130, region = "Sub-Saharan Africa")

fict <- predict(ols_adj, newdata = fict,
                interval = "confidence", level = 0.95)$fit %>%
  as_tibble() %>% bind_cols(fict, .)

ggplot(fict, aes(x = gdp_k, y = fit)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "GDP per person (thousands of $, PPP)\n(region held at Sub-Saharan Africa)",
       y = "Predicted tonnes of CO2 per person") +
  theme_minimal()

Write one sentence (as a # comment) describing what the prediction line shows, then compare:

Within a given region, predicted emissions rise steadily with national wealth — about 0.12 extra tonnes per person for every additional $1,000 of GDP per person — which is why controlling for wealth absorbs so much of the raw regional gaps.

Discuss with your neighbour — three steps:
1. Your two figures make the same point differently. Which better tells a non-statistician that "the regional gaps are mostly wealth" — the before/after coefficient plot or the prediction line?
2. Can we say wealth causes emissions here? Name one thing, besides GDP, that differs between rich and poor countries and could drive emissions.
3. Middle East & North Africa stays high even at equal GDP. What does a large surviving coefficient, after controlling, tell you — and what does a vanishing one tell you?

(a) The prediction line reads like a graph and needs no reference-category vocabulary; the coefficient plot is compact but assumes the reader knows what "vs. Sub-Saharan Africa" means. (b) No — these are descriptive cross-country comparisons; energy mix, industry structure, climate, and fossil-fuel endowments all differ and could drive emissions (all potential confounders). (c) A vanishing coefficient means the raw gap was "just" the control (wealth); a surviving one (MENA) means there is something regional beyond wealth — here, oil and gas.

Bonus (for the fast): swap the continuous control gdp_k for the categorical income_level. Do the region gaps shrink in the same way?

Dat <- Dat %>% mutate(income_level = fct_relevel(income_level, "Low income"))

modelsummary(
  list("CO2 (+ GDP)" = ols_adj,
       "CO2 (+ income group)" = lm_robust(co2 ~ region + income_level, data = Dat)),
  stars = TRUE, gof_map = c("nobs", "r.squared"))

	CO2 (+ GDP)	CO2 (+ income group)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	0.165	-0.272
	(0.256)	(0.325)
regionEast Asia & Pacific	0.799	-0.244
	(0.856)	(0.875)
regionEurope & Central Asia	-0.141	-1.172
	(1.471)	(1.132)
regionLatin America & Caribbean	-0.171	-2.513**
	(0.820)	(0.928)
regionMiddle East, North Africa, Afghanistan & Pakistan	5.803***	5.329**
	(1.610)	(1.828)
regionNorth America	1.001	2.369
	(5.405)	(3.564)
regionSouth Asia	-0.008	-0.238
	(0.352)	(0.755)
gdp_k	0.121***
	(0.035)
income_levelHigh income		8.525***
		(1.480)
income_levelLower middle income		0.960+
		(0.520)
income_levelUpper middle income		4.452***
		(0.879)
Num.Obs.	190	190
R2	0.465	0.434

Both controls tell the same story — wealth, measured either way, absorbs most of the regional differences. Good findings should not hinge on one particular way of measuring the control.

Exercise 2 (20 minutes)

You practice: visualizing the adjustment — a before/after coefficient plot and model predictions, on the carbon divide. When you're done, assist some of the others …