Skip to main content

Module 8 lab: what belongs in the model

Work through this at your own pace. Each step ends with a check – something you should be able to see on your screen before moving on. If a check fails, that is the point to put your hand up rather than pushing on.

This lab is all in R. Excel will fit a multiple regression through the Analysis ToolPak, but it will not let you add a variable and watch the other coefficients move, which is the whole point of today.

The shape of the day: add predictors one at a time and watch what happens to the ones already there, then work out which direction a variable you don’t have would push your answer, then stop trusting the fit and measure the model on data it has never seen.

Plan on about an hour and a quarter for the eleven steps, and half an hour more for the Practice boxes. It covers 26  Multivariate Linear Models, 27  Omitted Variable Bias and 28  Predictions, Errors and Model Quality. If you run out of time, Step 8 onwards is the part to finish at home.

Step 1: Get the data and set up the script

Open Positron and open your AREC_261 folder, the one from the Module 2 lab. Every read_csv("data/...") below is relative to the folder Positron has open.

No new packages today. Download these into data/:

Make a new script in code/ called module8_lab.R:

# ---
# Title: Module 8 lab
# Author: Your Name
# Date: 2026-11-04
# Description:
#   Building a multiple regression, signing omitted variable bias,
#   and judging a model out of sample.
# ---

# Load packages
library(tidyverse)
rm_weather <- read_csv("data/rm_canola_weather.csv")

dim(rm_weather)

Check: 2,950 rows and 20 columns.

Step 2: Look at the data before modelling it

# Canola is blank where the crop was not grown in that RM that year.
canola <- filter(rm_weather, !is.na(Canola))

nrow(canola)

Check: 2,904 rows, so 46 RM-years dropped.

Those 46 are not an error and they are not missing at random – they are RM-years where nobody grew canola, which is itself informative. Everything below is a statement about RM-years where canola was grown.

The soil zone needs to be a factor, and the order matters because R compares every level against the first one:

canola <- mutate(canola,
                 soil_zone = factor(soil_zone,
                                    levels = c("Brown", "Dark Brown", "Grey/Black")))

table(canola$soil_zone)

Check: 860 Brown, 940 Dark Brown, 1,104 Grey/Black.

Set the levels deliberately. Left alone R sorts them alphabetically, which here would make Brown the baseline anyway – but on a variable like c("High", "Low", "Medium") alphabetical order is nonsense and you would be comparing everything against “High” without noticing.

Step 3: Start with one predictor

rain_only <- lm(Canola ~ Precip_May_Aug_mm, data = canola)
summary(rain_only)

Check: the rainfall slope is 0.057, its p-value is below 0.001, and the adjusted \(R^2\) is 0.17.

Read that slope in units before going on. It is bushels per acre per millimetre, so 100 mm more rain goes with about 5.7 bu/ac more canola. A coefficient you cannot say out loud in units is a coefficient you are about to misinterpret.

Step 4: Add a second predictor and watch the first one move

both <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30, data = canola)
summary(both)

Check: rainfall has dropped to 0.005 and hot days comes in at -0.874. Adjusted \(R^2\) rises from 0.17 to 0.51.

Rainfall just lost about nine tenths of its slope because a second variable joined the model. Nothing about the rain data changed. Here is why:

cor(canola$Precip_May_Aug_mm, canola$Days_Above_30)

Check: -0.54.

Wet years are cool years. With hot days left out, the rainfall coefficient was being paid for two jobs – the moisture itself, and the heat that was absent whenever it rained. Adding hot days takes the second job away from it.

This is the mechanism behind the whole module, and it is worth being clear that neither coefficient is “the” answer. Which one you want depends on the question. If you are asking what a wet year looks like, 0.057 is the honest number – wet years really do come with 5.7 bu/ac per 100 mm, heat and all. If you are asking what water does, 0.005 is closer, because it holds the heat fixed.

Practice

Fit Canola ~ Days_Above_30 on its own. Compare the hot-day slope with the -0.874 you just got. Which direction did it move, and does the sign of the correlation above explain it?

