OLS Wisdoms

Assumptions, categorical predictors & how to show your results

Merlin Schaeffer · Department of Sociology

2026-10-26

By the end of today you can …

1. check the OLS assumptions that matter: spot outliers with Cook’s D, judge the linearity assumption, and decide what to do about a suspicious case;

2. use a categorical predictor: dummy-code it, choose a reference category, and read the coefficients as differences between groups;

3. visualize regression results: a coefficient plot for categorical predictors, and predicted values for a continuous one.

Today is craft, not new theory: they can already fit a line, so now they learn to distrust one. Frame it as three habits — interrogate the line, handle categories, show the result. The running example is the carbon divide, which is worth a sentence: this is their generation’s justice question, and we are doing it with the same OLS they already know.

Does the line fit?

Part 1 of 3

Regression is easy to run and easy to fool.

Research question of the day

The climate crisis is the defining North–South injustice of your generation: the countries that emit most are rarely the ones that suffer most.

The gap is staggering: the average Qatari emits hundreds of times more CO₂ than the average person in Chad or the DR Congo (Chancel, 2022).

Research question of the day: across the world’s countries, how does CO₂ emitted per person relate to a country’s wealth?

Earth at night, a satellite composite. Where the lights are, the carbon is. Image: NASA / NOAA, public domain (Wikimedia Commons).

Name the injustice before the statistics: the countries that emit most are rarely the ones that suffer most. The night map says it without numbers: the lights trace wealth and energy use, and most of Africa is dark. Then set today’s much blunter question — do richer countries simply emit more? Be explicit that this is a descriptive question; nothing today is causal, and they should notice we are not claiming it is.

Preparation

• Packages for today’s session
• Get the data
• Build one tibble

pacman::p_load( # Load (and install if needed) several R packages
  tidyverse,    # Data manipulation and visualization
  wbstats,      # Download data from the World Bank API
  estimatr,     # OLS with robust standard errors
  modelsummary  # Nicely formatted regression tables
)

# Both indicators in ONE call: name them, and wbstats names the columns for us
wb_raw <- wb_data(c(co2 = "EN.GHG.CO2.PC.CE.AR5", # CO2, tonnes per person
                    gdp = "NY.GDP.PCAP.PP.KD"),   # GDP p.c., PPP (const. 2021 $)
                  start_date = 2000, end_date = 2026)

# Region & income group are part of a different API, so a second call
wb_meta <- wb_countries()

(Dat <- wb_raw %>%
   drop_na(co2, gdp) %>%     # Keep country-years where we observe both
   group_by(country) %>%
   filter(date == max(date)) %>% # ...and of those, the most recent one
   ungroup() %>%
   inner_join(wb_meta %>% select(iso3c, region, income_level), by = "iso3c") %>%
   select(iso3c, country, date, co2, gdp, region, income_level) %>%
   mutate(
     gdp_k  = gdp / 1000, # GDP in THOUSANDS of $, for readable slopes
     region = fct_relevel(region, "Sub-Saharan Africa") # reference group
   ))
# # A tibble: 191 × 8
#    iso3c country               date    co2    gdp region                                         income_level gdp_k
#    <chr> <chr>                <dbl>  <dbl>  <dbl> <fct>                                          <chr>        <dbl>
#  1 ABW   Aruba                 2024  5.15  45428. Latin America & Caribbean                      High income  45.4 
#  2 AFG   Afghanistan           2024  0.283  1964. Middle East, North Africa, Afghanistan & Paki… Low income    1.96
#  3 AGO   Angola                2024  0.742  8753. Sub-Saharan Africa                             Lower middl…  8.75
#  4 ALB   Albania               2024  1.79  21641. Europe & Central Asia                          Upper middl… 21.6 
#  5 ARE   United Arab Emirates  2024 18.3   69702. Middle East, North Africa, Afghanistan & Paki… High income  69.7 
#  6 ARG   Argentina             2024  4.01  26772. Latin America & Caribbean                      Upper middl… 26.8 
#  7 ARM   Armenia               2024  2.47  20079. Europe & Central Asia                          Upper middl… 20.1 
#  8 ATG   Antigua and Barbuda   2024  3.80  29299. Latin America & Caribbean                      High income  29.3 
#  9 AUS   Australia             2024 14.1   60310. East Asia & Pacific                            High income  60.3 
# 10 AUT   Austria               2024  6.32  63788. Europe & Central Asia                          High income  63.8 
# # ℹ 181 more rows

