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

  1. Write it yourself: a coefficient plot. Fit co2 ~ income_level, tidy() it into a tibble, drop the intercept, and plot each coefficient with its 95% confidence interval (geom_pointrange() + coord_flip()).

tidy() gives you estimate, conf.low, and conf.high. Filter out "(Intercept)", then map estimate to y, the term to x, and add geom_hline(yintercept = 0, ...) as a reference line.

plotdata <- lm_robust(co2 ~ income_level, data = Dat) %>%
  tidy() %>%
  mutate(term = str_remove(term, "income_level")) %>%
  filter(term != "(Intercept)")

ggplot(data = plotdata,
       aes(y = estimate, x = reorder(term, estimate))) +
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  coord_flip() +
  labs(x = "Income group",
       y = "Average difference in CO2 (t per person)\nvs. low-income countries") +
  theme_minimal()

  1. Read your plot: compared with low-income countries, upper-middle-income countries emit about tonnes more COâ‚‚ per person (two digits).

  2. Which statement matches your coefficient plot?

  3. Write it yourself: a prediction plot. Fit the bivariate model co2 ~ state_ownership, predict across a synthetic range of the socialism index, and plot the fitted line with its confidence ribbon. Mark two benchmarks with dashed red lines: a market-based economy (-1) and a state-run one (+1).

The index runs from about −3 to +2, so build fict_dat <- tibble(state_ownership = seq(-3, 2, by = 0.1)), then predict(ols_state, newdata = fict_dat, interval = "confidence")$fit, turn it into a tibble and bind_cols() it back. Plot with geom_ribbon() + geom_line().

ols_state <- lm_robust(co2 ~ state_ownership, data = Dat)

# Fictional countries spanning the observed range of the index
fict_dat <- tibble(state_ownership = seq(-3, 2, by = 0.1)) #<<
fict_dat <- predict(ols_state, newdata = fict_dat,
                    interval = "confidence", level = 0.95)$fit %>%
  as_tibble() %>%
  bind_cols(fict_dat, .)

ggplot(data = fict_dat, aes(y = fit, x = state_ownership)) +
  geom_vline(xintercept = c(-1, 1), color = "#901A1E", lty = "dashed") + #<<
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "State ownership of the economy (V-Dem, higher = more)",
       y = "Predicted tonnes of CO2 per person") +
  theme_minimal()

  1. A prediction does not have to be a plot. You can also ask the model for a single number. Using the bivariate model co2 ~ state_ownership from Exercise 1, predict COâ‚‚ for a strongly market-based economy (state_ownership = -1, about where Denmark and Finland sit) and for a strongly state-run one (state_ownership = +1, about Russia and Angola). Build a tibble() with those two values and hand it to predict().

Same three steps as the plot, just with two rows instead of 120: tibble(state_ownership = c(-1, 1)), then predict(ols_state, newdata = ..., interval = "confidence")$fit.

ols_state <- lm_robust(co2 ~ state_ownership, data = Dat)

fict_two <- tibble(state_ownership = c(-1, 1)) # market-based vs state-run #<<

predict(ols_state, newdata = fict_two,
        interval = "confidence", level = 0.95)$fit %>%
  as_tibble() %>%
  bind_cols(fict_two, .)
# # A tibble: 2 × 4
#   state_ownership   fit   lwr   upr
#             <dbl> <dbl> <dbl> <dbl>
# 1              -1  3.84  3.17  4.51
# 2               1  5.90  3.30  8.50

The market-based country is predicted at tonnes per person, the state-run one at (two digits each).

  1. Look at the two confidence intervals you just produced. Which is wider, and why?

A prediction is only as good as the data behind it. Most countries in Dat sit on the market-based side, so the model has seen plenty of them and its interval there is narrow. Very few countries score high on state ownership, so out there the same model is guessing from a handful of cases — and says so, honestly, with a much wider interval. This is exactly the widening ribbon you drew in the prediction plot, read as numbers instead.

  1. Write one sentence (as a # comment) that a reader of your term paper could follow, describing what the prediction plot shows. Then compare:

Predicted CO₂ emissions rise with state ownership of the economy: a strongly market-based country is predicted to emit about 3.8 tonnes per person, a strongly state-run one about 5.9 — but the confidence ribbon fans out sharply on the right, where very few countries actually sit, so the model is far less sure about that end.


  1. Discuss with your neighbour — three steps toward the carbon divide:

    1. Your two figures tell the same story two ways. Which is clearer for a reader who has never seen a regression table — the coefficient plot or the prediction plot?
    2. The prediction plot looks clean and confident. But in Exercise 1 that same slope reversed once the petro economies were dropped. Could you honestly publish this figure on its own?
    3. The countries emitting the least are often the most exposed to climate damage. What does that mean for reading a "positive slope" as good or bad news?

(a) The prediction plot needs no statistics vocabulary — it reads like a line on a graph; the coefficient plot is compact but assumes the reader knows what a reference category is. (b) No. The figure is perfectly correct and still misleading on its own: it is drawn from a slope Exercise 1 showed is not robust, because the high-state-ownership end is populated by petro economies. A plot cannot show you what it excluded — which is why the robustness check belongs in the paper, not just in your script. (c) A positive slope is neither good nor bad by itself — it is the injustice that the low emitters bear the highest costs that turns the description into a research (and policy) agenda.

Bonus (for the fast): swap the coefficient plot's predictor from income_level to region. Which reference category does R pick, and how would you change it to "Sub-Saharan Africa"?

# R uses the first factor level (alphabetical for a character vector) as the
# reference — so set it explicitly with fct_relevel(), as we did for income.
Dat <- Dat %>% mutate(region = fct_relevel(region, "Sub-Saharan Africa"))

lm_robust(co2 ~ region, data = Dat) %>%
  tidy() %>%
  mutate(term = str_remove(term, "region")) %>%
  filter(term != "(Intercept)") %>%
  ggplot(aes(y = estimate, x = reorder(term, estimate))) +
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  coord_flip() +
  labs(x = "", y = "Difference in CO2 vs. Sub-Saharan Africa") +
  theme_minimal()