Week 8: Models and Their Limits

Linear regression, diagnostics, and knowing what your model can’t tell you

Introduction

A model is a deliberate simplification: it captures some of reality and ignores the rest. The skill is not fitting one — lm() does that in eight characters — but knowing what yours has thrown away.

Today you will fit two models. One is honest about being incomplete. The other is HolmesCo’s, and it is about to be taken apart.

Goals:

Fit a linear regression and say what the slope means in words
Read residual plots for the structure the model missed
Rebuild a bad model and measure how much was being left out
See why a model that fits better can predict worse

How these exercises work

The code boxes start empty, except for comments to guide you through the code that’s expected.

Work one comment at a time. If you get stuck:

Hint 1 restates what you are trying to do.
Hint 2 sketches the shape the code might take.
Hint 3 gives you the code with the key pieces blanked out.
Solution gives you the lot.

Take a real attempt before opening a hint: you will learn by doing, not by reading. The demonstrators are here to help.

Three things worth knowing:

The boxes share one R session, exactly like a script. Anything you create in Exercise 1 is still there in Exercise 2.
Variable names are up to you. The last thing your block prints is what gets checked, so end each block with the answer.
Exercises have extensions underneath them. If you finish early, give them a stab. These are designed to equip you to tackle real coding problems you might face, using the documentation. We can help – but check the manuals first!

Stuck?

Ask on the class board, in Code Q&A. Someone else is likely stuck on the same thing, so asking in public helps them too. Demonstrators check it most weekdays, but they give classmates a chance to answer first.

The board is private to this class, so sign in to GitHub first. If the link shows “404 – page not found”, you’re not signed in, or you haven’t accepted the Classroom 50 invitation yet.

Exercise 1: Does heat help or hurt?

solar holds 50 observations of panel temperature (°C) and electrical output_w (watts).

Draw the relationship: the points, plus a straight-line fit through them, properly labelled. Decide from the picture which way it slopes before you fit anything.

Hint 2

The shape of it:

ggplot(data, aes(x = predictor, y = outcome)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = "...", y = "...", title = "...")

The grey band around the line is a confidence interval for the mean at each temperature — not a range the next panel is likely to fall in. Those are two very different things, and the plot does not distinguish them.

Hint 3

ggplot(solar, aes(x = ______, y = ______)) +
  geom_point() +
  geom_smooth(method = "______") +
  labs(x = "Panel temperature (°C)",
       y = "Output (W)",
       title = "______")

Solution

ggplot(solar, aes(x = temperature, y = output_w)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = "Panel temperature (°C)",
       y = "Output (W)",
       title = "Hotter panels produce less power")

The line slopes downward: output falls as temperature rises. That surprises most people, and it is real — crystalline silicon loses roughly 0.4–0.5% of its efficiency per °C above 25 °C. A hot day brings more sunlight and a less efficient panel.

Look at the scatter around the line before you get attached to it. The points are spread widely at every temperature, so whatever this line is capturing, it is not most of what is going on. The next exercise puts a number on that.

Extension: what the line is hiding

Optional. ?geom_smooth, ?geom_abline, ?predict.lm, ?stat_summary.

Turn the confidence band off, then widen it to 99%. Which version would you publish, and which would a manufacturer’s marketing department prefer?
Fit a loess curve instead of a straight line. Does the curvature it finds look like physics or like noise? Try changing span and watch the “finding” move.
Draw the fitted line without using geom_smooth(). You will need the model’s coefficients and a geom that takes a slope and an intercept. Doing it this way is more work and makes it much harder to plot a line that is not the model you reported.
Add a prediction interval — the range a new panel would fall in — alongside the confidence interval. ?predict.lm and its interval argument. The difference between the two bands is usually startling, and it is the band your reader thinks they are looking at.

Exercise 2: What the slope means

Fit the model properly, read its summary, and finish by printing the slope — the coefficient on temperature.

Then say the slope out loud as a sentence with units in it. If you cannot, you have not finished.

Hint 2

The shape of it:

fit <- lm(outcome ~ predictor, data = df)
summary(fit)

coef(fit)["predictor"]