CO₂ emissions across the world

The extremes are oil-and-gas economies (Qatar, Kuwait, Brunei); at the other end sit some of the world’s poorest states.

Let the bars do the talking — ask them to find Denmark before you point at it. The shape of the distribution is the lesson: petro-states and rich anglophone countries at one end, the world’s poorest at the other. Plant Qatar in their minds now; it comes back as our outlier in ten minutes.

Does wealth drive emissions?

Ask every scatter plot four questions:

1. Direction of the relationship?
2. Form; is it a straight line?
3. Spread around the trend?
4. Any outliers?

Drill the four questions — direction, form, spread, outliers — and make them answer all four out loud before any model is fitted. This is a habit, not a slide: every scatter, every time. They should already be uneasy about the top-right corner.

A first regression

# OLS with robust standard errors
ols_1 <- lm_robust(
  co2 ~ gdp_k,
  data = Dat
)

modelsummary(
  list("CO2 per person" = ols_1),
  stars = TRUE,
  gof_map = c("nobs", "r.squared"),
  output = "kableExtra"
)

	CO2 per person
(Intercept)	1.445+
	(0.742)
gdp_k	0.112***
	(0.026)
Num.Obs.	191
R2	0.151
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Discuss: put the gdp_k coefficient into one plain sentence. How much of what goes with how much of what — and in which units?

Do not read the slope out for them — make them try it. Reading a coefficient aloud, with units, is the habit worth more than the number itself, and it is the thing that fails in exams. Listen for the three classic slips: dropping the units (“emissions go up by 0.12”), forgetting we rescaled to thousands of dollars, and sliding into causal language (“makes”, “causes”, “leads to”). Collect two or three attempts, then turn the page.

Reading the slope

# OLS with robust standard errors
ols_1 <- lm_robust(
  co2 ~ gdp_k,
  data = Dat
)

modelsummary(
  list("CO2 per person" = ols_1),
  stars = TRUE,
  gof_map = c("nobs", "r.squared"),
  output = "kableExtra"
)

	CO2 per person
(Intercept)	1.445+
	(0.742)
gdp_k	0.112***
	(0.026)
Num.Obs.	191
R2	0.151
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Each extra $1,000 of GDP per person comes with about 0.11 more tonnes of CO₂ per person.

Put their best attempt next to this one. Three things to make explicit: the $1,000 (because we divided by 1,000 ourselves — the slope would be a near-invisible 0.0001 in raw dollars), the units on both sides, and the words “comes with” rather than “causes”. Then point at the low \(R^2\) and ask the question the whole first part of today answers: should we trust this line? The honest answer is that we do not know yet, because we have not looked at how it was produced.

OLS makes assumptions

For OLS to describe the data honestly:

1. No influential outliers
2. Linearity (a straight line fits)
3. Homoscedasticity; don’t worry: lm_robust() handles it.
4. Independent observations.

\(\rightarrow\) The first two we check with plots.

Beware: a single strange country, or a curved relationship, can bend the OLS line away from the pattern in the bulk of the data, and quietly change your conclusion.

A model is a reduced representation of reality. It should capture the general pattern; not be dictated by one or two singular cases.

Four assumptions, but be honest that only two need their attention today — lm_robust() already handles heteroscedasticity, and independence is a design question. The sentence to leave them with is the last one: a model should describe the bulk of the data, not be dictated by one or two strange cases.

Assumption 1: outliers

• Standardized residuals: how many far (in standard deviations) a country sits from the OLS line.

  • Qatar (146): CO₂ 47.3 t: extremely rich and extremely polluting.
  • Palau (140): CO₂ 82.8 t: only middling GDP, yet off the chart polluting.

