1. Do the Setup, if you have not yet done so.

  2. Start a new R script and copy in this starter code. It is the lecture's World Bank pipeline, joined to V-Dem — the same join you practised in Lecture 2:

pacman::p_load(
  tidyverse,    # Data manipulation & figures
  wbstats,      # World Bank data via API
  vdemdata,     # Varieties of Democracy data
  estimatr,     # OLS with robust standard errors
  modelsummary  # Nicely formatted regression tables
)

# World Bank: CO2 and GDP per person, most recent year we observe both
Dat_co2 <- wb_data(c(co2 = "EN.GHG.CO2.PC.CE.AR5",
                     gdp = "NY.GDP.PCAP.PP.KD"),
                   start_date = 2000, end_date = 2026) %>%
  rename(country_text_id = iso3c, year = date) %>% # Rename to V-Dem's key
  select(country_text_id, country, year, co2, gdp) %>%
  drop_na(co2, gdp) %>%
  group_by(country) %>%
  filter(year == max(year)) %>%
  ungroup() %>%
  mutate(gdp_k = gdp / 1000)

# V-Dem: state ownership of the economy (our "how socialist is it?" measure)
Dat_vdem <- vdem %>%
  as_tibble() %>%
  select(country_text_id, year, state_own_raw = v2clstown) %>%
  mutate(state_ownership = -state_own_raw) # Reverse: higher = MORE state ownership

# Region & income group, as in the lecture
Dat_meta <- wb_countries() %>%
  select(country_text_id = iso3c, region, income_level)

# The join: a country is identified by its ISO code AND the year
Dat <- Dat_co2 %>%
  inner_join(Dat_vdem, by = c("country_text_id", "year")) %>% #<<
  inner_join(Dat_meta, by = "country_text_id")

