This bank holds 30 questions, in six sections. The first three sections are conceptual (no coding) and the next three involve coding.
The module test will contain two questions from the first three sections, and two questions from the last three sections. The first two questions will be answered on Canvas (with no help from AI). Then the last two questions can be answered using AI to help you. I will change the wording of these questions slightly – a different crop, a different weather column, or a different file – though the concepts will be exactly the same.
The following is output from a model with house sale prices (dollars) as the response variable and floor area (square feet), age of the house (years) and whether it has a garage as explanatory variables.
(a) Interpret the coefficients on sqft and garageYes, each in a sentence with units.
(b) Two houses are the same age and both have a garage. One is 200 sq ft larger. How much more would you expect it to sell for?
(c) The intercept is $52,012. Does this mean a brand new house with no floor area and no garage is worth $52,012?
(d) The \(R^2\) is 0.709. Explain what this means.
Answer
(a) Comparing two houses of the same age with the same garage status, one square foot more floor area goes with $94.50 more in sale price. Comparing two houses of the same age and size, the one with a garage sold for $12,872 more.
(b)\(200 \times 94.504 = \$18,901\) (to the nearest dollar, $18,900.80).
(c) No. No house in the data has zero floor area, so the intercept is far outside the range of the data – it is where the fitted plane happens to cross zero, not a measurement of anything real. It is extrapolation.
(d) The model accounts for about 71% of the variation in sale price across these houses, leaving 29% unexplained.
What it does not tell you is how far out a price would typically be, because it is a proportion rather than an amount. Two models with the same \(R^2\) can be wrong by very different numbers of dollars, depending on how much prices varied to begin with. For pricing a house you would want the error in dollars – the MAE or the RMSE – which is the measure in the units of the decision.
Question 2
The following is output from a model with milk yield per cow (litres per day) as the response variable and feed intake (kg of dry matter per day), the age of the cow (years) and days in milk as explanatory variables.
(a) Interpret each of the three slope coefficients in a sentence with units.
(b) A cow is 100 days further into her lactation (e.g., days in milk) than another, and the two match on feed and age. What does the model say about the difference in their milk yield?
(c) The farm manager reads the feed coefficient as “feeding a cow one more kg of dry matter will produce one more litre of milk.” Is there anything wrong with this logic?
(d) Predict the milk yield for a cow eating 24 kg of dry matter per day, aged 4 years, at 150 days in milk.
Answer
(a) Comparing cows of the same age at the same stage of lactation, one more kg of dry matter per day goes with 0.936 more litres of milk per day. Comparing cows on the same feed at the same days in milk, each additional year of age goes with 1.022 more litres per day. Comparing cows of the same age on the same feed, each additional day into the lactation goes with 0.022 fewer litres per day.
(b)\(100 \times (-0.022205) = -2.22\) litres per day, so the cow further into her lactation is giving about 2.2 litres per day less.
(c)Will is a causal claim about what happens if you intervene, and the model only describes a comparison between cows that already differ. Cows that eat more may differ in other ways that also raise milk – better health, better genetics, a stage of lactation not fully captured by days in milk. The model says cows eating more give more; it does not promise that feeding a given cow more would produce that result.
The following is output from a model with grain truck fuel use (litres per 100 km) as the response variable and load weight (tonnes), average speed (km/h) and whether the route is gravel or paved as explanatory variables (i.e., surfacePaved equals one if the route is paved and zero if it is gravel).
(a) Interpret the coefficient on surfacePaved in a sentence with units.
(b) The speed coefficient is negative. Does this mean driving faster always saves fuel? What would you want to check before believing it?
(c) A fully loaded truck carries 40 tonnes. Using the model, how much more fuel per 100 km does a full load use than an empty one, holding speed and surface fixed?
(d) Name a variable that is missing from this model. Would leaving it out worry you?
Answer
(a) Comparing two hauls with the same load at the same average speed, a paved route used 3.22 fewer litres per 100 km than a gravel route.
(b) It means that within the range of speeds in this data, and holding load and surface constant, faster hauls used slightly less fuel per 100 km. Two things to check. First, the range of speeds observed – if they run from 60 to 95 km/h, the model says nothing about 130 km/h, and fuel use per km rises again at high speed. Second, whether speed is standing in for something else: faster hauls are probably highway hauls with fewer stops, so the coefficient may be picking up road type and traffic rather than speed itself.
(c)\(40 \times 0.8412 = 33.6\) litres more per 100 km.
(d) Answers will vary. Reasonable candidates: truck age or engine size – older trucks burn more, and if older trucks also haul lighter loads this would bias the load coefficient, so it would worry me. Wind or weather – it matters for fuel use, but there is no obvious reason it would correlate with load weight or surface, so it would mostly just add noise rather than bias. The distinction to draw is between a variable that affects the response and is correlated with a predictor (a problem) and one that affects the response but is unrelated to the predictors (not a bias problem).
Question 4
The following is output from a model with the price received for a load of canola ($/tonne) as the response variable and the protein content, the dockage percentage, and the month of delivery (from 1 to 12 with January =1 and December = 12) as explanatory variables (fitted on 240 loads).
(a) Interpret the dockage coefficient in a sentence with units.
(b)month is entered as a number from 1 to 12. What does this imply about how price moves through the year, and why is that a questionable way to model a month?
(c) The \(R^2\) is 0.29. Is this a bad model?
(d) Before relying on this model to price a load, what else would you want to know?
Answer
(a) Comparing two loads delivered in the same month with the same protein, one percentage point more dockage went with $9.87 per tonne less in price.
(b) Entered as a number, the coefficient claims price rises by a constant $1.04 per tonne for every month later in the year – a straight line from January to December. That is a strong assumption and almost certainly wrong: prices move with harvest and storage cycles, not linearly with the calendar, and it forces December and January to be 11 months apart when they are adjacent. Entering month as a categorical variable would let each month have its own effect.
(c) Not necessarily. An \(R^2\) of 0.29 means most of the variation in price is unexplained, which is unsurprising – price depends heavily on the world market on the day of delivery, which is not in the model. Whether it is a bad model depends on what it is for. For understanding how dockage is discounted, it may be fine. For forecasting the price of the next load, it is close to useless on its own.
(d) The typical size of a miss in dollars per tonne – the MAE or the RMSE – since that is the number in the units of the decision, and the \(R^2\) on its own does not give it. Also the range of protein and dockage in the data, so a prediction is not an extrapolation beyond it; and whether the model has been tested on loads it was not fitted on, since the fit reported here is on the rows that chose the coefficients.
Question 5
The following is output from a model with crop insurance claim size ($ per acre) as the response variable and three explanatory variables: drought index, acres insured, and years of farming experience (fitted on 1,450 claims).
(b) The coefficient on acres_insured is 0.0021 with a standard error of 0.0094 and a p-value of 0.823. What does this tell you about the relationship between farm size and claims per acre?
(c) Would you be comfortable telling an insurer that taking on more experienced farmers as customers would reduce claims? Explain.
Answer
(a) Comparing two claims with the same drought index on the same insured acreage, each additional year of farming experience goes with $0.63 per acre less in claim size.
(b) The estimate is about a fifth of its own standard error, and the p-value of 0.823 says zero is entirely consistent with what we observe. The data cannot distinguish this coefficient from zero or even establish its sign.
That is not the same as showing the effect is zero – a non-significant coefficient is a statement about the evidence, not about the world. The estimate is also tiny in dollar terms, two tenths of a cent per acre insured, so it would not matter much either way.
(c) No. The model describes a comparison between farmers who already differ, and experience is bundled with many other things – better land acquired over time, more conservative cropping, better equipment. Any of those could be doing the work. The coefficient is a fact about the comparison, not a prediction of what would happen if a given farm gained experience.
Section 2 — Omitted variable bias
Question 1
(a) We want to estimate \(\beta_1\) in \(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon\) but we lack data on \(x_2\). Instead we estimate \(y = \tilde{\beta}_0 + \tilde{\beta}_1 x_1 + v\) instead. Write down the expression that relates \(\tilde{\beta}_1\) to \(\beta_1\) and explain what the bias term is in this expression.
(b) State the two conditions under which \(\tilde{\beta}_1\) equals \(\beta_1\).
Answer
(a)\(\tilde{\beta}_1 = \hat{\beta}_1 + \hat{\beta}_2 \tilde{\delta}_1\). The bias is the product \(\hat{\beta}_2 \tilde{\delta}_1\). The first piece, \(\hat{\beta}_2\), is the effect of the omitted variable on the response. The second, \(\tilde{\delta}_1\), is the slope from regressing the omitted variable on the included one – how strongly the two predictors move together.
(b) Either \(\hat{\beta}_2 = 0\), meaning the omitted variable has no effect on the response; or \(\tilde{\delta}_1 = 0\), meaning the omitted variable is uncorrelated with the included one. If either is zero the product is zero and there is no bias.
Question 2
For each of the following, state how omitting the stated variable would bias the estimate of the included variable.
(a) Estimating the effect of a new seed variety on yield, omitting the quality of the land it was sown on. The new variety is expensive, so growers tried it on their best fields first.
(b) Estimating the effect of a herbicide application on yield, omitting weed pressure.
Answer
(a) Better land raises yield, so \(\beta_2 > 0\). The variety went on the best fields, so \(\text{Corr}(\text{new variety},\ \text{land quality}) > 0\). Both terms are positive: positive bias, and the variety looks better than it is.
(b) Weed pressure cuts yield, so \(\beta_2 < 0\). Herbicide is sprayed where the weeds are, so the correlation is positive. A positive times a negative is negative: negative bias – spraying can look like it reduces yield when it is helping.
Question 3
A regression of farm income on years of education alone gives a slope of 7.81 ($000s). Adding grade 12 average gives a slope of 4.38 on education and 1.93 on grade 12 average. Regressing grade 12 on years of education gives a slope of 1.78.
(a) Verify that these four numbers are consistent with the omitted variable bias formula. Show the arithmetic.
(b) A reader says “so the real effect of education is $4,380, and the other $3,430 is fake.” Is that the right way to describe it?
Answer
(a)\(\tilde{\beta}_1 = \hat{\beta}_1 + \hat{\beta}_2 \tilde{\delta}_1 = 4.38 + 1.93 \times 1.78 = 4.38 + 3.43 = 7.82\), which reproduces the education-only slope of 7.81 to within rounding. All three inputs are rounded to two decimals; the full fitted coefficients give 7.811.
(b) Not quite. The $3,430 is not fake – it is a real difference in income between people with more and less education. It is simply not attributable to the education itself: it is the part that reflects the fact that people with more education also had better grade 12 marks, and those marks are worth $1,928 a point. The education-only regression correctly describes the raw comparison; it just cannot be read as the effect of education.
The $4,380 should not be called the “real effect” either. It is the estimate adjusted for grade 12 marks, and that is all adding one control buys. Family finances, access to a university, and what someone intended to do with the degree are all still missing, and any of them could be doing part of the work.
Question 4
A newspaper reports: “Farms that hired an agronomist averaged 6 bu/ac more canola than those that did not.”
(a) Write down the regression the newspaper has effectively run.
(b) Identify a variable that causes bias by its omission and work out the direction of the bias.
(c) The agronomy firm quotes this result in its advertising. Is this false advertising?
Answer
(a) The newspaper’s comparison is \(\text{yield} = \tilde{\beta}_0 + \tilde{\beta}_1 \text{agronomist} + v\), with \(\tilde{\beta}_1 = 6\). What it needed is \(\text{yield} = \beta_0 + \beta_1 \text{agronomist} + \beta_2 x_2 + \varepsilon\), where \(x_2\) is whatever else differs between the two groups of farms.
(b) Management quality is the most likely candidate – farms that hire an agronomist are the ones already paying attention to agronomy. Good management raises yield (\(\beta_2 > 0\)) and is positively correlated with hiring one, so the bias is positive and the 6 bu/ac overstates what the agronomist contributed. Farm size and soil quality would work similarly.
(c) Not false, but easy to state in a misleading way. What the data support is that farms using an agronomist yielded 6 bu/ac more on average – stated as a comparison, not an effect. They could not honestly claim that hiring one would raise a given farm’s yield by 6 bu/ac, because the farms being compared differ in ways the figure does not account for.
Question 5
A provincial report finds that fields enrolled in a soil-conservation programme had 4% more organic matter than fields that were not enrolled, and concludes that the programme built soil.
(a) Name two things that could differ between enrolled and unenrolled fields other than the programme itself, and say for each whether it would push the 4% up or down.
(b) The report’s authors add soil zone and field size as controls and the gap falls to 1.5%. What does that tell you about the original figure?
(c) A reviewer says the remaining 1.5% must be the true effect, since the obvious confounders have been controlled for. Why is that too strong?
Answer
(a) Answers will vary. Operator attitude to soil – growers who already care about soil enrol, and they were already managing for organic matter, so this pushes the 4% up. Land type – if the programme targets marginal or erodible land, which tends to be lower in organic matter, this pushes the figure down. A good answer names the variable, says which way it affects organic matter, and says which way it relates to enrolling.
(b) Soil zone and field size together accounted for well over half of the original gap, so enrolment was substantially confounded with where and what the land is. It also tells you the remaining variables are likely to matter too, since the two easiest ones to measure moved the estimate this much.
(c) Controlling for two variables removes the bias from those two and nothing else. The most worrying confounder here – the grower’s own management – is unmeasured, and it is correlated with enrolling and with organic matter. Controlling for what you happen to have is not the same as controlling for what matters.
(b)\(\hat{\varepsilon} = 38.2 - 43.19 = -4.99\). The RM yielded about 5 bu/ac less than the model predicted – a negative residual means the actual value came in below the fitted one.
(c) The fitted value \(\hat{y}_i\) is the model’s prediction for observation \(i\), obtained by putting that observation’s predictor values into the fitted equation. The residual \(\hat{\varepsilon}_i\) is the part the model did not capture. They are connected by \(y_i = \hat{y}_i + \hat{\varepsilon}_i\), or equivalently \(\hat{\varepsilon}_i = y_i - \hat{y}_i\).
Question 2
(a) Define the total sum of squares (SST) and the sum of squared errors (SSE), in words and in formulas.
(b) Write the formula for \(R^2\) in terms of SSE and SST, and explain why it equals 0 when SSE equals SST.
(c) A model has SSE = 3,187 and SST = 7,802. Compute \(R^2\) and state what it means in a sentence.
Answer
(a) SST is the total sum of squares, the total variability in the response measured against its own mean: \(\text{SST} = \sum (y_i - \bar{y})^2\). It is the error you would make predicting every observation with the mean. SSE is the sum of squared errors from the model: \(\text{SSE} = \sum (y_i - \hat{y}_i)^2\).
(b)\(R^2 = 1 - \text{SSE}/\text{SST}\). When SSE equals SST the model’s errors are exactly as large as the errors from guessing the mean, so the model has added nothing and the ratio is 1, giving \(R^2 = 0\).
(c)\(R^2 = 1 - 3187/7802 = 1 - 0.4085 = 0.5915\). The model accounts for about 59% of the variation in the response across these observations.
Question 3
A grain company models the cost of trucking a load to the elevator, in dollars per tonne, on the distance hauled (km), the weight of the load (tonnes) and whether the route is paved. Fitted on 240 hauls.
paved is 1 for a paved route and 0 for gravel. Distances run from 16 to 180 km and loads from 8.2 to 39.8 tonnes.
(a) Predict the cost per tonne for a 28-tonne load hauled 120 km over gravel. Show the arithmetic term by term.
(b) Predict how the cost would change for the same haul on a paved route. Even if you used a calculator to give you this answer, explain why you didn’t need one.
(c) Predict the cost for a 36-tonne load hauled 45 km over gravel. Explain the two reasons this is different from your answer in part a.
(b)\(11.15\) per tonne. Everything except the surface is identical, so the whole difference is the paved coefficient: switching the indicator from 0 to 1 subtracts 2.3698 and nothing else moves.
(c)\(9.5243 + 0.0849 \times 45 - 0.2213 \times 36 - 2.3698 \times 0 = 5.38\) per tonne. The gap from (a) is \(13.52 - 5.38 = 8.14\). Of that, 75 fewer kilometres is worth \(0.0849 \times 75 = 6.37\), and 8 more tonnes on the load is worth \(0.2213 \times 8 = 1.77\). Together \(6.37 + 1.77 = 8.14\).
Question 4
A canola model is fitted on RM-years, and three RMs from 2023 come out as follows.
RM
Actual yield
Fitted value
Bone Creek
21.4
33.9
Eldon
44.8
37.2
Mayfield
35.1
35.4
(a) Compute the residual for each RM.
(b) Which RM did the model predict the best, and which worst?
(c) After reviewing the data, Bone Creek’s weather that year looked ordinary. Give two reasons an RM’s yield can come in far below its fitted value even when the weather inputs are unremarkable.
(d) Across all 2,904 RM-years the residuals sum to approximately zero. Does that mean the model is unbiased? Explain.
Answer
(a) Bone Creek \(21.4 - 33.9 = -12.5\); Eldon \(44.8 - 37.2 = +7.6\); Mayfield \(35.1 - 35.4 = -0.3\).
(b) Mayfield best, missed by 0.3 bu/ac. Bone Creek worst, by 12.5 bu/ac – the crop came in 12.5 bushels an acre below what the weather and location predicted, which on a 33.9 bu/ac prediction is a miss of more than a third.
(c) Answers will vary. Anything that affects yield and is not a model input: hail or a localised storm; disease or an insect outbreak; late seeding; a frost outside the measured window; the weather station being far from the RM, so the recorded rainfall is not what the fields received. The common thread is that the residual collects everything the model does not measure.
(d) No. The residuals summing to zero is arithmetic, not evidence – least squares forces it whenever the model contains an intercept, on the rows it was fitted to, and it would hold just as exactly for a badly confounded model. (On held-out rows, or in a model without an intercept, there is no such guarantee.) Unbiasedness in the omitted variable bias sense is about whether a coefficient can be read as the effect of that variable, which the residuals say nothing about. A model can fit the data closely, have residuals that sum to zero, and still have every coefficient misleading.
Question 5
A researcher is building a model to predict canola yields for next season. She has a model with two weather variables (rainfall and heat) and two location variables (latitude and soil zone). She is debating adding two variables to the model: the average soil class in the RM and the number of school age children in the RM.
(a) Explain what adding each of these variables will do to the in-sample \(R^2\), MAE and RMSE?
(b) Do you have a hypothesis about what these variables will do to the out of sample MAE and RMSE?
(c) Would you suggest adding these variables to the model?
Answer
(a) Both variables behave the same way in sample. \(R^2\) cannot fall and RMSE cannot rise, because least squares could set the new coefficient to zero and recover the old fit. The guarantee allows equality – these cannot get worse, but they need not get better.
MAE carries no such guarantee and can increase, since least squares minimises squared error rather than absolute error.
(b) Out of sample is where the two variables differ. Soil class has a mechanism – soil productivity affects yield, and soil zone is only a coarse three-category version of it – so expect a modest improvement in MAE and RMSE. School age children have no mechanism, so expect nothing: a small move either way, meaning nothing.
(c) Add the soil class, leave the children out.
The children should stay out whatever the in-sample \(R^2\) does, since that improvement was guaranteed. Even a small out-of-sample gain would more likely be luck than signal on one test set. The cost of keeping it is a coefficient someone will eventually try to interpret.
Section 4 — Running a multiple regression in R
Question 1
Using RM yields and weather, one row per Saskatchewan RM per year from 2015 to 2024. Yields are in bu/ac, Precip_May_Aug_mm is total May-August rainfall in mm, and Days_Above_30 counts days hotter than 30 °C.
(a) Read the file, keep the rows with a spring wheat yield, and fit a model of spring wheat yield on growing-season rainfall (Precip_May_Aug_mm) alone. Report the slope and the \(R^2\).
(b) Add the count of days above 30 °C (Days_Above_30) and refit. Report both slopes and the new \(R^2\).
(c) Explain how and why the coefficient on rainfall changes from part a to part b.
Answer
library(tidyverse)wheat <-read_csv("data/rm_canola_weather.csv", show_col_types =FALSE) |>filter(!is.na(`Spring Wheat`))## (a) rainfall alonerain_only <-lm(`Spring Wheat`~ Precip_May_Aug_mm, data = wheat)c(slope =coef(rain_only)[2], R2 =summary(rain_only)$r.squared)## (b) adding the hot daysboth <-lm(`Spring Wheat`~ Precip_May_Aug_mm + Days_Above_30, data = wheat)c(coef(both)[2:3], R2 =summary(both)$r.squared)## how the two predictors move togethercor(wheat$Precip_May_Aug_mm, wheat$Days_Above_30)
> library(tidyverse)
>
> wheat <- read_csv("data/rm_canola_weather.csv", show_col_types = FALSE) |>
+ filter(!is.na(`Spring Wheat`))
>
> ## (a) rainfall alone
> rain_only <- lm(`Spring Wheat` ~ Precip_May_Aug_mm, data = wheat)
> c(slope = coef(rain_only)[2], R2 = summary(rain_only)$r.squared)
slope.Precip_May_Aug_mm R2
0.06250651 0.11659422
>
> ## (b) adding the hot days
> both <- lm(`Spring Wheat` ~ Precip_May_Aug_mm + Days_Above_30, data = wheat)
> c(coef(both)[2:3], R2 = summary(both)$r.squared)
Precip_May_Aug_mm Days_Above_30 R2
-0.001149179 -1.065533812 0.404113133
>
> ## how the two predictors move together
> cor(wheat$Precip_May_Aug_mm, wheat$Days_Above_30)
[1] -0.5441098
(a) 0.0625 bu/ac per mm, \(R^2\) = 0.117.
(b) Rainfall -0.0011, hot days -1.066, \(R^2\) = 0.404.
(c) The coefficient does not just shrink – it changes sign, from +0.0625 to -0.0011.
Rainfall and hot days are strongly negatively correlated: the correlation is -0.544, because wet seasons in Saskatchewan are cool seasons. With hot days left out, the rainfall coefficient absorbs their effect, so part of what looked like the benefit of rain was really the absence of heat. Once heat is measured separately almost nothing is left for rainfall to explain, and the little that remains comes out marginally negative.
This is omitted variable bias large enough to reverse a sign, not merely shrink one. The two-variable model is the one to quote – though -0.0011 should not be read as “rain slightly harms wheat” either, since it is small and other things are still omitted.
Question 2
Using RM yields and weather, one row per Saskatchewan RM per year from 2015 to 2024. Yields are in bu/ac, Precip_May_Aug_mm is total May-August rainfall in mm, and Days_Above_30 counts days hotter than 30 °C.
(a) Run two regressions, one with barley yields as the dependent variable and one with canola yields as the dependent variable. Include the following explanatory variables in the model: rainfall, days above 30 °C, latitude and soil_zone.
(b) Which of the crops is more heat sensitive, according to your results? Compare the hot-day coefficients in bushels, and then as a share of each crop’s average yield.
(c) Is one crop better suited to a certain growing region than another? Explain.
Answer
library(tidyverse)rm_data <-read_csv("data/rm_canola_weather.csv", show_col_types =FALSE) |>mutate(soil_zone =factor(soil_zone,levels =c("Brown", "Dark Brown", "Grey/Black")))## one model per crop, each on the rows where that crop has a yieldbarley <- rm_data |>filter(!is.na(Barley))canola <- rm_data |>filter(!is.na(Canola))barley_fit <-lm(Barley ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone, data = barley)canola_fit <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone, data = canola)round(coef(barley_fit), 4)round(coef(canola_fit), 4)## average yields, for the comparison in part (b)c(barley_mean =round(mean(barley$Barley), 1),canola_mean =round(mean(canola$Canola), 1))
> library(tidyverse)
>
> rm_data <- read_csv("data/rm_canola_weather.csv", show_col_types = FALSE) |>
+ mutate(soil_zone = factor(soil_zone,
+ levels = c("Brown", "Dark Brown", "Grey/Black")))
>
> ## one model per crop, each on the rows where that crop has a yield
> barley <- rm_data |> filter(!is.na(Barley))
> canola <- rm_data |> filter(!is.na(Canola))
>
> barley_fit <- lm(Barley ~ Precip_May_Aug_mm + Days_Above_30 + latitude +
+ soil_zone, data = barley)
> canola_fit <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude +
+ soil_zone, data = canola)
>
> round(coef(barley_fit), 4)
(Intercept) Precip_May_Aug_mm Days_Above_30 latitude
108.0865 0.0115 -1.2523 -0.8721
soil_zoneDark Brown soil_zoneGrey/Black
9.7146 14.0241
> round(coef(canola_fit), 4)
(Intercept) Precip_May_Aug_mm Days_Above_30 latitude
26.4649 0.0077 -0.7296 0.2612
soil_zoneDark Brown soil_zoneGrey/Black
2.7251 4.0360
>
> ## average yields, for the comparison in part (b)
> c(barley_mean = round(mean(barley$Barley), 1),
+ canola_mean = round(mean(canola$Canola), 1))
barley_mean canola_mean
59.3 35.0
(a) Barley: rainfall 0.0115, hot days -1.2523, latitude -0.8721, Dark Brown +9.715, Grey/Black +14.024, adjusted \(R^2\) 0.508. Canola: 0.0077, -0.7296, +0.2612, +2.725, +4.036, adjusted \(R^2\) 0.530.
(b) In bushels, barley – it loses 1.25 bu/ac per hot day against canola’s 0.73. But that comparison is unfair, because barley averages 59.3 bu/ac and canola only 35.0, so a bushel is a smaller share of the crop for barley. As a proportion of average yield a hot day costs barley \(1.25/59.3\) and canola \(0.73/35.0\), both about 2.1 per cent – essentially identical. The honest answer is that these two crops are about equally heat-sensitive once you account for how much they yield.
(c) Yes, and the soil zone and latitude coefficients show it. Barley gains far more from moving into the better soil zones (+9.7 and +14.0 bu/ac against canola’s +2.7 and +4.0), but its latitude coefficient is negative (-0.87) while canola’s is positive (+0.26).
So the two crops respond differently to where they are grown, and on these estimates barley gains more from the better soil zones while canola holds up better further north.
Be careful how far that is taken. Both crops are grown in essentially the same places here – 294 RMs for barley against 295 for canola, over the same range of latitudes – so this is not a story about different growing areas. What the coefficients describe is how each crop’s yield varies with zone and latitude, holding the weather variables fixed. Turning that into advice about which crop to plant would need the prices of both, the rest of the rotation, and some reason to think these associations would hold for a crop moved somewhere it is not currently grown.”
Question 3
Use Canola trial, which reports yield, rainfall and nitrogen use (fertilizer_kg_ha) across 117 fields.
(a) Run a regression with yield as the response variable and nitrogen rate, rainfall and variety as the explanatory variables. Report the nitrogen coefficient and interpret it in a sentence with units.
(b) Fit the same model with nitrogen alone. How much does the nitrogen coefficient move?
(c) Does the omission of rainfall bias the nitrogen coefficient? Explain why or why not.
Answer
library(tidyverse)trial <-read_csv("data/canola_trial_complete.csv", show_col_types =FALSE)## (a) the full modelfull <-lm(yield_bu_acre ~ fertilizer_kg_ha + rainfall_mm + variety, data = trial)round(coef(full), 4)## (b) nitrogen on its ownnitrogen_only <-lm(yield_bu_acre ~ fertilizer_kg_ha, data = trial)round(c(full =coef(full)[2], nitrogen_only =coef(nitrogen_only)[2]), 4)## (c) do the two predictors move together?cor(trial$fertilizer_kg_ha, trial$rainfall_mm)
> library(tidyverse)
>
> trial <- read_csv("data/canola_trial_complete.csv", show_col_types = FALSE)
>
> ## (a) the full model
> full <- lm(yield_bu_acre ~ fertilizer_kg_ha + rainfall_mm + variety, data = trial)
> round(coef(full), 4)
(Intercept) fertilizer_kg_ha rainfall_mm varietyDEKALB
23.3141 0.1345 0.0548 0.7230
varietyInVigor
1.5414
>
> ## (b) nitrogen on its own
> nitrogen_only <- lm(yield_bu_acre ~ fertilizer_kg_ha, data = trial)
>
> round(c(full = coef(full)[2], nitrogen_only = coef(nitrogen_only)[2]), 4)
full.fertilizer_kg_ha nitrogen_only.fertilizer_kg_ha
0.1345 0.1255
>
> ## (c) do the two predictors move together?
> cor(trial$fertilizer_kg_ha, trial$rainfall_mm)
[1] -0.0902854
(a) Comparing two fields with the same rainfall growing the same variety, one more kg/ha of nitrogen goes with about 0.134 bu/ac more yield.
(b) From 0.1345 with the controls to 0.1255 without them, a movement of 0.009 – about 7% of the coefficient.
(c) Very little movement, and the reason is the weak correlation rather than the small change itself.
A small movement means the product\(\hat{\beta}_2\tilde{\delta}_1\) is small, and that can happen three ways: the correlation is weak, the omitted variable matters little to yield, or two omitted variables pull in opposite directions. Rainfall and variety were dropped together here, so the movement alone cannot separate them.
The correlation is -0.09, a weak linear relationship, which is what you would expect from a designed trial: nitrogen rates were assigned rather than chosen by farmers, so there is no reason for the wetter fields to have received more fertilizer. The weaker that correlation, the smaller the bias from leaving rainfall out – the second of the two conditions from Defining omitted variable bias.
Note this is a statement about size, not a clean zero. Holding variety in the model and dropping only rainfall moves nitrogen from 0.1345 to 0.1226, about 9% of the coefficient. Whether 9% matters depends on what the number is for.
Question 4
Using RM yields and weather. The column Population_2021 is the number of people living in the RM at the 2021 census.
(a) Run a regression with barley yield as the response variable and rainfall, days above 30 °C, latitude and soil_zone as the explanatory variables. Then add Population_2021 to it.
(b) Report the coefficient on population and its p-value. Is it statistically significant?
(c) Should population be in a model of barley yield? Give your reasoning.
(a) Both models are in the output above: barley yield on the weather and place variables, then the same model with Population_2021 added. The \(R^2\) rises from 0.5084 to 0.5093.
(b) The coefficient is 0.00086 bu/ac per person with a p-value of 0.024, so it is statistically significant at the 5% level, though not overwhelmingly so. An RM with 1,000 more people is predicted to yield about 0.9 bu/ac more barley.
(c) Generally no. There is no mechanism by which the number of people living in an RM changes what a barley crop yields. The coefficient is picking up something else that varies with population – RMs with towns differ systematically in soil, in how much land is farmed, and in proximity to services – so it cannot be read as an effect of population. While it may improve in-sample prediction, it is unlikely to improve out of sample prediction.
Question 5
Using RM yields and weather, one row per Saskatchewan RM per year from 2015 to 2024. Yields are in bu/ac, Precip_May_Aug_mm is total May-August rainfall in mm, and Days_Above_30 counts days hotter than 30 °C.
(a) Keep the rows with an oat yield. Run a regression of oat yield on rainfall, days above 30 °C and Year entered as a number. Report the coefficient on Year and say what it claims.
(b) Run the same regression with factor(Year) instead, and compare the two \(R^2\) values.
(c) Compute the mean oat yield in each year, and use it to explain why the two models differ. Which form would you use?
Answer
library(tidyverse)oats <-read_csv("data/rm_canola_weather.csv", show_col_types =FALSE) |>filter(!is.na(Oats))## (a) year as a number -- one straight-line trendnumeric_year <-lm(Oats ~ Precip_May_Aug_mm + Days_Above_30 + Year, data = oats)c(Year =coef(numeric_year)["Year"], R2 =summary(numeric_year)$r.squared)## (b) year as a category -- each year gets its own coefficientfactor_year <-lm(Oats ~ Precip_May_Aug_mm + Days_Above_30 +factor(Year),data = oats)summary(factor_year)$r.squared## (c) what actually happened, year by yearoats |>group_by(Year) |>summarise(mean_yield =round(mean(Oats), 1))
> library(tidyverse)
>
> oats <- read_csv("data/rm_canola_weather.csv", show_col_types = FALSE) |>
+ filter(!is.na(Oats))
>
> ## (a) year as a number -- one straight-line trend
> numeric_year <- lm(Oats ~ Precip_May_Aug_mm + Days_Above_30 + Year, data = oats)
> c(Year = coef(numeric_year)["Year"], R2 = summary(numeric_year)$r.squared)
Year.Year R2
1.924589 0.466151
>
> ## (b) year as a category -- each year gets its own coefficient
> factor_year <- lm(Oats ~ Precip_May_Aug_mm + Days_Above_30 + factor(Year),
+ data = oats)
> summary(factor_year)$r.squared
[1] 0.4918503
>
> ## (c) what actually happened, year by year
> oats |>
+ group_by(Year) |>
+ summarise(mean_yield = round(mean(Oats), 1))
# A tibble: 10 × 2
Year mean_yield
<dbl> <dbl>
1 2015 71.7
2 2016 84.3
3 2017 79.8
4 2018 77.9
5 2019 89.4
6 2020 89.1
7 2021 45.1
8 2022 85.8
9 2023 71.1
10 2024 84
(a) 1.92, with \(R^2\) = 0.466. Entered as a number, the coefficient claims oat yields rose by a constant 1.92 bu/ac in every year from 2015 to 2024 – a single straight line through the decade.
(b)\(R^2\) rises to 0.492 with year as a factor.
(c) The yearly means are nothing like a straight line: between about 71 and 89 bu/ac in most years, and 45.1 in 2021, the drought year. The numeric term can only tilt a line, so it forces one constant step per year onto a decade that did not move that way. The factor version gives each year its own coefficient, which is where the extra fit comes from.
Note what the 1.92 is and is not. It is the trend after adjusting for rainfall and heat – 2021 was hot and dry, so the model attributes that year’s low yield to the weather and the leftover year effect comes out positive. Regressed on year alone the slope is -0.22, in the opposite direction. The two numbers answer different questions and neither is wrong.
Use the factor version to describe what happened in these ten years. The numeric version is for when you want a single trend estimate, or need a prediction for a year outside the data, which a factor cannot give you.
Section 5 — Measuring omitted variable bias in R
These questions ask you to fit a model without a control variable and then with it, run the auxiliary regression, and verify Equation 27.1 numerically. Several use data from outside agriculture, because the algebra is the same everywhere and the classic examples are clearer.
Question 1
Using Wages and IQ, 935 working men. lwage is the log of monthly earnings, educ is years of schooling, IQ is a test score.
(a) Run a regression of log monthly earnings (lwage) on years of schooling (educ), and report the coefficient.
(b) Now add IQ to the model as an additional explanatory variable and report both coefficients.
(c) Fit the auxiliary regression of the omitted variable on the included one, IQ ~ educ, and report the coefficient. Does this coefficient make sense to you?
(d) Explain why, and by how much, the education-only coefficient is biased by the exclusion of IQ.
Answer
library(tidyverse)wages <-read_csv("data/wage_iq.csv", show_col_types =FALSE)## (a) education onlyeduc_only <-lm(lwage ~ educ, data = wages)## (b) adding IQwith_iq <-lm(lwage ~ educ + IQ, data = wages)## (c) the auxiliary regression: the omitted variable on the included oneaux <-lm(IQ ~ educ, data = wages)b_educ_only <-coef(educ_only)["educ"]b1 <-coef(with_iq)["educ"]b2 <-coef(with_iq)["IQ"]d1 <-coef(aux)["educ"]c(educ_only = b_educ_only, b1 = b1, b2 = b2, delta1 = d1)## (d) the identityc(bias = b2 * d1, reconstructed = b1 + b2 * d1, actual = b_educ_only)
> library(tidyverse)
>
> wages <- read_csv("data/wage_iq.csv", show_col_types = FALSE)
>
> ## (a) education only
> educ_only <- lm(lwage ~ educ, data = wages)
>
> ## (b) adding IQ
> with_iq <- lm(lwage ~ educ + IQ, data = wages)
>
> ## (c) the auxiliary regression: the omitted variable on the included one
> aux <- lm(IQ ~ educ, data = wages)
>
> b_educ_only <- coef(educ_only)["educ"]
> b1 <- coef(with_iq)["educ"]
> b2 <- coef(with_iq)["IQ"]
> d1 <- coef(aux)["educ"]
>
> c(educ_only = b_educ_only, b1 = b1, b2 = b2, delta1 = d1)
educ_only.educ b1.educ b2.IQ delta1.educ
0.059839204 0.039119899 0.005863131 3.533829187
>
> ## (d) the identity
> c(bias = b2 * d1, reconstructed = b1 + b2 * d1, actual = b_educ_only)
bias.IQ reconstructed.educ actual.educ
0.0207193 0.0598392 0.0598392
(a) 0.05984. Since the response is logged, a year more schooling goes with about 6% higher monthly earnings.
(b)educ drops to 0.03912; IQ is 0.005863.
(c) 3.53383 – each additional year of schooling goes with about 3.5 more IQ points.
Yes, that makes sense, in both directions. People who test well tend to stay in school longer, and time spent in school plausibly raises a measured test score as well. More to the point for this question, a coefficient this far from zero says the two predictors are substantially related, which is exactly the condition under which omitting one distorts the other.
(d)\(0.03912 + 0.005863 \times 3.53383 = 0.05984\), reproducing the education-only coefficient exactly. The bias is \(0.005863 \times 3.53383 = 0.02072\).
The sign is positive because both pieces are positive: IQ raises earnings (\(\hat{\beta}_2 > 0\)) and IQ rises with schooling (\(\tilde{\delta}_1 > 0\)). Leaving IQ out therefore makes schooling look better paid than it is – the education coefficient absorbs part of what IQ was doing.
The size is \(0.02072 / 0.05984 = 34.6\%\)of the education-only estimate – against the IQ-adjusted estimate it is 53.0%, so say which denominator you mean.
The identity is algebra, and it does not establish that 0.0391 is the true return to schooling. Family background, school quality and access to credit are still missing.
Question 2
In the data Job training there are 1,130 people observed. train is 1 if the person took a training programme in 1997, earn98 is 1998 earnings in $000s, and earn96 is 1996 earnings.
(a) Fit earn98 ~ train. What does the coefficient appear to say about the training programme?
(b) Add earn96 and refit. What happens to the coefficient on training?
(c) Explain why the coefficient changed between part a and part b, using the formula we have learned for omitted variable bias.
Answer
library(tidyverse)jt <-read_csv("data/job_training.csv", show_col_types =FALSE)## (a) training onlytrain_only <-lm(earn98 ~ train, data = jt)## (b) controlling for what they earned before the programmewith_prior <-lm(earn98 ~ train + earn96, data = jt)## (c) the auxiliary regressionaux <-lm(earn96 ~ train, data = jt)round(c(train_only =coef(train_only)["train"],with_prior =coef(with_prior)["train"],b2 =coef(with_prior)["earn96"], delta1 =coef(aux)["train"]), 4)## mean prior earnings in each groupjt |>group_by(train) |>summarise(mean_earn96 =round(mean(earn96), 2))
> library(tidyverse)
>
> jt <- read_csv("data/job_training.csv", show_col_types = FALSE)
>
> ## (a) training only
> train_only <- lm(earn98 ~ train, data = jt)
>
> ## (b) controlling for what they earned before the programme
> with_prior <- lm(earn98 ~ train + earn96, data = jt)
>
> ## (c) the auxiliary regression
> aux <- lm(earn96 ~ train, data = jt)
>
> round(c(train_only = coef(train_only)["train"],
+ with_prior = coef(with_prior)["train"],
+ b2 = coef(with_prior)["earn96"], delta1 = coef(aux)["train"]), 4)
train_only.train with_prior.train b2.earn96 delta1.train
-2.0501 2.5120 0.4212 -10.8314
>
> ## mean prior earnings in each group
> jt |>
+ group_by(train) |>
+ summarise(mean_earn96 = round(mean(earn96), 2))
# A tibble: 2 × 2
train mean_earn96
<dbl> <dbl>
1 0 15.8
2 1 4.97
(a) The coefficient is -2.05, so trainees earned about $2,050 less in 1998 than non-trainees. Taken at face value the programme looks harmful.
(b) It flips sign, to +2.51. Controlling for what people earned before the programme, trainees earned about $2,510 more.
(c) The two terms of the bias formula are \(\hat{\beta}_2 = 0.4212\) (prior earnings strongly predict later earnings) and \(\tilde{\delta}_1 = -10.831\) (trainees earned $10,831 less in 1996 than non-trainees). The bias is their product: \(0.4212 \times (-10.831) = -4.562\), and \(2.512 + (-4.562) = -2.050\), which is exactly the training-only coefficient.
The bias is negative because a positive times a negative is negative: prior earnings raise later earnings, but are strongly negatively related to taking the training. In plain terms, trainees averaged $4,970 in 1996 against $15,800 for non-trainees – the programme went to people who were already earning far less. The raw comparison is mostly measuring who was selected into the programme rather than what it did, and the bias is large enough to flip the sign.
Question 3
Using Spray and yield, 500 simulated fields. spray_passes is the number of fungicide passes, pest_index is a 0-10 measure of pest pressure, yield_bu_ac is the yield.
(a) Run a regression with yield as the response variable and spray passes as the only explanatory variable. Report the coefficient.
(b) Add pest_index as an additional explanatory variable. Report both coefficients.
(c) Explain why the coefficient on spray passes changed from part a to part b, using the formula we have learned for omitted variable bias.
Answer
library(tidyverse)spray <-read_csv("data/spray_yield.csv", show_col_types =FALSE)## (a) spray passes onlyspray_only <-lm(yield_bu_ac ~ spray_passes, data = spray)## (b) adding pest pressurewith_pest <-lm(yield_bu_ac ~ spray_passes + pest_index, data = spray)## (c) the auxiliary regressionaux <-lm(pest_index ~ spray_passes, data = spray)round(c(spray_only =coef(spray_only)["spray_passes"],with_pest =coef(with_pest)["spray_passes"],b2 =coef(with_pest)["pest_index"],delta1 =coef(aux)["spray_passes"]), 4)
(c) The two terms are \(\hat{\beta}_2 = -1.6688\) (pests cut yield) and \(\tilde{\delta}_1 = 0.5764\) (spraying happens where the pests are). The bias is their product:
\[-1.6688 \times 0.5764 = -0.962\]
and \(2.871 + (-0.962) = 1.909\), which is the spray-only coefficient.
The bias is negative – a negative times a positive – so leaving pest pressure out makes spraying look less effective than it is, not more. The fields that got sprayed most were the ones with the worst pest problems, and their lower yields get charged to the spraying.
This data was simulated with a true spray effect of exactly 3.0 bu/ac per pass. The model including pest pressure recovers 2.871, close to the truth; the spray-only model gives 1.909, understating the benefit by more than a third.
Question 4
Using Housing and air pollution, 506 census tracts in the Boston area in 1970. price is the median house price in the tract, nox is the nitrogen oxide concentration in the air in parts per 100 million, and rooms is the average number of rooms per house. The file also holds lprice and lnox, the logs of price and pollution, which are the forms this question uses.
(a) Run a regression of lprice on lnox, and report the coefficient. Both variables are logged, so read it as the percentage change in price from a one percent increase in pollution.
(b) Add rooms as an additional explanatory variable and refit. Report the pollution coefficient and the rooms coefficient.
(c) Run the auxiliary regression of rooms on lnox. Report the coefficient and say what it tells you about where the larger houses are.
(d) Verify Equation 27.1 numerically, and explain the sign of the bias – why does leaving house size out make pollution look worse for prices than it is?
Answer
library(tidyverse)housing <-read_csv("data/housing_pollution.csv", show_col_types =FALSE)## (a) pollution alonenox_only <-lm(lprice ~ lnox, data = housing)## (b) adding house sizewith_rooms <-lm(lprice ~ lnox + rooms, data = housing)## (c) the auxiliary regressionaux <-lm(rooms ~ lnox, data = housing)round(c(nox_only =coef(nox_only)["lnox"],with_rooms =coef(with_rooms)["lnox"],b2 =coef(with_rooms)["rooms"],delta1 =coef(aux)["lnox"]), 4)## (d) the identityb2 <-coef(with_rooms)["rooms"]d1 <-coef(aux)["lnox"]round(c(bias = b2 * d1,reconstructed =coef(with_rooms)["lnox"] + b2 * d1,actual =coef(nox_only)["lnox"]), 4)
(a) -1.0431. A one percent increase in nitrogen oxide goes with about a one percent fall in house prices.
(b) The pollution coefficient becomes less negative, moving from -1.0431 to -0.7177. rooms is 0.3059: since the response is logged, a tract with one more average room has prices higher by a factor of \(e^{0.3059}\), about 36%. (The usual rule of thumb, reading 0.3059 as 31%, is noticeably off at a coefficient this large.)
(c) -1.0639. Tracts with dirtier air have smaller houses. That is the key fact for this question: the two explanatory variables are related, so the condition for bias is met.
(d) The bias is the product:
\[0.3059 \times (-1.0639) = -0.3255\]
and \(-0.7177 + (-0.3255) = -1.0431\), which is the pollution-only coefficient.
The bias is negative – a positive times a negative. Big houses sell for more (\(\hat{\beta}_2 > 0\)), and big houses are in the cleaner tracts (\(\tilde{\delta}_1 < 0\)). With house size left out, the pollution coefficient has to stand in for both things at once: part of what it reports as the cost of dirty air is really the cost of a smaller house. Pollution ends up looking about 45% more damaging to prices than it does once size is held constant.
The adjusted estimate of -0.72 is a better comparison than -1.04, but it is still not the causal effect of air quality. Tracts also differ in distance to employment, crime and property taxes, all of which are in this file and none of which are in the model.
Question 5
Using College and wages, 6,763 people surveyed in the 1990s. lwage is the log of the hourly wage, jc is years of schooling completed at a junior (two-year) college, and univ is years completed at a four-year university. Of the 6,763, 824 attended a junior college only, 2,627 a university only, 829 both, and 2,483 neither.
(a) Run a regression of lwage on jc alone, and report the coefficient.
(b) Add univ and refit. Report both coefficients. What happens to the junior-college coefficient?
(c) Run the auxiliary regression of univ on jc. The coefficient is negative. Explain why that is what you would expect from these two variables.
(d) Verify Equation 27.1 numerically, and explain what a student who ran only the regression in part a would wrongly conclude about going to a junior college.
Answer
library(tidyverse)college <-read_csv("data/college_wages.csv", show_col_types =FALSE)## (a) junior college alonejc_only <-lm(lwage ~ jc, data = college)## (b) adding universitywith_univ <-lm(lwage ~ jc + univ, data = college)## (c) the auxiliary regressionaux <-lm(univ ~ jc, data = college)round(c(jc_only =coef(jc_only)["jc"],with_univ =coef(with_univ)["jc"],b2 =coef(with_univ)["univ"],delta1 =coef(aux)["jc"]), 5)## (d) the identityb2 <-coef(with_univ)["univ"]d1 <-coef(aux)["jc"]round(c(bias = b2 * d1,reconstructed =coef(with_univ)["jc"] + b2 * d1,actual =coef(jc_only)["jc"]), 5)
(a) 0.04450. A year at a junior college goes with about 4.5% higher wages.
(b)jc rises to 0.07012 and univ is 0.06903. The junior-college coefficient gets larger, by more than half.
(c) -0.37116. Someone with more junior-college years has, on average, fewer university years – the two substitute for each other in this sample. That is the opposite of the usual textbook setup, where the omitted variable rises with the included one. Whether one causes the other is not something the auxiliary regression can tell you; what matters here is that they move together, which is the condition for bias.
(d) The bias is:
\[0.06903 \times (-0.37116) = -0.02562\]
and \(0.07012 + (-0.02562) = 0.04450\), the junior-college-only coefficient.
Because \(\tilde{\delta}_1\) is negative while \(\hat{\beta}_2\) is positive, the bias is negative and the junior-college-only estimate is too low. It is understated by 58% of its own value.
The wrong conclusion is that junior college is associated with about a third less earnings per year than it is. The part-a regression compares junior-college years against the university years they substitute for rather than against no schooling at all. Once both are in the model the two coefficients are almost identical – 0.0701 and 0.0690 – so in this sample a year at either is associated with about the same wage. That is an adjusted association, not a return: selection into schooling is still uncontrolled.
The direction of the bias needs both signs. Here the omitted variable raises wages (\(\hat{\beta}_2 > 0\)) and falls as the included one rises (\(\tilde{\delta}_1 < 0\)), so the product is negative and the estimate is pulled down. Flip either sign and the direction flips with it.
Section 6 — Measuring prediction accuracy in R
These questions ask you to compute fitted values, residuals and the three accuracy measures, and to compare performance in sample against a hold-out sample.
Question 1
Using RM yields and weather, run a regression with canola yield as the response variable and rainfall, days above 30 °C, latitude and soil_zone as the explanatory variables.
(a) Compute the \(R^2\) from the sum of squared errors and the total sum of squares, rather than reading it off summary(). Check that the two agree.
(b) Compute the residuals and use them to compute the MAE and RMSE.
(c) If you were reporting the predictive accuracy of this model to a group of producers, which of these measures (\(R^2\), MAE, RMSE) would you report, and why?
Answer
library(tidyverse)canola <-read_csv("data/rm_canola_weather.csv", show_col_types =FALSE) |>filter(!is.na(Canola)) |>mutate(soil_zone =factor(soil_zone,levels =c("Brown", "Dark Brown", "Grey/Black")))model4 <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,data = canola)## (a) R-squared from its partssse <-sum(resid(model4)^2)sst <-sum((canola$Canola -mean(canola$Canola))^2)c(SSE = sse, SST = sst, R2 =1- sse / sst,check =summary(model4)$r.squared)## (b) the two error measures, in bushelsc(MAE =mean(abs(resid(model4))),RMSE =sqrt(mean(resid(model4)^2)),sd_y =sd(canola$Canola))
(a) SSE = 150,138.2 and SST = 319,972.5, so \(R^2 = 1 - 150138.2/319972.5 = 0.5308\) – matching summary() exactly.
(b) MAE 5.698 bu/ac and RMSE 7.190 bu/ac.
(c)MAE or RMSE, not \(R^2\) – because they are in bushels per acre, the units the audience actually thinks in. “This model is typically off by about 5.7 bu/ac” means something to a producer; “the model explains 53% of the variation” does not tell them whether to use it.
Between the two, MAE is easier to explain – it is the average size of a miss. RMSE weights large misses more heavily, and the gap between 5.70 and 7.19 says some RM-years are badly missed.
Give them a benchmark too. Guessing the average every time gives an RMSE of 10.50 bu/ac (the standard deviation of yield) and an MAE of 8.36. Against those, the model’s 7.19 and 5.70 cut a typical error by about a third. Say that these are in-sample figures: the model is being judged on the same RM-years it was fitted to.
Question 2
For this question, the RM canola data has been split into a training set, canola_yield_training.csv (2,032 RM-years), and a test set, canola_yield_test.csv (872 RM-years). The question will ask you to run a regression on the training set and evaluate it on the test set.
(a) Run a regression of canola yield on rainfall, days above 30 °C, latitude and soil_zone. Report the in-sample and out-of-sample RMSE.
(b) Suppose we add Population_2021 to the model. Would you expect that to improve the in-sample RMSE, the out-of-sample RMSE, both, or neither? Explain before you run it.
(c) Run the same regression with Population_2021 added. Did the two RMSEs change the way you expected?
Answer
library(tidyverse)train <-read_csv("data/canola_yield_training.csv", show_col_types =FALSE)test <-read_csv("data/canola_yield_test.csv", show_col_types =FALSE)## (a) the weather-and-place model, fitted on the training rowsbase <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,data = train)## in-sample errors are the residuals of the fitted modelbase_in <-sqrt(mean(resid(base)^2))## out-of-sample errors: predict the test rows, subtract the actual yieldsbase_out_err <- test$Canola -predict(base, test)base_out <-sqrt(mean(base_out_err^2))round(c(in_sample = base_in, out_of_sample = base_out), 4)## (c) the same model with population addedwith_pop <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone + Population_2021, data = train)pop_in <-sqrt(mean(resid(with_pop)^2))pop_out_err <- test$Canola -predict(with_pop, test)pop_out <-sqrt(mean(pop_out_err^2))round(c(in_sample = pop_in, out_of_sample = pop_out), 4)
> library(tidyverse)
>
> train <- read_csv("data/canola_yield_training.csv", show_col_types = FALSE)
> test <- read_csv("data/canola_yield_test.csv", show_col_types = FALSE)
>
> ## (a) the weather-and-place model, fitted on the training rows
> base <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,
+ data = train)
>
> ## in-sample errors are the residuals of the fitted model
> base_in <- sqrt(mean(resid(base)^2))
>
> ## out-of-sample errors: predict the test rows, subtract the actual yields
> base_out_err <- test$Canola - predict(base, test)
> base_out <- sqrt(mean(base_out_err^2))
>
> round(c(in_sample = base_in, out_of_sample = base_out), 4)
in_sample out_of_sample
7.1769 7.2285
>
> ## (c) the same model with population added
> with_pop <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude +
+ soil_zone + Population_2021, data = train)
>
> pop_in <- sqrt(mean(resid(with_pop)^2))
> pop_out_err <- test$Canola - predict(with_pop, test)
> pop_out <- sqrt(mean(pop_out_err^2))
>
> round(c(in_sample = pop_in, out_of_sample = pop_out), 4)
in_sample out_of_sample
7.1738 7.2165
(a) In-sample RMSE 7.177, out-of-sample RMSE 7.229. The two are almost identical: four predictors fitted on 2,032 rows leaves little scope to chase noise.
(b) The in-sample RMSE must improve, or at worst stay the same: least squares could set the population coefficient to zero and recover the old fit. Expect the improvement to be tiny.
The out-of-sample RMSE has no such guarantee and could move either way. Population has no agronomic mechanism, so it is close to a coin flip – a small gain if it proxies for something real like land type, a small loss if it is noise the model has fitted.
(c) In-sample RMSE falls from 7.1769 to 7.1738, as it had to. Out-of-sample RMSE also falls slightly, from 7.2285 to 7.2165.
So population helped a little on both – but note how little. The out-of-sample gain is 0.012 bu/ac against a typical error of 7.2, and the p-value is 0.19. This variable neither helps nor hurts in any way you could act on, so leave it out: it adds nothing and invites a reader to interpret a coefficient with no mechanism behind it.
(a) Run two regressions with spring wheat yield as the response variable. One with rainfall as the only explanatory variable, and one with rainfall, days above 30 °C, latitude and soil_zone.
(b) Which model would you expect to have the higher in-sample \(R^2\)? Which would you expect to have the higher out-of-sample RMSE? Explain before you look.
(c) Report both the in-sample \(R^2\) and the out-of-sample RMSE for both models.
Answer
library(tidyverse)train <-read_csv("data/rm_weather_training.csv", show_col_types =FALSE)test <-read_csv("data/rm_weather_test.csv", show_col_types =FALSE)## (a) one predictor, then fourrain_only <-lm(`Spring Wheat`~ Precip_May_Aug_mm, data = train)full <-lm(`Spring Wheat`~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone, data = train)## (c) in-sample R-squared, read off each modelround(c(rain_only =summary(rain_only)$r.squared,full =summary(full)$r.squared), 4)## out-of-sample RMSE: predict the test rows, subtract the actual yieldsrain_err <- test$`Spring Wheat`-predict(rain_only, test)full_err <- test$`Spring Wheat`-predict(full, test)round(c(rain_only =sqrt(mean(rain_err^2)),full =sqrt(mean(full_err^2))), 4)
> library(tidyverse)
>
> train <- read_csv("data/rm_weather_training.csv", show_col_types = FALSE)
> test <- read_csv("data/rm_weather_test.csv", show_col_types = FALSE)
>
> ## (a) one predictor, then four
> rain_only <- lm(`Spring Wheat` ~ Precip_May_Aug_mm, data = train)
> full <- lm(`Spring Wheat` ~ Precip_May_Aug_mm + Days_Above_30 +
+ latitude + soil_zone, data = train)
>
> ## (c) in-sample R-squared, read off each model
> round(c(rain_only = summary(rain_only)$r.squared,
+ full = summary(full)$r.squared), 4)
rain_only full
0.1224 0.5065
>
> ## out-of-sample RMSE: predict the test rows, subtract the actual yields
> rain_err <- test$`Spring Wheat` - predict(rain_only, test)
> full_err <- test$`Spring Wheat` - predict(full, test)
>
> round(c(rain_only = sqrt(mean(rain_err^2)),
+ full = sqrt(mean(full_err^2))), 4)
rain_only full
12.7483 9.5582
(a) Both regressions are fitted on the training rows in the output above – spring wheat on rainfall alone, then on rainfall, days above 30 °C, latitude and soil_zone.
(b) The full model must have an equal or higher in-sample \(R^2\). That is guaranteed by construction: adding predictors cannot reduce \(R^2\) on the rows the model was fitted on.
The out-of-sample RMSE is a genuine question rather than a guarantee, and here the full model should win it too. Heat, latitude and soil zone all have obvious mechanisms for affecting wheat yield, and real predictors usually help out of sample as well as in. You would expect the opposite from variables with no mechanism, or from many of them in a small sample.
Both expectations hold, and the out-of-sample gap is large – the full model’s typical miss is 3.2 bu/ac smaller on data neither model has seen. This is the ordinary case, and it is worth seeing alongside the overfitting examples: adding predictors that genuinely belong improves prediction on new data, not just the fit on old data. The two rankings agreeing is what you hope for; when they disagree, the out-of-sample figure is the one to act on.
Question 4
The data Indian Head yields and weather 36 years of crop yields for RM 156 (Indian Head) one row per year, 1990 to 2025. It contains the following weather data: Precip_May_Aug_mm is total May-August rainfall, GDD_base5 is growing degree days, Days_Above_30 counts days hotter than 30 °C. It also contains columns for the yields of various crops
(a) Split the data in time: the 30 years from 1990 to 2019 are the training set, and the 6 years from 2020 to 2025 are the test set.
(b) Using the training set only, run a regression with canola yields as the response variable and Year, rainfall, growing degree days and days above 30 °C as explanatory variables. Report the in-sample \(R^2\), MAE and RMSE.
(c) Report the MAE and RMSE for the same model on the test years.
(d) Does the model predict better in sample or out of sample? Explain why that is the expected direction, and say whether the gap here looks like a model that has learned something real.
Answer
library(tidyverse)indian_head <-read_csv("data/indian_head_yields_weather.csv",show_col_types =FALSE)## (a) split on the year, not at randomtrain <- indian_head |>filter(Year <=2019)test <- indian_head |>filter(Year >=2020)c(train_years =nrow(train), test_years =nrow(test))## (b) fit on the training years onlymodel <-lm(Canola ~ Year + Precip_May_Aug_mm + GDD_base5 + Days_Above_30,data = train)## in-sample: the residuals of the fitted modelin_err <-resid(model)round(c(R2 =summary(model)$r.squared,MAE =mean(abs(in_err)),RMSE =sqrt(mean(in_err^2))), 4)## (c) out-of-sample: predict the test years, subtract the actual yieldsout_err <- test$Canola -predict(model, test)round(c(R2 =1-sum(out_err^2) /sum((test$Canola -mean(test$Canola))^2),MAE =mean(abs(out_err)),RMSE =sqrt(mean(out_err^2))), 4)
> library(tidyverse)
>
> indian_head <- read_csv("data/indian_head_yields_weather.csv",
+ show_col_types = FALSE)
>
> ## (a) split on the year, not at random
> train <- indian_head |> filter(Year <= 2019)
> test <- indian_head |> filter(Year >= 2020)
>
> c(train_years = nrow(train), test_years = nrow(test))
train_years test_years
30 6
>
> ## (b) fit on the training years only
> model <- lm(Canola ~ Year + Precip_May_Aug_mm + GDD_base5 + Days_Above_30,
+ data = train)
>
> ## in-sample: the residuals of the fitted model
> in_err <- resid(model)
>
> round(c(R2 = summary(model)$r.squared,
+ MAE = mean(abs(in_err)),
+ RMSE = sqrt(mean(in_err^2))), 4)
R2 MAE RMSE
0.6650 4.3815 5.3893
>
> ## (c) out-of-sample: predict the test years, subtract the actual yields
> out_err <- test$Canola - predict(model, test)
>
> round(c(R2 = 1 - sum(out_err^2) / sum((test$Canola - mean(test$Canola))^2),
+ MAE = mean(abs(out_err)),
+ RMSE = sqrt(mean(out_err^2))), 4)
R2 MAE RMSE
0.3049 4.8378 5.7867
(a) 30 training years and 6 test years.
(b) In sample: \(R^2\) 0.665, MAE 4.38 bu/ac, RMSE 5.39 bu/ac.
(c) Out of sample: MAE 4.84 bu/ac, RMSE 5.79 bu/ac. (The out-of-sample \(R^2\) is 0.305, which is computed against the mean of the test years.)
(d)In sample, on every measure. That is the usual direction: the coefficients were chosen to make the errors on these rows as small as possible, so training performance flatters the model. Any other set of rows is predicted by coefficients that were not tuned to it. It is not a guarantee – an easy test sample can come out better – but it is what you should expect.
The size of the gap is what matters, and here it is small – the typical miss grows from 4.38 to 4.84 bu/ac, about a 10% deterioration. That is what a model that has learned something real looks like. A model that had mostly memorised the training years would show a much larger gap, and the overfitting examples in Predicting in sample and out of sample are far worse than this.
Question 5
The data Indian Head yields and weather 36 years of crop yields for RM 156 (Indian Head) one row per year, 1990 to 2025. It contains the following weather data: Precip_May_Aug_mm is total May-August rainfall, GDD_base5 is growing degree days, Days_Above_30 counts days hotter than 30 °C. It also contains columns for the yields of various crops
(a) Split the data in time, training on 1990 to 2019 and testing on 2020 to 2025.
(b) Run two regressions using spring wheat yield as the response variable. The first uses Year, rainfall and growing degree days as the explanatory variables. The second adds Days_Above_30. Report the in-sample MAE and RMSE for both.
(c) Report the out-of-sample MAE and RMSE for both models.
(d) Which model would you use, and what did adding the extra variable buy you? Explain why the in-sample figures on their own could not have told you that.
Answer
library(tidyverse)indian_head <-read_csv("data/indian_head_yields_weather.csv",show_col_types =FALSE)## (a) split on the year, not at randomtrain <- indian_head |>filter(Year <=2019)test <- indian_head |>filter(Year >=2020)## (b) without the heat variable, then with itwithout_heat <-lm(`Spring Wheat`~ Year + Precip_May_Aug_mm + GDD_base5,data = train)with_heat <-lm(`Spring Wheat`~ Year + Precip_May_Aug_mm + GDD_base5 + Days_Above_30, data = train)## in sample, the errors are the residualsround(c(MAE_without =mean(abs(resid(without_heat))),RMSE_without =sqrt(mean(resid(without_heat)^2)),MAE_with =mean(abs(resid(with_heat))),RMSE_with =sqrt(mean(resid(with_heat)^2))), 4)## (c) out of sample, predict the test years and subtract the actual yieldserr_without <- test$`Spring Wheat`-predict(without_heat, test)err_with <- test$`Spring Wheat`-predict(with_heat, test)round(c(MAE_without =mean(abs(err_without)),RMSE_without =sqrt(mean(err_without^2)),MAE_with =mean(abs(err_with)),RMSE_with =sqrt(mean(err_with^2))), 4)
> library(tidyverse)
>
> indian_head <- read_csv("data/indian_head_yields_weather.csv",
+ show_col_types = FALSE)
>
> ## (a) split on the year, not at random
> train <- indian_head |> filter(Year <= 2019)
> test <- indian_head |> filter(Year >= 2020)
>
> ## (b) without the heat variable, then with it
> without_heat <- lm(`Spring Wheat` ~ Year + Precip_May_Aug_mm + GDD_base5,
+ data = train)
> with_heat <- lm(`Spring Wheat` ~ Year + Precip_May_Aug_mm + GDD_base5 +
+ Days_Above_30, data = train)
>
> ## in sample, the errors are the residuals
> round(c(MAE_without = mean(abs(resid(without_heat))),
+ RMSE_without = sqrt(mean(resid(without_heat)^2)),
+ MAE_with = mean(abs(resid(with_heat))),
+ RMSE_with = sqrt(mean(resid(with_heat)^2))), 4)
MAE_without RMSE_without MAE_with RMSE_with
5.0983 6.2034 4.9909 6.0902
>
> ## (c) out of sample, predict the test years and subtract the actual yields
> err_without <- test$`Spring Wheat` - predict(without_heat, test)
> err_with <- test$`Spring Wheat` - predict(with_heat, test)
>
> round(c(MAE_without = mean(abs(err_without)),
+ RMSE_without = sqrt(mean(err_without^2)),
+ MAE_with = mean(abs(err_with)),
+ RMSE_with = sqrt(mean(err_with^2))), 4)
MAE_without RMSE_without MAE_with RMSE_with
8.2255 9.6502 7.3506 9.4155
(a) 30 training years (1990 to 2019) and 6 test years (2020 to 2025).
(b) Without the heat variable, MAE 5.098 and RMSE 6.203. With it, MAE 4.991 and RMSE 6.090. Both improve a little.
(c) Without, MAE 8.226 and RMSE 9.650. With, MAE 7.351 and RMSE 9.416. Both improve again, and the MAE improves by much more than it did in sample.
(d) Use the model with Days_Above_30. It is better on every measure on both sets of rows, which is the clearest case you can get for keeping a variable.
What it bought is mostly out of sample. In sample it lowered the MAE by 0.11 bu/ac; on the test years it lowered it by 0.87, about eight times as much.