• Leverage: influence on the OLS line; \(x_i\) far from \(\bar{x}\) (Qatar’s huge GDP).

• Cook’s D: how much all predictions shift if case \(i\) were removed \((\hat{\beta} - \hat{\beta_{-i}})\).

  • Beyond the dashed grey lines \(\rightarrow\) influential.

# Re-estimate with base lm(), then the best outlier plot:
lm(co2 ~ gdp_k, data = Dat) %>%
  plot(which = 5)

Three ideas, and they are genuinely different — keep them apart. Residual = far from the line. Leverage = far out on the x-axis, so the line has to pay attention. Cook’s D combines both: how much everything moves if this case vanished. Qatar and Palau are the two names past the dashed lines; ask which is which before you read them out.

Assumption 2: linearity

• Fitted values \(\hat{Y}\) on the X-axis; residuals on the Y-axis.

• Red line: the smoothed relationship between the two.

• Ideal: the red line is flat \(\rightarrow\) OLS captures the relationship well.

• Our case: the red line bends slightly; the linearity assumption is hardly strained.

# Same model, the best linearity plot: #<<
lm(co2 ~ gdp_k, data = Dat) %>%
  plot(which = 1)

Flat red line means the straight line caught everything; a bend means it did not. Ours bends — and note why: the same rich, high-emitting countries cause both problems, so the outlier plot and this plot are pointing at the same cases. Resist the temptation to fix the curve with a log here; that is Lecture 12’s job. Today we only diagnose.

Not every outlier is an error

Discuss: the plots flag two countries. Palau, which is a tiny Pacific island whose per-person figure is distorted by a tanker-refuelling port and ~18,000 residents. Qatar, which is a genuine petro-state. Should we delete them?

Palau — Koror, the main island. Photo: Luka Peternel, CC BY-SA 4.0.

Qatar — the Ras Laffan LNG terminal. Photo: Matthew Smith, CC BY 2.0.

The most important slide of the lecture, and the least statistical. Put faces to the two names, then make them vote on “delete or keep?” before you turn the page — the verdict is the next slide.

So: keep or drop?

Tough call! Palau is a real case, but a measurement artefact for our question (national wealth → emissions): one can argue to drop it. Qatar is a real, important case; we keep it, but stay aware it has high leverage.

Never delete a case just because it is an inconvenient outlier by DFbeta standards.

Their vote first, then this. The distinction: Palau’s number is an artefact of how it was measured — a refuelling port divided by 18,000 residents — so it does not answer our question. Qatar’s number is true and substantively interesting. Delete the artefact, keep the inconvenient truth. Say it flatly: dropping cases because they spoil your result is misconduct, and the defence is that you state what you dropped and why.

Removing the artefact

Dat <- Dat %>%
  filter(country != "Palau")

With Palau gone, Qatar becomes a clear outlier. Sometimes, removing an outlier starts a cascade.

• Outliers again
• Linearity again
• The scatter
• Petro vs. rest
• Log–log

Two lines is an interaction — wealth buying a different amount of CO₂ depending on where the wealth comes from. We only put that in a model in Lecture 11; today, just look. The story: wealth drives emissions above all when the wealth is fossil fuel.

On log axes the bend disappears and the cloud tightens around a straight line — the same data, a fairer description. Re-scaling a curved relationship this way is Lecture 12’s job; today we only notice that it would help.

The line, re-fitted

ols_1 <- lm_robust(co2 ~ gdp_k, data = Dat)

modelsummary(
  list("CO2 per person" = ols_1),
  stars = TRUE,
  gof_map = c("nobs", "r.squared"),
  output = "kableExtra"
)

Same specification, one artefact removed. The slope is now more stable and \(R^2\) has roughly doubled: dropping one distorting case sharpened the whole picture.

One row out of 190-odd, and \(R^2\) roughly doubles — that is the number to dwell on. It should feel slightly alarming, and it should: a single case had that much say over the fit. Good moment to ask what would have happened if we had never looked at the diagnostic plots.

	CO2 per person
(Intercept)	0.888+
	(0.486)
gdp_k	0.117***
	(0.026)
Num.Obs.	190
R2	0.343
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Comparing groups