Stuck?

  1. First, check what the join cost you. The World Bank gives us ~190 countries, but Dat has rows. Nothing went wrong — V-Dem simply does not code every microstate. (Our Lecture 4 artefact, Palau, is one of the countries it drops.) An inner_join() is a silent filter: always look at the row count.

  2. Now write it yourself. The World Bank sorts countries into four income groups. Make income_level an ordered factor with "Low income" as the reference, then regress co2 on it with lm_robust(). (This is today's categorical-predictor skill.)

Use fct_relevel() inside mutate() to set the reference level, then lm_robust(co2 ~ income_level, data = Dat).

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

ols_inc <- lm_robust(co2 ~ income_level, data = Dat)

modelsummary(list("CO2 per person" = ols_inc), stars = TRUE,
             coef_rename = c("income_levelLower middle income" = "Lower middle",
                             "income_levelUpper middle income" = "Upper middle",
                             "income_levelHigh income" = "High income"),
             gof_map = c("nobs", "r.squared"))
CO2 per person
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 0.237***
(0.061)
High income 8.275***
(1.024)
Lower middle 0.833***
(0.143)
Upper middle 3.892***
(0.418)
Num.Obs. 165
R2 0.339

  1. Read off your table: on average, high-income countries emit tonnes of CO₂ per person more than low-income countries (two digits; point or comma is fine).

  2. Write one sentence in your R script (as a # comment) interpreting that coefficient — descriptively, as in the lecture. Then compare:

Across the world's countries, high-income countries emit on average about 8.3 tonnes more CO₂ per person than low-income countries — a difference of group means, not (yet) a causal effect.

(Name the direction, the rough size, and the units — and use associational, not causal, language.)

  1. Now the continuous predictor, and today's real question: does state control of the economy go with higher emissions? Make a scatter plot with co2 on the Y-axis and state_ownership on the X-axis. Ask it the four questions: direction, form, spread, outliers.
ggplot(data = Dat, aes(y = co2, x = state_ownership)) +
  geom_text(aes(label = country), size = 3) +
  geom_smooth(method = "lm", formula = y ~ x, color = "#901A1E") +
  labs(y = "Tonnes of CO2 per person",
       x = "State ownership of the economy (V-Dem, higher = more)") +
  theme_minimal()

  1. Before you compute anything: judging only from your scatter plot, the correlation between co2 and state_ownership is …

  2. Fit lm_robust(co2 ~ state_ownership, data = Dat). The slope is positive — tempting headline: state-run economies pollute more. But look at how big its standard error is, then re-fit with base lm() and run the two diagnostic plots. Which statement fits what you see? 

ols_state <- lm_robust(co2 ~ state_ownership, data = Dat)
modelsummary(list("CO2 per person" = ols_state), stars = TRUE,
             gof_map = c("nobs", "r.squared"))
CO2 per person
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 4.869***
(0.722)
state_ownership 1.030
(0.634)
Num.Obs. 165
R2 0.036 
mod <- lm(co2 ~ state_ownership, data = Dat)
plot(mod, which = 5) # Outliers (Cook's D)

plot(mod, which = 1) # Linearity

Bonus (for the fast): those Gulf states nationalised their oil, which is a very different thing from a socialist economy. Flag the petro economies the way the lecture did — fuel rents above 10% of GDP — drop them, and re-fit. Does the headline survive?

# Oil + gas rents, same definition as the lecture's "Petro vs. rest" panel
petro_flag <- wb_data(c(oil = "NY.GDP.PETR.RT.ZS",
                        gas = "NY.GDP.NGAS.RT.ZS"),
                      start_date = 2000, end_date = 2026) %>%
  drop_na(oil, gas) %>%              # Years with both rent series
  mutate(rents = oil + gas) %>%      # Oil + gas rents, in % of GDP
  group_by(iso3c) %>% filter(date == max(date)) %>% ungroup() %>%
  select(country_text_id = iso3c, rents)

Dat <- Dat %>%
  left_join(petro_flag, by = "country_text_id") %>% # keeps countries without rents data
  mutate(petro = case_when(
    rents > 10 ~ "Petro economy",       # fuel rents above 10% of GDP
    TRUE ~ "Other country"))

m_all   <- lm_robust(co2 ~ state_ownership, data = Dat)
m_noOil <- lm_robust(co2 ~ state_ownership, data = Dat %>% filter(petro == "Other country"))

modelsummary(list("All countries" = m_all, "Without petro economies" = m_noOil),
             stars = TRUE, gof_map = c("nobs", "r.squared"))
All countries Without petro economies
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 4.869*** 2.773***
(0.722) (0.348)
state_ownership 1.030 -0.654*
(0.634) (0.309)
Num.Obs. 165 141
R2 0.036 0.032

Look at what happens. With every country in, the slope is positive but its confidence interval comfortably includes zero — there was never a solid finding to report. Drop the petro economies and the slope flips sign and becomes significant: among the other 141 countries, more state ownership goes with less CO₂ per person.

So the "state-run economies pollute" headline was never about socialism. It was a couple of dozen petro economies that happen to keep their oil and gas in state hands — exactly the countries the diagnostic plots flagged in step 9.

Two footnotes worth knowing. (1) Plain lm() would have called the first slope significant (\(p\) = 0.015); lm_robust() does not (\(p\) = 0.11), because those same petro outliers make the errors wildly heteroscedastic. That is why we use robust standard errors. (2) Dropping cases you dislike is not a licence to hunt for the answer you want — it is a robustness check, and the honest write-up reports both models. Remember that for your term paper.

Special (for the very fast): you have only ever controlled for nothing so far. Wealth is the obvious rival explanation — so put gdp_k in the model alongside state_ownership and see what happens to the slope. (This is Lecture 8's skill; today just run it and look.)

m_biv <- lm_robust(co2 ~ state_ownership, data = Dat)
m_mlt <- lm_robust(co2 ~ state_ownership + gdp_k, data = Dat) #<<

modelsummary(list("Bivariate" = m_biv, "Controlling for wealth" = m_mlt),
             stars = TRUE, gof_map = c("nobs", "r.squared"))
Bivariate Controlling for wealth
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 4.869*** 1.457**
(0.722) (0.449)
state_ownership 1.030 1.627**
(0.634) (0.533)
gdp_k 0.143***
(0.031)
Num.Obs. 165 165
R2 0.036 0.475

Surprise: controlling for wealth does not shrink the slope, it makes it bigger (about 1.03 → 1.63) and now clearly significant, while \(R^2\) jumps from about 0.04 to 0.48. That happens because state ownership is mildly negatively correlated with wealth (r ≈ −0.16): state-run economies tend to be poorer, and poverty holds emissions down. Wealth was masking the association, not inflating it — you will meet this as a suppressor in Lecture 8.

Now put the two bonuses together. Fit all four models — with and without the petro economies, with and without the wealth control — and read the pattern.

modelsummary(
  list("All"            = lm_robust(co2 ~ state_ownership, data = Dat),
       "All + wealth"   = lm_robust(co2 ~ state_ownership + gdp_k, data = Dat),
       "No petro"       = lm_robust(co2 ~ state_ownership, data = Dat %>% filter(petro == "Other country")),
       "No petro + wealth" = lm_robust(co2 ~ state_ownership + gdp_k, data = Dat %>% filter(petro == "Other country"))),
  stars = TRUE, gof_map = c("nobs", "r.squared"))
All All + wealth No petro No petro + wealth
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 4.869*** 1.457** 2.773*** 1.082***
(0.722) (0.449) (0.348) (0.265)
state_ownership 1.030 1.627** -0.654* 0.142
(0.634) (0.533) (0.309) (0.220)
gdp_k 0.143*** 0.093***
(0.031) (0.012)
Num.Obs. 165 165 141 141
R2 0.036 0.475 0.032 0.496

Read the state_ownership row across the four columns: +1.03 (n.s.), +1.63 (p < 0.01), −0.65 (p < 0.05), +0.14 (n.s.). Every apparently solid result disappears the moment you take the petro economies out. The honest conclusion is not "socialism pollutes", nor "socialism is green" — it is that these data cannot tell us, because the countries with the most state ownership are largely the countries sitting on oil.

That is a perfectly publishable finding, and a much better answer than picking whichever of the four columns you liked best.