On its own the hot-day slope is -0.898, against -0.874 alongside rainfall. Slightly more negative, and the sign of the correlation does explain the direction: hot years are dry years, so with rainfall omitted the hot-day coefficient picks up a little of the drought as well.

What is worth noticing is how small the move is – about 3%, against the 90% collapse rainfall suffered in the other direction. The asymmetry is the lesson. Omitting a variable biases what remains by its own effect times how strongly the two move together, so a confounder only does real damage when both parts are substantial. Rainfall has almost no effect of its own once heat is held fixed, so leaving it out costs the hot-day coefficient almost nothing. Reversing the two does not reverse the size of the problem.

Step 5: Add latitude

with_lat <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude,
               data = canola)
summary(with_lat)

Check: latitude is 1.01, hot days has softened from -0.874 to -0.791, and rainfall has risen to 0.008.

Latitude says an RM a degree further north yields about a bushel more. The growing area spans about five degrees, so taken literally that is 5 bu/ac from the bottom of the province to the top.

Notice the hot-day coefficient softening. The south is hotter and worse for canola for reasons that are not heat, and the hot-day count had been carrying some of that. Latitude takes part of it back.

Step 6: Add a categorical predictor

with_both <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,
                data = canola)
summary(with_both)

Check: two soil-zone rows, Dark Brown at 2.73 and Grey/Black at 4.04. Latitude has collapsed from 1.01 to 0.26 and its p-value has risen to 0.099.

Three things to take from this output.

There is no row for Brown. R turned the three-level factor into two columns and held Brown out as the baseline, so every zone coefficient is a comparison against Brown. Dark Brown yields 2.73 bu/ac more than Brown; Grey/Black yields 4.04 more. Neither number is that zone’s yield.

The Grey/Black-to-Dark-Brown comparison is not in the table. You get it by subtracting: \(4.036 - 2.725 = 1.31\) bu/ac.

And latitude has almost vanished. That is the finding, not a failure. Northern RMs are mostly dark brown and grey/black RMs, so what looked like an effect of being further north was largely an effect of the soil underneath. With soil in the model, latitude has little left to explain.

Check: you can say in one sentence why latitude’s p-value went from below 0.001 to 0.099 without any new data arriving.

Step 7: Predict, and notice what prediction will do for you

new_rm <- data.frame(Precip_May_Aug_mm = 250,
                     Days_Above_30     = 8,
                     latitude          = 52,
                     soil_zone         = "Dark Brown")

predict(with_both, newdata = new_rm)

Check: 38.85 bu/ac.

Several rows at once is the usual way to compare scenarios:

scenarios <- data.frame(
  Precip_May_Aug_mm = c(150, 250, 250),
  Days_Above_30     = c( 20,   8,   8),
  latitude          = c(49.5, 52, 52),
  soil_zone         = c("Brown", "Dark Brown", "Grey/Black"))

scenarios$predicted <- predict(with_both, newdata = scenarios)
scenarios

Check: 25.95, 38.85 and 40.17 bu/ac.

The last two rows differ only in soil zone, and the 1.31 bu/ac between them is exactly the difference of the two zone coefficients from Step 6.

Now ask the model something it has no business answering:

silly <- data.frame(Precip_May_Aug_mm = 1000,
                    Days_Above_30     = 8,
                    latitude          = 52,
                    soil_zone         = "Dark Brown")

predict(with_both, newdata = silly)
max(canola$Precip_May_Aug_mm)

Check: a confident number around 44.6 bu/ac, against a wettest-ever RM-year of 470 mm.

No warning, no NA. The model extrapolated to more than double the wettest observation in the data and returned a number to four decimal places. predict() will not tell you when you have left the data behind – that check is yours.

Practice

Try soil_zone = "dark brown" in lower case, then drop latitude from the data frame entirely. What does each do?

The lower-case label errors – R matches it against the factor levels and does not find it. Dropping a predictive column errors too: predict() needs a column for every variable in the model and will not assume a value.

Both are good failures. The bad failure is the one above, where the answer is wrong and nothing complains.

Step 8: Sign the bias from a variable you do not have

Switch datasets. This one is about whether education raises farm income, and it is the clearest demonstration of omitted variable bias because here you do have the omitted variable.