Part 2 of 3

Wealth is continuous. But much of what we compare is categorical.

CO₂ by world region

Read the gradient across regions out loud, then plant the diamonds: those are the region means, and the means are what a regression on region will predict. Point out that they sit above the median line in most regions, because a few heavy emitters pull the mean up; OLS predicts the mean, not the median. Do not explain why yet — the next two slides make them work it out. Once they have, dummy coding stops being a trick and becomes obvious.

A dummy is a regression line

Take just two regions and code them 0/1:

1. What is \(\hat{Y}\) for a Sub-Saharan African country?
2. What is \(\hat{Y}\) for a European one?
3. How does the difference between them relate to \(\hat{\beta}\)?

This is the slide that makes dummies click, so do not rush it, and do not answer it yourself. Nothing new is happening here: it is the same OLS line as this morning, drawn over an X that only ever takes the values 0 and 1. Work the blue box in order and make them commit to a number for each of the first two questions before you turn the page. The diamonds are the same mean markers they just saw on the boxplot — point at them, say nothing more, and wait.

\(\hat{\beta}\) is a difference in means

1. \(\hat{Y}\) at X = 0 is the intercept: 0.92 t — Sub-Saharan Africa’s mean.

2. \(\hat{Y}\) at X = 1 is intercept + slope: 5.48 t — Europe’s mean.

3. So \(\hat{\beta}\) = 4.56 is simply the difference between two group means.

Now confirm it, and make the geometry explicit: the line is pinned through both diamonds. With only two possible X values there is nowhere else for a straight line to go — it must hit the mean of each column, so the slope cannot be anything other than the gap between those means. That is the whole idea of a dummy coefficient, and it is why “average difference from the reference” is a description of a slope rather than a new kind of quantity. Then generalise out loud: with more than two categories R just draws one of these lines per dummy, all sharing the same reference. The next slide is only the bookkeeping for doing that.

Dummy coding: from labels to 0/1

A regression needs numbers, not words. R turns each category into a dummy: 1 if the case is in that category, 0 otherwise — leaving one category out as the reference.

\[x = \begin{cases} 1, & \text{if in the category} \\ 0 & \text{otherwise} \end{cases}\]

Country	N. America	Europe	…
USA	1	0	0
Canada	1	0	0
Denmark	0	1	0
Germany	0	1	0
…	0	0	1
Reference (Sub-Saharan Africa, see slide 5)	0	0	0

Each coefficient = that region’s average difference in CO₂ from the reference.

Walk one row of the table so they see it: the USA gets a 1 in North America and 0 everywhere else. Then the bit that always confuses people — one category has to be left out, and it becomes the reference. Every other coefficient is a comparison to that one. Ask what would happen if we included all of them; you do not need the algebra, just the intuition that the last column is redundant.

Categorical predictors in R

• The code
• The table
• Interpretation

# R dummy-codes a factor (or character) automatically,
# using its FIRST level as the reference category.
ols_2 <- lm_robust(co2 ~ region, data = Dat)

modelsummary(
  list("CO2 per person" = ols_2),
  stars = TRUE,
  # Shorten the region labels for the table
  coef_rename = c("regionEast Asia & Pacific" = "East Asia & Pacific",
                  "regionEurope & Central Asia" = "Europe & Central Asia",
                  "regionLatin America & Caribbean" = "Latin America",
                  "regionMiddle East, North Africa, Afghanistan & Pakistan" = "Middle East & N. Africa",
                  "regionNorth America" = "North America",
                  "regionSouth Asia" = "South Asia"),
  gof_map = c("nobs", "r.squared"),
  output = "kableExtra"
)

	CO2 per person
(Intercept)	0.922***
	(0.224)
East Asia & Pacific	3.320***
	(0.926)
Europe & Central Asia	4.557***
	(0.463)
Latin America	2.109***
	(0.573)
Middle East & N. Africa	8.617***
	(2.571)
North America	9.667**
	(3.236)
South Asia	0.756
	(0.517)
Num.Obs.	190
R2	0.241
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

