gss <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels()4 Prediction and error
Most of the rest of this book is some version of the same thing. We come up with a rule for predicting an outcome, we see how wrong the rule is, and we ask whether a more complicated rule is enough better to be worth it. In this chapter we start with a very simple model and build up from there.
We use haven::zap_labels() to get rid of the labels that come from the Stata import that is the origin of the data in the gssr package. These can be useful sometimes, but they can also get in the way of certain functions.
4.1 A model is a description of a process
We can write a model for something as simple as the observed values of a single variable, with no predictors at all. The GSS asks respondents to place themselves on a seven-point scale of political views, from 1 (extremely liberal) to 7 (extremely conservative).
pol <- gss |>
select(polviews, sex, wordsum, degree) |>
filter(!is.na(polviews))
nrow(pol)[1] 3160If we want to write a model for each individual observation, we could write something like this:
\[\text{polviews}_i = \text{mean} + \text{weirdness}_i\]
This means that everyone starts at the same place and then departs from it in some way due to their “weirdness” (i.e., difference from what is typical).
Note
The mean doesn’t have a subscript because everyone shares the same one; “weirdness” does because everyone is weird in their own way.
We usually write such a model like this:
\[y_i = \mu + \epsilon_i\]
We think of this model as the data-generating process. Of course it’s very simple, but that’s why it’s a model.
ImportantError and residual are not the same thing
Error (\(\epsilon_i\)) is the part of a person that our model doesn’t capture. It’s a real thing in the world, and we never get to see it.
Residuals are our estimates of the errors, based on a particular model: \(e_i = y_i - \hat{y}_i\). If we fit a different model, the residuals change. The errors don’t.
People use the two words interchangeably all the time. In this book, “error” means the real, unobserved thing and “residual” means the thing we compute.1
1 We often add a second line to the model, \(\epsilon_i \sim \mathcal{N}(0, \sigma)\). This means that the “weirdness” (\(\epsilon_i\)) is normally distributed, meaning that less weirdness is more common than more weirdness! Nothing in this chapter needs that assumption, because sums of squared residuals (and the comparisons we make with them) are just descriptive. We’ll need it starting in Chapter 5, where it gives us the sampling distributions we use for tests and intervals.
4.2 The simplest possible prediction
Suppose we have to guess one number for everyone. The scale runs from 1 to 7, so 4 is the midpoint (“moderate”). Guessing 4 amounts to saying that Americans are, on average, moderate.
pol <- pol |>
mutate(pred_c = 4,
resid_c = polviews - pred_c)
head(pol$resid_c, 10) [1] 0 -1 -3 0 -2 0 0 0 -2 -3Some residuals are positive and some are negative. To see how wrong we are overall, we need to combine them somehow. We can’t just add them up, because the positives and negatives cancel out (a terrible model could have residuals that add up to zero!). So we square them first and then add them up:
\[\text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]
sse_c <- sum(pol$resid_c^2)
sse_c[1] 7431We’ll rely on the sum of squared errors (SSE) to quantify ERROR in this framework, at least for now.2 Squaring gets rid of the cancellation problem. It also means that big misses count a lot more than small ones: being off by 4 is sixteen times as bad as being off by 1, not four times.
2 This is also sometimes called the sum of squared residuals (SSR) or the residual sum of squares (RSS). These are all the same quantities. Given what I said above, “sum of squared residuals” is really the more accurate name. But SSE is what you’ll see most often, so that’s what we’ll use.
4.3 A model that learns from the data
Our first model guessed 4 without looking at the data at all. We can do better by estimating the typical value from the data instead.
pol <- pol |>
mutate(pred_a = mean(polviews),
resid_a = polviews - pred_a)
sse_a <- sum(pol$resid_a^2)
c(sse_c = sse_c, sse_a = sse_a) sse_c sse_a
7431.000 7394.202 The sample mean is 4.108. That’s close to 4, so the improvement is small. But it is an improvement.
The mean is the value that minimizes SSE. No other single guess does better. We can see this by trying a bunch of different guesses:
candidates <- seq(3.5, 4.7, by = 0.01)
sse_for <- map_dbl(candidates, \(g) sum((pol$polviews - g)^2))plt(sse_for ~ candidates, type = "l",
xlab = "Guess", ylab = "Sum of squared errors")
abline(v = mean(pol$polviews), lty = 3, lwd = 2)
ggplot(data.frame(candidates, sse_for), aes(x = candidates, y = sse_for)) +
geom_line() +
geom_vline(xintercept = mean(pol$polviews), linetype = "dotted", linewidth = 1) +
labs(x = "Guess", y = "Sum of squared errors")
The SSE makes a bowl shape, and the bottom of the bowl is at the mean. Fitting a model means choosing parameter values that make the SSE as small as possible. Later on we’ll have more than one parameter to choose, but the idea is exactly the same.
4.4 Comparing two models
Now we have two models and their SSEs. How much better is the second one?
The raw SSE isn’t very useful for this, because its size depends on the units and on how many people we have. A proportion is easier to interpret:
\[\text{PRE} = \frac{\text{SSE}_{\text{C}} - \text{SSE}_{\text{A}}}{\text{SSE}_{\text{C}}}\]
This is the proportional reduction in error. It tells us what proportion of the simpler model’s error the more complex model gets rid of.
The terms here (the compact model C, the augmented model A, and PRE) come from Correll and colleagues (Correll et al. 2025). The idea of comparing models is much older, but I like their terminology.
pre <- (sse_c - sse_a) / sse_c
pre[1] 0.004951929Estimating the mean rather than assuming 4 gets rid of about 0.5% of the error.
This isn’t a lot, but it’s not nothing, either. It just means that model C (Americans are moderate on average) was very nearly right to begin with, so there wasn’t much for model A to improve on. A small PRE can tell you something about the world. It doesn’t have to mean there’s something wrong with the model.
4.4.1 A second look, with countries
The same logic works for other kinds of data, and it works even when we don’t estimate the predictions from the data at all. Here’s the share of each country’s population using the internet in 2021, from the World Bank’s World Development Indicators.
internet <- readRDS(here::here("data", "WDI.rds")) |>
filter(region != "Aggregates") |>
select(country, iso = iso3c, intpct = IT.NET.USER.ZS, income) |>
tidyr::drop_na()
nrow(internet)[1] 180Suppose we just guess that 70% of people are online in every country, a number I picked by eye rather than calculating. Then suppose we make a conditional guess instead: 90% in high-income countries, 70% everywhere else.
internet <- internet |>
mutate(high_income = income == "High income",
guess_flat = 70,
guess_cond = if_else(high_income, 90, 70))
c(flat = sum((internet$intpct - internet$guess_flat)^2),
conditional = sum((internet$intpct - internet$guess_cond)^2)) flat conditional
114666.11 89166.05 sse_flat <- sum((internet$intpct - internet$guess_flat)^2)
sse_cond <- sum((internet$intpct - internet$guess_cond)^2)
(sse_flat - sse_cond) / sse_flat[1] 0.2223853Notice that making separate guesses makes the ERROR (SSR or RSS or SSE) go down. That’s an improvement. Conditioning on income gets rid of about 22% of the error, and neither number was estimated. We made both guesses up!
So a model doesn’t have to be fitted to be a model. It just has to make predictions that we can compare to the data. Fitting is just a way of choosing good predictions. Our made-up guesses did OK because one of them happened to be pretty close:
internet |>
group_by(high_income) |>
summarize(mean_intpct = mean(intpct), n = n())# A tibble: 2 x 3
high_income mean_intpct n
<lgl> <dbl> <int>
1 FALSE 57.3 117
2 TRUE 90.1 63Guessing 90 for high-income countries was almost exactly right. Guessing 70 for the rest was way too high, and that’s where most of the remaining error comes from.
4.5 The same thing with lm()
It’s good to do this by hand once, but from here on we’ll let R fit the models with lm(). We can get R to estimate model A, our one parameter “augmented” estimate, with the formula y ~ 1, which asks it to estimate an intercept and nothing else.
mod_a <- lm(polviews ~ 1, data = pol)
coef(mod_a)(Intercept)
4.107911 deviance(mod_a)[1] 7394.202deviance() gives us the SSE, which matches what we got by hand. The estimated intercept is the sample mean.
Model C, our zero parameter guess, is a little weirder. First let’s make a version of polviews “centered” on 4, and then estimate a model with no constant at all (that is, estimating no parameters), written ~ 0:
mod_c <- lm(I(polviews - 4) ~ 0, data = pol)
deviance(mod_c)[1] 7431This is the same SSE as before. (Because of the way we transformed the outcome variable (deviating it from 4), fitting ~ 1 to it instead would answer “how different is the estimated mean from 4.”)
Now we can get the PRE from the model objects:
(deviance(mod_c) - deviance(mod_a)) / deviance(mod_c)[1] 0.0049519294.6 Predictions that depend on something
So far every model has made the same prediction for everyone. But we can also use the data to make different predictions for different groups. In Chapter 3 we compared conditional probabilities across groups. This is the same idea with a continuous outcome: instead of one mean, we have a mean for each group.
Let’s look at vocabulary scores and education. degree runs from 0 (less than high school) to 4 (graduate degree).
voc <- gss |>
filter(!is.na(wordsum), !is.na(degree))
voc |>
group_by(degree) |>
summarize(mean_wordsum = mean(wordsum), n = n())# A tibble: 5 x 3
degree mean_wordsum n
<dbl> <dbl> <int>
1 0 4.48 187
2 1 5.91 992
3 2 6.08 185
4 3 7.21 477
5 4 7.63 315Scores go up pretty steadily with education. Now let’s compare two models: one that ignores education, and one that gives each degree category its own mean.
uncond <- lm(wordsum ~ 1, data = voc)
cond <- lm(wordsum ~ factor(degree), data = voc)
c(sse_uncond = deviance(uncond), sse_cond = deviance(cond))sse_uncond sse_cond
10711.473 8994.675 pre_cond <- (deviance(uncond) - deviance(cond)) / deviance(uncond)
pre_cond[1] 0.1602765Knowing someone’s degree gets rid of about 16% of the error we made by predicting the overall mean for everyone. That’s a lot more than half a percent!
TipPRE has another name you already know
Look at what the standard model summary reports:
| Unconditional | Conditional on degree | |
|---|---|---|
| (Intercept) | 6.340 | 4.481 |
| factor(degree)1 | 1.432 | |
| factor(degree)2 | 1.600 | |
| factor(degree)3 | 2.724 | |
| factor(degree)4 | 3.144 | |
| Num.Obs. | 2156 | 2156 |
| R2 | 0.000 | 0.160 |
| RMSE | 2.23 | 2.04 |
\(R^2\) for the conditional model is 0.16, which is exactly the PRE we just computed by hand.
When the compact model is the intercept-only model, PRE is \(R^2\). These are two names for the same number, and you’ll see both. PRE is more general, since we can use it to compare any pair of nested models. \(R^2\) is just the name for this particular comparison.
The model summary also reports the RMSE (the square root of the mean squared error). This is roughly the size of a typical residual, in the units of the outcome. It goes down from 2.23 to 2.04 words. So PRE shows the improvement as a proportion and RMSE shows it in words correct.3
3 broom::glance() is also a useful function to look at models. The sigma column (\(\sigma\)) is the RMSE.
4.7 Is the improvement worth it?
This question will keep us busy for the next several chapters.
Adding parameters to a model always reduces the SSE (or at least doesn’t make it bigger). If we take twenty respondents and give each one their own parameter, the SSE goes to zero. That model predicts the data perfectly, but it hasn’t told us anything about the world.
twenty <- voc |> slice_head(n = 20)
perfect <- lm(wordsum ~ factor(seq_len(20)), data = twenty)
round(deviance(perfect), 10)[1] 0So the fact that a more complex model has a lower SSE doesn’t tell us much by itself. So how do we know if it’s “worth it” to adopt the more complex model? We need to know whether the reduction is bigger than we’d expect if the extra parameters were just fitting noise. To answer that, we need to know about sampling distributions, which are the subject of Chapter 5.
4.8 Recap
- a model is a claim about how the data were generated; the simplest is \(y_i = \mu + \epsilon_i\), a common value plus each person’s own departure from it
- error is real and unobserved; residuals are our estimates of it, conditional on a model
- change the model and the residuals change, but the errors don’t
- SSE squares residuals before summing, so positives and negatives can’t cancel
- the mean minimizes SSE, and fitting a model means finding the parameter values that make the SSE as small as possible
- PRE expresses one model’s improvement over another as a proportion of the simpler model’s error
- when the simpler model is intercept-only, PRE is \(R^2\)
- a small PRE isn’t necessarily a problem with the model
- more parameters always reduce SSE, so a reduction by itself proves nothing