farm_income_data <- read_csv("data/farm_income.csv")

short <- lm(farm_income ~ years_educ, data = farm_income_data)
coef(short)

Check: 7.81 – a year of post-secondary education goes with about $7,811 more farm income.

Students who go further in school were, on average, stronger students beforehand, and that would have shown up in their income regardless. grade12 is a measure of them before university:

long <- lm(farm_income ~ years_educ + grade12, data = farm_income_data)
coef(long)

Check: education has fallen to 4.38 and grade 12 average comes in at 1.93.

Nearly half the apparent return to education was really the students who were going to earn more anyway. Now close the loop with the auxiliary regression – the omitted variable on the included one:

aux <- lm(grade12 ~ years_educ, data = farm_income_data)
coef(aux)

Check: 1.78 – each extra year of schooling goes with a grade 12 average 1.78 points higher.

Those three numbers are not independent. Equation 27.1 says the short slope equals the long slope plus the bias, where the bias is the omitted variable’s own effect times the auxiliary slope:

# the short regression's slope, reconstructed from the long one
coef(long)["years_educ"] + coef(long)["grade12"] * coef(aux)["years_educ"]

coef(short)["years_educ"]

Check: both print 7.8114. \(4.38 + 1.93 \times 1.78 = 7.81\), exactly.

That is an algebraic identity, not an approximation – it holds to every decimal place R will print. The short regression was never a different estimate of the same thing; it was the long one plus a bias term you can now name.

Check: before moving on, say which cell of the Which direction? table this is. Grade 12 raises income, so \(\beta_2 > 0\); stronger students go further in school, so the correlation is positive; top-left, positive bias, the estimate is too big. You could have called the direction before fitting anything.

Practice

Fertilizer and yield, omitting soil quality. Farmers put more fertilizer on their better land. Which way is the fertilizer coefficient biased, and would you quote it as an upper or a lower bound on the real effect?

Better land gets more fertilizer (positive correlation) and better land yields more on its own (\(\beta_2 > 0\)). Both terms positive, so the bias is positive and the coefficient overstates the effect – quote it as an upper bound.

That bounding move is most of what you can do when the omitted variable is genuinely unavailable. It is weaker than an unbiased estimate and much stronger than nothing: “the effect is at most this large” often settles the decision on its own.

Step 9: Stop trusting the fit

Back to the canola model. Every number so far came from the same 2,904 rows the model was fitted on.

# in-sample accuracy
err <- resid(with_both)

mean(abs(err))          # MAE
sqrt(mean(err^2))       # RMSE

Check: MAE 5.70 and RMSE 7.19 bu/ac.

RMSE above MAE is the normal relationship – squaring punishes the big misses harder, so the gap between them tells you how much of your error is concentrated in a few bad years.

Now watch the fit improve for no reason at all. Take 150 rows to fit on, hold back 50, and bolt on 100 columns of pure random noise:

set.seed(261)
rows  <- sample(nrow(canola), 200)
train <- canola[rows[1:150], ]
test  <- canola[rows[151:200], ]

# 100 columns of random numbers, unrelated to anything
add_junk <- function(df, seed) {
  set.seed(seed)
  junk <- as.data.frame(matrix(rnorm(nrow(df) * 100), nrow = nrow(df)))
  names(junk) <- paste0("junk", 1:100)
  cbind(df, junk)
}
train <- add_junk(train, 1)
test  <- add_junk(test, 2)

real_model <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,
                 data = train)
junk_model <- lm(as.formula(paste("Canola ~ Precip_May_Aug_mm + Days_Above_30 +",
                                  "latitude + soil_zone +",
                                  paste0("junk", 1:100, collapse = "+"))),
                 data = train)

accuracy <- function(model, data) {
  err <- data$Canola - predict(model, data)
  c(R2   = 1 - sum(err^2) / sum((data$Canola - mean(data$Canola))^2),
    MAE  = mean(abs(err)),
    RMSE = sqrt(mean(err^2)))
}

round(rbind(real = accuracy(real_model, train),
            junk = accuracy(junk_model, train)), 3)