• Reference (Sub-Saharan Africa): average CO₂ is the intercept, 0.92 t per person; the value when all region dummies are 0.

• North America vs. Sub-Saharan Africa: 9.67 t per person more, on average → an average of (0.92 + 9.67 =) 10.59 t.

• Every coefficient is a difference of group means from the reference. No line, no linearity assumption to check.

Read the intercept first — it is Sub-Saharan Africa’s average, the case where every dummy is zero. Then North America as a difference from that, and add the two to get its actual average. Flag the last bullet explicitly: with a categorical predictor there is no line, so the linearity assumption from Part 1 simply does not apply. Students routinely try to check it anyway.

Break

Your turn: exercise 1

You practice today’s diagnostics and a categorical predictor. This time with income groups instead of regions.

Open exercise 1 in a new tab ↗

Break

Showing your results

Part 3 of 3

Figures are often much better for your reader. Two go-to plots: coefficient plots, and predictions.

(1) Coefficient plots

• Preparation
• Plotting
• Plot

(plotdata <- lm_robust(co2 ~ region, data = Dat) %>%
   tidy() %>% # Turn the results into a tibble #<<
   mutate(term = term %>% # Tidy up the region labels
            str_remove("region") %>%
            str_replace(", North Africa, Afghanistan & Pakistan", " & N. Africa")) %>%
   filter(term != "(Intercept)")) # Drop the reference
#                        term estimate std.error statistic  p.value conf.low conf.high  df outcome
# 1       East Asia & Pacific    3.320     0.926      3.59 4.31e-04    1.493      5.15 183     co2
# 2     Europe & Central Asia    4.557     0.463      9.85 1.27e-18    3.644      5.47 183     co2
# 3 Latin America & Caribbean    2.109     0.573      3.68 3.07e-04    0.978      3.24 183     co2
# 4   Middle East & N. Africa    8.617     2.571      3.35 9.75e-04    3.545     13.69 183     co2
# 5             North America    9.667     3.236      2.99 3.20e-03    3.281     16.05 183     co2
# 6                South Asia    0.756     0.517      1.46 1.46e-01   -0.265      1.78 183     co2

ggplot(data = plotdata,
       aes(y = estimate,
           x = reorder(term, estimate))) + # Order by size
  # Reference line at "no difference"
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  # A point with its 95% confidence interval
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + #<<
  coord_flip() + # Flip the axes so labels are readable #<<
  labs(x = "World region",
       y = "Average difference in CO2 (t per person)\nvs. Sub-Saharan Africa") +
  theme_minimal(base_size = 15)

The pitch: nobody reads a table of seven coefficients, but everyone reads this picture in two seconds. Point out the three deliberate choices — the dashed line at zero for “no difference”, ordering by size, and coord_flip() so the labels are legible. They rebuild this in exercise 2, so linger on the code tab.

(2) Model predictions

• The model
• Step 1: synthetic data
• Step 2: predict
• Step 3: the code
• Step 4: the plot

# Back to the continuous predictor
(ols_1 <- lm_robust(co2 ~ gdp_k, data = Dat))

#             Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper  DF
# (Intercept)    0.888     0.4859    1.83 6.91e-02  -0.0702    1.847 188
# gdp_k          0.117     0.0255    4.60 7.90e-06   0.0669    0.168 188

# We want the line for a RANGE of wealth levels, so we invent countries to
# predict FOR. They are fictional on purpose: no country has to exist at
# exactly $37k - we only need X values to feed the model.
(fict_dat <- tibble(
  # One fictional country per $1,000 step, from $1k up to $120k.
  # seq() builds 1, 2, 3, ... 120 - stay inside the OBSERVED range,
  # because a model says nothing about wealth levels it never saw.
  gdp_k = seq(from = 1, to = 120, by = 1)
))
# # A tibble: 120 × 1
#    gdp_k
#    <dbl>
#  1     1
#  2     2
#  3     3
#  4     4
#  5     5
#  6     6
#  7     7
#  8     8
#  9     9
# 10    10
# # ℹ 110 more rows