In the summary, read along the temperature row: estimate, standard error, t value, p-value. Then find Multiple R-squared lower down. A small p-value and a small R-squared together are entirely normal, and they mean “there is definitely something here, and it is not much”.

Solution

solar_model <- lm(output_w ~ temperature, data = solar)
summary(solar_model)

coef(solar_model)["temperature"]

The slope is −0.354: each additional degree costs about a third of a watt. Over the 37 °C range in the data, that is about 13 W out of roughly 250 — a 5% effect, and consistent with the published temperature coefficients for silicon panels.

Three things in that summary are worth more than the slope:

p = 0.0015. The relationship is not a fluke.
R² = 0.19. Temperature explains 19% of the variation. Cloud cover, panel age, angle, dust and inverter losses are all in the other 81%, and none of them is in this model.
The intercept, 257.5 W, is the prediction at 0 °C — a temperature that appears nowhere in these data. Intercepts are routinely reported as though they meant something. Usually they are an artefact of where zero happens to be.

“Significant but small” is the most common honest result in science and the least often reported, because it does not make a headline.

Extension: the same model, interrogated

Optional. ?confint, ?predict.lm, ?scale, ?cor.

Get a confidence interval for the slope. Does it exclude zero — and does it exclude values you would consider practically irrelevant? Those are different questions and only the first has a standard test.
Re-fit with temperature centred on 25 °C — temperature - 25 inside the formula. The slope is unchanged and the intercept now means something. Which version would you report?
Predict the output at 60 °C. R will answer without complaint. Say what is wrong with the answer, and find the widest temperature for which you would defend a prediction.
R² for a simple regression is the square of the correlation. Confirm it with cor(), then explain why that stops being true the moment you add a second predictor.

Quick check

The solar model’s R² is 0.19. What percentage of the variation in output is left unexplained? Print the number.

Exercise 3: HolmesCo’s dominant control

HolmesCo monitored groundwater level and rainfall near a quarry for five years, 60 monthly readings. They fitted one model and wrote:

“Rainfall is the dominant control on groundwater levels (R² = 0.49, p < 0.001).”

Reproduce their model, then look at what it left behind. Residuals — the gap between each observation and the model’s prediction — should look like noise. If they have a pattern, the pattern is something the model does not know about.

Hint 2

The shape of it:

fit <- lm(outcome ~ predictor, data = df)

df$resid <- residuals(fit)

ggplot(df, aes(x = time, y = resid)) +
  geom_point() +
  geom_line(alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "...", y = "Residual (m)")

summary(fit)$r.squared

Joining the points with a faint line makes a cycle far easier to see than the points alone. That is legitimate here because the x axis is time and consecutive points really are adjacent.

Hint 3

gw_model <- lm(______ ~ ______, data = groundwater)

groundwater$resid <- ______(gw_model)

ggplot(groundwater, aes(x = ______, y = resid)) +
  geom_point() +
  geom_line(alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Month", y = "Residual (m)",
       title = "What HolmesCo's model threw away")

summary(gw_model)$______

Solution

gw_model <- lm(gw_level_m ~ rainfall_mm, data = groundwater)

groundwater$resid <- residuals(gw_model)

ggplot(groundwater, aes(x = month, y = resid)) +
  geom_point() +
  geom_line(alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Month", y = "Residual (m)",
       title = "What HolmesCo's model threw away")

summary(gw_model)$r.squared

R² = 0.489. HolmesCo transcribed their own number accurately, and everything else about the sentence is wrong.

The residual plot shows two clear structures:

a 12-month oscillation — recharge is seasonal, and lags rainfall by a few months, which a same-month model cannot represent;
a slow downward drift across the five years, which in the full report is a “temporary pumping phase” mentioned nowhere near the conclusion.

So 51% of the variation is unexplained, and that 51% is not noise. It is two named physical processes. “Dominant control” is a claim about rank, and this model has never compared rainfall against the two competitors sitting in its own residuals.

This is the most useful habit in the week: a model’s residuals are where its next version lives.

Extension: the diagnostics R gives you for free

Optional. ?plot.lm, ?acf, ?cooks.distance.

