In the previous module, we examined a univariate regression with a response variable and a single explanatory variable. But most response variables cannot be explained by just one factor; yield does not depend on just precipitation, but also on temperature, on moisture, on soil, on the variety and seeding date. We could run separate regressions for each of these, but there are several issues with this approach. First, it takes time. Second – as we learned in the last chapter – if these factors are correlated with each other, then the resulting coefficients will be confounded. Finally, we would end up with a much worse prediction of yield than if we had used all of the factors together. The solution is to use a multivariate regression that includes all of the factors at once.
The data for this chapter is rm_canola_weather.csv. Each row is one rural municipality (RM) in one year – 295 RMs over the ten years from 2015 to 2024, so 2,950 rows. Saskatchewan’s RMs are the local government units the province’s crop statistics are reported by, which makes them a convenient unit for this kind of question even though they are nobody’s field.
The file has already been assembled for you out of four sources: reported crop yields, a lookup table naming each RM, weather from the nearest station, and the coordinates and soil zone of each RM. Putting it together is ordinary join work, and doing it again here would be a distraction from the modelling. The columns are:
Column
What it holds
Year, RM, RM_Name
which RM, which year
Census_Division
the larger region the RM sits in
latitude, longitude
the centre of the RM, in degrees
soil_zone
Brown, Dark Brown or Grey/Black
Population_2021
people living in the RM
Weather_Station, Station_km
the station the weather came from, and how far away it is
Precip_May_Aug_mm
growing-season rainfall, in millimetres
GDD_base5
growing degree days above 5 °C, accumulated over the season
> rm_weather <- read_csv("data/rm_canola_weather.csv")
Rows: 2950 Columns: 20
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): RM_Name, Census_Division, soil_zone, Weather_Station
dbl (16): Year, RM, latitude, longitude, Population_2021, Station_km, Precip...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
>
> dim(rm_weather)
[1] 2950 20
Two things about this data are worth knowing before we use it. The weather is measured at the nearest station rather than in the RM itself, a median of 62 km away, so it describes the RM approximately. And a yield column is blank where a crop was not grown in that RM that year, which is why canola has 2,904 usable rows rather than 2,950.
26.1 Adding additional predictors
A multiple regression puts several predictors in the same model:
\[
y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \ldots + \varepsilon_i
\] where \(y_i\) is the response for observation \(i\), \(x_{1i}\) is the first predictor, \(x_{2i}\) is the second predictor, and so on. The \(\beta\) coefficients are the slopes for each predictor, and \(\varepsilon_i\) is the error term.
For example, let’s consider trying to explain the average canola yield in rural municipalities (RMs) in Saskatchewan, using weather and other data. Let’s start with just including growing-season rainfall as a predictor in a univariate regression.
# canola yield against May-August rainfallcanola <-filter(rm_weather, !is.na(Canola))rain_only <-lm(Canola ~ Precip_May_Aug_mm, data = canola)summary(rain_only)
> # canola yield against May-August rainfall
> canola <- filter(rm_weather, !is.na(Canola))
> rain_only <- lm(Canola ~ Precip_May_Aug_mm, data = canola)
> summary(rain_only)
Call:
lm(formula = Canola ~ Precip_May_Aug_mm, data = canola)
Residuals:
Min 1Q Median 3Q Max
-34.698 -5.647 0.680 6.768 27.638
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.573471 0.501642 46.99 <2e-16 ***
Precip_May_Aug_mm 0.057027 0.002332 24.45 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.562 on 2902 degrees of freedom
Multiple R-squared: 0.1708, Adjusted R-squared: 0.1705
F-statistic: 597.8 on 1 and 2902 DF, p-value: < 2.2e-16
Table 26.1: Canola yield on rainfall alone.
Model 1: rainfall only. (adjusted \(R^2\) = 0.1705, n = 2904)
Term
Estimate
Std. Error
t value
p value
Intercept
23.5735
0.5016
46.99
<0.001
May to Aug precip (mm)
0.0570
0.0023
24.45
<0.001
The slope is 0.057. Taken literally, that says an RM-year with 100 mm more rain yields about 5.7 bu/ac more canola, and the model (including only rainfall) accounts for about 17% of the variation. To see the size of the effect, consider that the average rainfall in the data is about 200 mm with a standard deviation of 75 mm. Meaning that we can easily see rainfall in RMs differing by as much as 150 mm (one standard deviation above and below the mean).
But consider what happens when we add the number of days above 30 °C to the model. We do this in R by just adding the second predictor to the formula with a plus sign:
# rainfall and the number of days above 30 degrees, togetherboth <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30, data = canola)summary(both)
> # rainfall and the number of days above 30 degrees, together
> both <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30, data = canola)
> summary(both)
Call:
lm(formula = Canola ~ Precip_May_Aug_mm + Days_Above_30, data = canola)
Residuals:
Min 1Q Median 3Q Max
-27.2977 -4.5641 0.1984 4.9660 23.7216
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.522568 0.607917 73.238 <2e-16 ***
Precip_May_Aug_mm 0.004996 0.002141 2.333 0.0197 *
Days_Above_30 -0.873475 0.019571 -44.632 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.364 on 2901 degrees of freedom
Multiple R-squared: 0.5084, Adjusted R-squared: 0.508
F-statistic: 1500 on 2 and 2901 DF, p-value: < 2.2e-16
Table 26.2: Canola yield on rainfall and the count of days above 30 degrees.
Model 2: rainfall and hot days. (adjusted \(R^2\) = 0.508, n = 2904)
Term
Estimate
Std. Error
t value
p value
Intercept
44.5226
0.6079
73.24
<0.001
May to Aug precip (mm)
0.0050
0.0021
2.33
0.020
Days above 30 degrees
-0.8735
0.0196
-44.63
<0.001
With the count of hot days added to the model, the rainfall coefficient drops to 0.005 – about one tenth of its previous size. Holding rainfall constant, one more day above 30 °C is associated with about 0.87 bu/ac lower canola yield.
What happened to the coefficient on rainfall? The simple regression we ran first was confounded by the fact that wet years are also cool years. So precipitation was not just capturing the effect of additional moisture, but also the fact that there were fewer hot days in those years. Once we added the number of days above 30 °C to the model, the effect of rainfall dropped precipitously.
Of course, we can add even more to our model. For example, both precipitation and heat might be affected by latitude. Perhaps it is just that northern RMs are both cooler and have better soil for canola. We can add latitude to the model and see what happens.
The latitude column is already in the data, but it is worth knowing where it came from, because an RM does not have a latitude – it is an area, not a point. The number in the file is the latitude of the RM’s centroid, the centre of its boundary polygon, computed from the sask_rms.geojson file used in Chapter 19:
Now we can add latitude to the model and see what happens.
# rainfall, hot days and latitudewith_lat <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude,data = canola)
> # rainfall, hot days and latitude
> with_lat <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude,
+ data = canola)
Table 26.3: Canola yield on rainfall, hot days and latitude.
Model 3: adding latitude. (adjusted \(R^2\) = 0.5202, n = 2904)
Term
Estimate
Std. Error
t value
p value
Intercept
-8.8835
6.2100
-1.43
0.153
May to Aug precip (mm)
0.0081
0.0021
3.77
<0.001
Days above 30 degrees
-0.7906
0.0216
-36.64
<0.001
Latitude
1.0113
0.1170
8.64
<0.001
Latitude has a coefficient of 1.01, so an RM one degree further north yields about a bushel more canola. The growing area of the province spans five-degrees of latitude. Taken literally, the coefficient says that yields in the northernmost RM are about 5 bu/acre greater than those in the southernmost RM, holding rainfall and hot days constant.
The inclusion of latitude also shrinks the coefficient on days above 30 a little bit, from -0.87 to -0.79. This is because the south is both hotter and worse for canola for reasons that are not heat, and the hot-day count had been carrying some of that. The rainfall coefficient goes the other way, from 0.005 to 0.008.
But maybe latitude itself is just a proxy for something else, such as soil. Saskatchewan’s agricultural soils are usually described in bands running from the dry southwest to the moister northeast: Brown in the south, Dark Brown through the middle, and the Grey/Black soils of the north and east. Organic matter rises as you move into darker soil zones.
The data has a soil_zone column. But how do we include it in a regression? It is not a number, but a text label. The answer is that R treats it as a categorical predictor. This creates a separate column for each soil zone – Brown, Dark Brown and Grey/Black. An RM in the Brown soil zone has a value of 1 in the Brown column and 0 in the other two. We then include these binary or dummy variables in the regression. However, we must exclude one of the soil zones. This becomes the omitted group against which we compare the other soil zones. This may be best illustrated by example. Here is the regression with both latitude and the soil zones:
# rainfall, hot days, latitude and the soil zone categorieswith_both <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,data = canola)
> # rainfall, hot days, latitude and the soil zone categories
> with_both <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone,
+ data = canola)
Table 26.4: Canola yield on rainfall, hot days, latitude and soil zone together.
Model 4: latitude and soil zone together. (adjusted \(R^2\) = 0.53, n = 2904)
Term
Estimate
Std. Error
t value
p value
Intercept
26.4649
8.0503
3.29
0.001
May to Aug precip (mm)
0.0077
0.0021
3.61
<0.001
Days above 30 degrees
-0.7296
0.0227
-32.10
<0.001
Latitude
0.2612
0.1582
1.65
0.099
Soil zone: Dark Brown
2.7251
0.3980
6.85
<0.001
Soil zone: Grey/Black
4.0360
0.5358
7.53
<0.001
You can see that R does not give us a coefficient for the brown soil zone – we can think of this coefficient as essentially zero. The positive coefficients for the dark brown and grey/black soil zones tell us that, holding rainfall, hot days and latitude constant, yields in these zones are higher than those in the brown zone. Yields in the dark brown zone are 2.7 bu/ac better than the brown zone, and yields in the grey/black zone are 4.0 bu/ac better.
We can also see that the latitude coefficient has dropped from 1.01 to 0.26, and its p-value has risen to 0.099. This suggests that the effect that we saw previously is not about latitude itself, but was rather more about the fact that northern RMs are more likely to be in the dark brown and grey/black soil zones.
I realize that the above discussion involved a lot of numbers, and it is easy to get lost in them. The key point is that the coefficients in the model are sensitive to what other predictors are included, as many of them are correlated with each other. A variable that looks useful can stop looking useful when a better measure of the same underlying cause joins the model. Hopefully, this highlights that it is very important to carefully think about what to include in a model and even more important to carefully interpret the results of a model.
A regression is essentially only telling us correlation – it is up to us to interpret this correlation.
Interpreting regression coefficients
The standard phrase for interpreting a coefficient in a multiple regression is holding all else constant, and it is worth being precise about what it means, because it is easy to read more into it than it says.
Take the rainfall coefficient in Model 4, which is 0.0077. The full reading is: comparing two RM-years that had the same number of hot days, at the same latitude, in the same soil zone, the one with a millimetre more rain yielded 0.0077 bu/ac more canola. Every predictor in the model gets held fixed; every predictor not in the model does not.
We should note that this is a model-based comparison. There are unlikely to be two RM-years in the data matching exactly on hot days, latitude and soil zone and differing only in rainfall; the model estimates the comparison from all the data at once rather than from matched pairs. That is ordinary and it is what regression is for. It becomes less credible when the data hold little information about the combination being described – a very wet RM-year in the Brown zone, say – because then the comparison rests on the straight-line form rather than on observations.
Or think about the latitude coefficient, which is 0.2612. The reading is: comparing two RM-years in the same soil zone, with the same rainfall and the same number of hot days, the one a degree further north yielded 0.2612 bu/ac more canola. Note that soil zones do span several degrees of latitude, so it is certainly possible to find two RM-years that differ in latitude and are in the same soil zone. But this is different than saying that all RMs that are further north are better for canola than those in the south – as there are some RMs in the black soil zone in the southeast of the province, and some brown soil zones that are somewhat more northern in the west of province.
To continue beating this point home: the phrase also says nothing about causation. A coefficient tells you about a comparison between observations, not about what would happen if you intervened. If it rained one more millimetre on a given RM, this model does not promise that if you applied another mm of water that you would get 0.0077 bu/ac more canola. There may be other factors that aren’t in the model that correlate with rainfall. The claim we can make is much more modest: RM-years that differ by a millimetre, and match on the other explanatory variables included in the model, differ 0.0077 bu/ac, on average.
Prediction in multivariate models
One of the benefits of including multiple predictors in a model is that it allows us to make better predictions about the response variable.
To predict yields based on our model we need values for every explanatory variable in the model. We then simply put these values into the model equation along with the estimated coefficients. For example, in our ultimate yield model, we would predict RM yields as:
where the zone is 0 for Brown, 2.73 for Dark Brown and 4.04 for Grey/Black.
Take an RM at 52 °N in the Dark Brown zone, in a season with 250 mm of rain and 8 days above 30 °C. Multiply each coefficient by its value and add:
Term
Calculation
Contribution
Intercept
26.46
Rainfall
\(0.0077 \times 250\)
1.92
Hot days
\(-0.7296 \times 8\)
-5.84
Latitude
\(0.2612 \times 52\)
13.58
Soil zone
Dark Brown
2.73
Predicted yield
38.85
In R the predict() function does this for you. Give it the model and a data frame holding one row per prediction, with a column for every predictor in the model, named exactly as it is named in the data:
predict() handles more than one row at a time, which is the usual way to compare scenarios. Here are a hot dry southern RM and two cooler northern ones that differ only in soil zone:
The southern RM comes out at 26.0 bu/ac and the northern ones at 38.9 and 40.2. The last two differ only in soil zone, so the 1.3 bu/ac between them is the difference between the two zone coefficients – \(4.036 - 2.725 = 1.311\) – and nothing else. Both coefficients are measured against Brown, the level R held out, so neither one on its own is the Grey/Black-to-Dark-Brown comparison; you get that by subtracting.
Three things will trip up a prediction. The data frame needs a column for every predictor, so leaving one out is an error rather than a default. The names must match the model’s, including capitalisation. And a categorical predictor must be spelled exactly as it appears in the data – "Dark Brown", not "dark brown" – because R matches it against the factor levels and an unrecognised label gives an error.
The arithmetic will also answer questions it has no business answering. Ask this model about an RM with 1,000 mm of rain and it will return a number, even though the wettest RM-year in the data had 470 mm.
26.2 How to build a model
Now that we know how to include multiple explanatory variables in a model, the question becomes: which ones should we include?
The tempting answer is everything available. It is the wrong answer, and seeing why is most of what this section is about.
Unfortunately, some practitioners just throw a bunch of explanatory variables into a model and see which lands as statistically significant. But this is not a good approach. To see why, consider what happens when we add an RM’s population to our yield model.
with_pop <-lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude + soil_zone + Population_2021, data = canola)summary(with_pop)$coefficients["Population_2021", ]summary(with_pop)$adj.r.squared
> with_pop <- lm(Canola ~ Precip_May_Aug_mm + Days_Above_30 + latitude +
+ soil_zone + Population_2021, data = canola)
>
> summary(with_pop)$coefficients["Population_2021", ]
Estimate Std. Error t value Pr(>|t|)
0.0004447167 0.0001935617 2.2975444714 0.0216586909
> summary(with_pop)$adj.r.squared
[1] 0.5306611
Population is statistically significant at the 5% level. But it likely is just correlated with yields, not a causal factor. There are many stories we could tell about why population might be correlated with yield. For example, land with good soil may have attracted early settlers, and those settlers may have built towns that grew into the cities we see today. But this doesn’t mean that adding population to RMs will help yields.
Generally, you should include the variables that you think could have an effect on the response variable. Sometimes we might also include things that are proxies for variables that we do not have a direct measure of. So for example, we might include fertilizer prices not because higher prices literally cause lower yields, but because we might not have data on fertilizer use, and higher prices might be correlated with lower use.
Of course, we don’t have access to all the variables that might affect yields (or other outcomes), so we sometimes have to make do with what we have, or spend time and money collecting other variables. One question is how hard should we work to find these variables.
If our goal is explaining the impact of explanatory variables, then we should work very hard to include the variables that might correlate with our explanatory variables of interest. For example, if we want to know the effect of nitrogen on yields, and we have a dataset of the amount of nitrogen producers have applied, then we should also try to get data on the other production practices of farmers and the soil quality of their fields, as these might also affect yields and be correlated with nitrogen use.
Conversely, if we want to know the impact of rainfall on yields, then these management practices may be less important. Not because they don’t affect yields, but because they are less likely to be correlated with rainfall.
If our goal is prediction, we want to work hard to include variables that carry useful information about yields – weather, management practices, soil quality, and so on. But note that “improves the fit” is not the test: as Section 28.4 shows, adding variables can only improve the fit on the data the model was built from, whether or not those variables help. A variable earns its place in a prediction model by improving performance on data the model has not seen.