(fict_dat <- predict(
  object   = ols_1,      # The fitted model: what we learned from real countries
  newdata  = fict_dat,   # The fictional countries: where we want the line
  # "confidence" = uncertainty about the AVERAGE prediction (a narrow ribbon),
  # NOT where individual countries fall. level = 0.95 -> the usual 95%.
  interval = "confidence", level = 0.95)$fit %>% # $fit = the fit/lwr/upr matrix
   as_tibble() %>%                # predict() hands back a matrix; make it a tibble
   bind_cols(fict_dat, .))        # Glue it beside the gdp_k it belongs to
# # A tibble: 120 × 4
#    gdp_k   fit    lwr   upr
#    <dbl> <dbl>  <dbl> <dbl>
#  1     1  1.01 0.0928  1.92
#  2     2  1.12 0.255   1.99
#  3     3  1.24 0.417   2.06
#  4     4  1.36 0.578   2.14
#  5     5  1.47 0.739   2.21
#  6     6  1.59 0.898   2.29
#  7     7  1.71 1.06    2.36
#  8     8  1.83 1.21    2.44
#  9     9  1.94 1.37    2.52
# 10    10  2.06 1.52    2.60
# # ℹ 110 more rows

ggplot(data = fict_dat,        # The FICTIONAL countries, not Dat #<<
       aes(y = fit,            # fit = the predicted CO2 #<<
           x = gdp_k)) +
  # Two benchmarks to talk about: a poor country and a rich one
  geom_vline(xintercept = c(5, 60), color = "#901A1E", lty = "dashed") +
  # The 95% interval, drawn FIRST so the line sits on top of it
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) + #<<
  geom_line(linewidth = 1) + # The predicted line itself #<<
  labs(x = "GDP per person (thousands of $, PPP)",
       y = "Predicted tonnes of CO2 per person") +
  theme_minimal(base_size = 14)

Four tabs now, but still three ideas, and step 1 is the one they forget: invent fictional countries spanning the GDP range, then predict for them. Walk the comments in step 2 — newdata, and the fact that predict() returns a matrix that has to be glued back beside its gdp_k, are the two places their own code will break. Stress that the ribbon is uncertainty about the average prediction, not a range where countries fall — that misreading is extremely common. The two dashed benchmarks make the divide concrete: a poor country and a rich one, and the gap between their predictions.

The carbon divide, in two pictures

Richer countries and richer regions emit far more per person. Yet the people who emit the least (across Sub-Saharan Africa and South Asia) are the ones most vulnerable to the global warming those emissions cause (Chancel, 2022).

Both plots together, then the point that is not statistical at all: the regions at the bottom of the coefficient plot are the ones that will suffer most from the warming. Our slope describes the divide; it does not explain or excuse it. Say plainly that this is where description ends and the research agenda begins.

Your turn: exercise 2

You practice visualizing regression: a coefficient plot and a prediction plot of your own.

Open exercise 2 in a new tab ↗

One last wisdom: whose carbon is it?

Everything so far counted CO2 where it is produced. But Denmark burns little at home and imports carbon-heavy goods, like steel, electronics, and gasoline.

Consumption-based accounting re-assigns them to whoever consumes the goods: “Green” Denmark’s footprint jumps from 4.8 to 8.3 t; factory economies (China, India) produce what we outsourced.

The ship is Danish on purpose: Maersk carries a large share of the world’s containers, and those containers are how carbon gets imported. This lands well with a Danish audience: everything so far counted carbon where it is burned, and by that measure Denmark looks rather good. Re-assign imported steel, electronics and meat to whoever consumes them and the Danish figure jumps sharply, while China and India fall. Read the dumbbells left to right. The wisdom is about the outcome variable, not about climate: how you choose to measure \(Y\) can reverse who looks responsible.

Whose carbon is in these containers? Photo: Quintin Soloviev, CC BY 4.0.

Run the fit again — with the consumption footprint

• Fit both models
• The comparison
• The picture

# Join each country's consumption footprint, then run the SAME regression
Dat_cons <- Dat %>%
  inner_join(
    readRDS("data/owid_consumption.rds") %>%
      select(iso3c = iso_code, consumption),
    by = "iso3c"
  )