Check: the junk model wins on all three. \(R^2\) 0.815 against 0.452, MAE 3.54 against 5.98, RMSE 4.44 against 7.63.

A hundred columns of rnorm() just beat the weather. Every measure you have says the junk model is better, and every one of those columns is random numbers.

Step 10: Judge it on data it has never seen

The 50 rows held back in Step 9 played no part in fitting either model.

round(rbind(real = accuracy(real_model, test),
            junk = accuracy(junk_model, test)), 3)

Check: the real model holds at RMSE 7.32. The junk model goes to 13.99, and its \(R^2\) is -1.211.

The verdict reverses completely. The real model’s out-of-sample RMSE (7.32) is close to its in-sample RMSE (7.63), which is what an honest model looks like. The junk model roughly triples.

A negative \(R^2\) is worth a moment, because the in-sample version cannot produce one. It is still \(1 - \text{SSE}/\text{SST}\); the model’s squared errors simply come to more than the total variation around the hold-out mean, so the ratio passes one. In plain terms the junk model predicts these 50 RM-years worse than someone who knew nothing but their average.

In a real forecast you would not know the test mean in advance, so here is the benchmark you could actually have used:

# guess the training mean for every held-out RM
sqrt(mean((test$Canola - mean(train$Canola))^2))

Check: 9.48 – and the junk model’s 13.99 loses to that too.

Check: you can state what the 100 junk columns bought. They did not teach the model about canola; they let it memorise 150 particular rows, and that memory is worth less than nothing anywhere else.

Practice

Re-run Step 9 and Step 10 with a different seed in sample() – try 1, then 7. Do the in-sample and out-of-sample verdicts ever disagree about which model is better?

The in-sample comparison always favours the junk model and the out-of-sample comparison always favours the real one. The exact numbers move; the direction does not.

That is the point of running it more than once. The gap is not a quirk of seed 261 – adding variables cannot worsen in-sample fit, so the junk model wins there by construction, every time.

Step 11: One model, written up honestly

Fit the model you would actually defend, on the full data, and say what it does and does not support:

final <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,
            data = canola)

summary(final)

Then answer these five, in writing, in your script as comments:

  1. What does each coefficient mean, in units, holding what constant?
  2. Which predictor would you drop first, and on what grounds other than its p-value?
  3. Name one variable that is not in this model, is likely to matter, and is likely correlated with rainfall or heat. Which way does omitting it bias those coefficients?
  4. Over what range of rainfall is this model usable? What happens outside it?
  5. If you had to predict 2027 canola yields with it, what would you want to check first?

Check: your answer to 3 names a direction, not just a variable.

Question 3 is the one with a real answer worth comparing. Soil moisture at seeding is a common one: it raises yield and it correlates with growing-season rain, so both terms are positive and the rainfall coefficient is biased up. Management – seeding date, variety, fertility – is the other obvious candidate, and there the correlation with weather is weaker, which is exactly the argument for worrying about it less.

Practice on your own

The test bank for this module is here, and it is the best use of the remaining time.

Two other things if you want a break from the bank.

Take the wheat data, sask_wheat_2025.csv, and regress yield_bu_ac on protein_pct. You will get a slope of about -3.64, which says that more protein goes with less yield. Then add soil_zone and variety: the protein slope falls to about -1.37 but stays significant. The protein-yield trade-off is real agronomy, and roughly two thirds of the raw slope was which varieties get grown where.

Then go back to the canola model and add Population_2021. It comes in statistically significant at the 5% level, with a p-value around 0.022, and the adjusted \(R^2\) barely moves – 0.5307 against 0.5300. Write down the story you would have to believe for RM population to cause canola yields, and then the story you actually believe. This is the clearest example in the course of significance and relevance being different things.

Before you leave

Save your script. You should have one that fits a regression in four stages and tracks what happens to each coefficient as predictors arrive, reconstructs a short regression’s slope from the long one plus a bias term, and judges two models both in and out of sample.

The habit to keep is the one from Step 10. A model’s fit on the data it was built from is not evidence that it predicts. Hold something back, or say plainly that you did not.