plot(gw_model) draws four standard diagnostics. Work out what each is for, and which of the four would have caught the seasonal cycle. (One of them will not, because it does not know the order of the observations — and that is exactly why a residuals-against-time plot has to be drawn by hand.)
acf(residuals(gw_model)) plots autocorrelation: how strongly each residual predicts the next. A spike at lag 12 is a year. Find it.
Cook’s distance measures how much each single observation moves the fit. Which month has the most influence, and would you have guessed it from the scatter plot?
Now draw residuals against fitted values instead of time. It looks much healthier. Write one sentence explaining how a model can pass that diagnostic and still be badly wrong.

Quick check

HolmesCo’s R² is 0.489. What percentage of the variation in groundwater level does their model fail to explain? Print it.

Exercise 4: Put the missing physics back

You can see a season and a trend in the residuals. Add them to the model and find out how much they were worth.

A trend is easy: month is already a column, and adding it lets the model drift. A 12-month cycle needs a pair of terms — sin(2 * pi * month / 12) and cos(2 * pi * month / 12) — which between them can represent a yearly wave of any size and any phase. One term alone fixes the phase, and the phase is exactly what you do not know.

Fit it and print the new R².

Hint 2

The shape of it:

fit <- lm(outcome ~ a + b + sin(2 * pi * t / 12) + cos(2 * pi * t / 12),
          data = df)
summary(fit)
summary(fit)$r.squared

Why a pair rather than one term: a sine alone peaks in a fixed month. Weighted together, a sine and a cosine can peak in any month, so the model can find the phase rather than being told it. That trick is worth remembering for any cyclical data — tides, diurnal temperature, seasonal demand.

Hint 3

gw_better <- lm(
  gw_level_m ~ rainfall_mm + ______ +
    sin(2 * pi * month / 12) + ______(2 * pi * month / 12),
  data = groundwater
)
summary(gw_better)

summary(gw_better)$______

Solution

gw_better <- lm(
  gw_level_m ~ rainfall_mm + month +
    sin(2 * pi * month / 12) + cos(2 * pi * month / 12),
  data = groundwater
)
summary(gw_better)

summary(gw_better)$r.squared

R² goes from 0.49 to 0.88. Most of what HolmesCo reported as unexplained variation was a seasonal cycle and a downward trend, both visible in their own residuals.

The coefficients matter more than the R². The month term is −0.017 with p = 0.0001: the water table is dropping by about 0.017 m per month, a metre over the five years of record. That is the single most consequential number in the dataset, it is nowhere in HolmesCo’s report, and their model could not have found it because they never gave it a way to.

Two cautions now that you have a better model.

A rising R² is not proof. Adding predictors can only increase it, so a higher R² is evidence only when the new terms have a reason to be there. Here they do — seasonal recharge and pumping are physical processes someone named in advance.

And “dominant control” is still not established. Rainfall is still significant in the improved model. What has changed is that it now has competitors, and you can see their sizes. The original claim was not false because rainfall does not matter; it was false because nothing had been compared with anything.

Extension: how would you know you had gone too far?

Optional. ?anova, ?AIC, ?step, ?summary.lm for adjusted R².

Compare the two models with anova(gw_model, gw_better). What is the null hypothesis of that test, and why can you only use it on models where one is nested inside the other?
Compare adjusted R² instead of R². Then add three columns of pure random noise as predictors and watch the two diverge.
AIC() penalizes complexity differently again. Rank all your models by AIC and by adjusted R². Do they agree — and what would you do if they did not?
step() will search for a model automatically. Run it, then write down why letting an algorithm choose your predictors and then reporting the p-values it produces is a form of the multiple comparisons problem from Week 7.

Exercise 5: The model that fits better and predicts worse

A more complex model always fits the data you already have at least as well. That is a mathematical fact about fitting, not evidence about the world.

Fit a 15-term polynomial to the solar data alongside the straight line, draw both, then ask each one what a panel does at 45 °C — just outside the range of the data.

Hint 2

The shape of it:

simple <- lm(y ~ x, data = df)
overfit <- lm(y ~ poly(x, 15), data = df)