modelsummary(
  list("Produced (WB)"  = lm_robust(co2 ~ gdp_k,         data = Dat_cons),
       "Consumed (GCP)" = lm_robust(consumption ~ gdp_k, data = Dat_cons)),
  stars = TRUE,
  coef_rename = c("gdp_k" = "GDP per person ($1,000s)"),
  gof_map = c("nobs", "r.squared"),
  output = "kableExtra"
)

	Produced (WB)	Consumed (GCP)
(Intercept)	0.846	0.590
	(0.809)	(0.588)
GDP per person ($1,000s)	0.136***	0.175***
	(0.035)	(0.025)
Num.Obs.	118	118
R2	0.358	0.599
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Same predictor, same countries, but a different outcome. Wealth explains far more of the consumption footprint (\(R^2 \approx 0.6\)) than of production (\(R^2 \approx 0.36\)), and the slope is steeper. What a country consumes tracks its wealth more tightly than what it happens to emit at home; and Qatar is no longer the lone outlier.

Same predictor, same countries, one different outcome column — and the fit improves markedly. Wealth tracks what a country consumes far more tightly than what it happens to burn at home. Note the bonus: Qatar stops being the lone outlier, because the production measure was partly an accounting artefact of where the gas is extracted. A nice closing echo of Part 1.

Today’s general lessons

1. Outliers can bend the OLS line. Cook’s D (plot(which = 5)) finds the influential ones.

2. OLS assumes a linear relationship for continuous predictors: check it (plot(which = 1)); it does not apply to categorical predictors.

3. Not every outlier should be deleted: separate artefacts (drop) from real, important cases (keep, but note their leverage).

4. Categorical predictors are dummy-coded: each coefficient is a group’s average difference from the reference category. R does this automatically.

5. Coefficient plots show categorical results; prediction plots show what a continuous model implies across a range of values.

Today’s important functions

• Post-lm() diagnostics:
  • plot(model, which = 5) — identify influential outliers (Cook’s D).
  • plot(model, which = 1) — check the linearity assumption.
• Building the figures:
  • tidy() — turn a model into a tibble of coefficients + confidence intervals.
  • geom_pointrange() + coord_flip() — the coefficient plot.
  • predict(model, newdata = ..., interval = "confidence") — \(\hat{Y}\) for synthetic data.
• fct_relevel() — set which category is the reference.

Check yourself: today’s goals

Look back at the goals from the start. Can you tick all three?

• Run plot(model, which = 5) and which = 1 and say, in words, what each one is telling you.
• Decide whether a flagged outlier is an artefact to drop or a real case to keep and defend your call.
• Read a coefficient on a dummy variable as a difference from the reference, and build a coefficient plot and a prediction plot for your reader.

Anything shaky? That is what this week’s Absalon quiz and the Friday exercise class are for.

References

Chancel, L. (2022). “Global Carbon Inequality over 1990-2019”. In: Nature Sustainability, pp. 931-938.

Appendix: binary outcomes

What if the outcome is yes/no? A short preview of the Linear Probability Model.

LPM vs. generalised linear models

Linear Probability Model (LPM): just OLS with a 0/1 outcome — e.g. is this a high-emitting country? (above 2 t per person).

• It predicts \(\hat{\text{P}}(y_i = 1 \mid x_i)\).
• Simple, and the coefficients read as changes in probability.
• But: it can predict probabilities below 0 or above 1.

The red line pierces 1 at high GDP — a probability above 100%. That is the LPM’s controversial side.

Logistic regression — a GLM for binary outcomes

\[\begin{aligned} \text{logit}(y_i) &= \alpha + \beta x, \\ y_i &= \text{logit}^{-1}(\alpha + \beta x) \\ &= \frac{1}{1 + e^{-(\alpha + \beta x)}}. \end{aligned}\]

A link function squeezes the linear model into the \([0, 1]\) range. The sigmoid-shaped logistic function does exactly that.

GLMs bring their own interpretation traps (Breen, Karlson, and Holm, 2018). In this course we use Linear Probability Models for binary outcomes — and 0/1 coding for categorical ones.