(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.)
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()
Read your plot: compared with low-income countries, upper-middle-income countries emit about tonnes more COâ‚‚ per person (two digits).
Which statement matches your coefficient plot?
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()
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).
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.
# 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.
Discuss with your neighbour — three steps toward the carbon divide:
(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()