grid <- data.frame(x = seq(3, 44, by = 0.2))
grid$simple <- predict(simple, newdata = grid)
grid$poly15 <- predict(overfit, newdata = grid)

ggplot(df, aes(x, y)) +
  geom_point() +
  geom_line(data = grid, aes(y = simple, colour = "Linear")) +
  geom_line(data = grid, aes(y = poly15, colour = "Polynomial")) +
  labs(colour = "Model")

predict(overfit, newdata = data.frame(x = 45))

The grid column has to be named exactly as the predictor is named in the model, or predict() cannot find it. You will also want coord_cartesian(ylim = c(200, 300)), because otherwise the curve leaves the plot and takes the data with it.

Hint 3

simple <- lm(output_w ~ temperature, data = solar)
overfit <- lm(output_w ~ ______(temperature, 15), data = solar)

temp_grid <- data.frame(temperature = seq(3, 44, by = 0.2))
temp_grid$simple <- predict(simple, newdata = temp_grid)
temp_grid$poly15 <- predict(______, newdata = temp_grid)

ggplot(solar, aes(temperature, output_w)) +
  geom_point() +
  geom_line(data = temp_grid, aes(y = simple, colour = "Linear")) +
  geom_line(data = temp_grid, aes(y = poly15, colour = "Polynomial (15)")) +
  coord_cartesian(ylim = c(200, 300)) +
  labs(x = "Temperature (°C)", y = "Output (W)", colour = "Model")

predict(overfit, newdata = data.frame(temperature = ______))

Solution

simple <- lm(output_w ~ temperature, data = solar)
overfit <- lm(output_w ~ poly(temperature, 15), data = solar)

temp_grid <- data.frame(temperature = seq(3, 44, by = 0.2))
temp_grid$simple <- predict(simple, newdata = temp_grid)
temp_grid$poly15 <- predict(overfit, newdata = temp_grid)

ggplot(solar, aes(temperature, output_w)) +
  geom_point() +
  geom_line(data = temp_grid, aes(y = simple, colour = "Linear")) +
  geom_line(data = temp_grid, aes(y = poly15, colour = "Polynomial (15)")) +
  coord_cartesian(ylim = c(200, 300)) +
  labs(x = "Temperature (°C)", y = "Output (W)",
       colour = "Model",
       title = "Which would you trust at 45 °C?")

predict(overfit, newdata = data.frame(temperature = 45))

The polynomial has R² = 0.60 against the straight line’s 0.19 — three times the fit by that measure, and you need coord_cartesian() to see the data at all, because the curve leaves the plot near the edges.

At 45 °C the straight line predicts 241.6 W. The polynomial predicts about 67,800 W, from a panel rated at 250. Nothing objected.

Three things to take from this:

Fit and prediction are different. R² measures agreement with data you have already seen — the one thing you never needed a model for.
Extrapolation fails silently. The polynomial is unreliable inside the data too; outside it, it is absurd. A matter of degree, not kind.
This is HolmesCo’s groundwater extrapolation in a costume: a line fitted to five years, extended forward to predict an empty aquifer by
1. The model does not know what it is modelling.

The straight line misses plenty. It is still the one you would stake a recommendation on.

Extension: find where it breaks

Optional.

Refit at degrees 2, 5, 8 and 15, predicting at 45 °C each time. Where does the answer stop being physically plausible? There is no sharp threshold, which is itself the point.
Split the data: fit on a random 35 panels and measure the error on the 15 you held back. Do that for each degree and plot error against complexity. The curve should fall and then rise — that U is the whole of model selection in one picture.
Repeat task 2 with a different random split. How much does the best degree move? With 50 observations, rather a lot — which is why serious model selection uses many splits rather than one.
Fit the model physics actually suggests: output against temperature - 25, so the intercept is the rated output at the reference temperature. Compare its R² with the polynomial’s, then say why “it fits worse” is not the objection it sounds like.

Save your work

Copy the code you’re most proud of into your week8.R file. Commit and push via GitHub Desktop. Write a commit message that describes what you learned — not just “week 8”.