ril_link <- "https://raw.githubusercontent.com/ybrandvain/datasets/refs/heads/master/clarkia_rils.csv"
rils <- readr::read_csv(ril_link) |>
dplyr::select(ril, prop_hybrid, petal_area_mm, asd_mm, location, petal_color) |>
na.omit()Motivating example: ASD is associated with the proportion of hybrid seed, but ASD is also associated with other floral traits. So, when we see that RILs with larger ASD have more hybrid seed, we want to know whether ASD itself is associated with hybrid seed proportion, or whether this association is partly explained by other traits that covary with ASD.
Learning goals: By the end of this section, you should be able to:
We aim to build a multiple regression model to ask if the association between ASD and proportion of hybrid seed holds up after accounting for associations between ASD and the other explanatory variables. But first, let’s review how linear models make predictions.
From a simple linear regression, we find the predicted proportion of hybrid seeds, \(\hat{Y_i}\), of a RIL with an ASD of \(\text{ASD}_{mm}\) is:
lm(prop_hybrid ~ asd_mm,
data = rils) | (Intercept) | asd_mm |
|---|---|
| 0.027 | 0.128 |
\[\hat{Y_i} = 0.027 + 0.128 \times \text{ASD}_{mm}\]
That is, we expect an additional 12.8% hybrid seed for each mm increase in anther-stigma distance. So, for a RIL with ASD of 1 mm, we expect a hybrid seed proportion of \(0.027 + 0.128 \times 1 = 0.155\), and we expect a hybrid seed proportion of \(0.027 + 0.128 \times 2 = 0.283\) for a RIL with ASD of 2 mm.
But, as we saw in the previous section, ASD is also associated with petal color and petal area:
A similarly simple linear model shows that white-petaled RILs have anther-stigma distances about 0.2 mm shorter than pink-petaled RILs.
lm(asd_mm ~ petal_color,
data = rils) | (Intercept) | petal_colorwhite |
|---|---|
| 0.968 | -0.203 |
\[ \widehat{ASD}_i = 0.968 - 0.203 \times \begin{cases} 0, & \text{if petal color is pink} \\ 1, & \text{if petal color is white} \end{cases} \]
That is, the average pink-petaled RIL has an ASD of 0.968, and the average white-petaled RIL has an ASD of 0.765.
From a simple linear regression, we find the ASD, \(\widehat{ASD_i}\), of a RIL with a petal area of \(\text{petal area}_{mm^2}\) is:
lm(asd_mm ~ petal_area_mm,
data = rils) | (Intercept) | petal_area_mm |
|---|---|
| 0.5281 | 0.00543 |
\[\widehat{ASD_i} = 0.5281 + 0.00543 \times \text{petal area}_{mm^2}\]
That is, we expect an additional 0.00543 ASD in mm for each \(mm^2\) increase in petal area. So, for a RIL with petal area of \(40\) \(mm^2\), we expect an ASD of \(0.5281 + 0.00543 \times 40 = 0.7453\) mm, and we expect an ASD of \(0.5281 + 0.00543 \times 100 = 1.0711\) mm for a RIL with a petal area of \(100\) \(mm^2\).
Now let’s see how adjusting for petal color or petal area influences the association between ASD and proportion of hybrid seed. To help us understand and interpret these models, we will initially consider each pair of traits before building a bigger model with more of the traits.
We can add petal color into the model with a + sign. We see that doing so nearly halves the increase in proportion of hybrid seed with increasing ASD:
lm(prop_hybrid ~ petal_color + asd_mm, data = rils)
Call:
lm(formula = prop_hybrid ~ petal_color + asd_mm, data = rils)
Coefficients:
(Intercept) petal_colorwhite asd_mm
0.16669 -0.16679 0.06332 We can use the coefficients to find the expected proportion of hybrid seed, \(\hat{Y_i}\), for a plant with a given petal color and anther-stigma distance (ASD). To do so, start with the intercept (0.167), then add the contribution of each explanatory variable. For example:
| ASD = 1 add 0.063 |
ASD = 2 add 0.126 |
|
|---|---|---|
| Pink add 0 |
\(\hat{Y} = 0.167 + 0.063\) \(= 0.23\) |
\(\hat{Y} = 0.167 + 0.126\) \(= 0.293\) |
| White subtract 0.167 |
\(\hat{Y} = 0.167 - 0.167 + 0.063\) \(= 0.063\) |
\(\hat{Y} = 0.167 - 0.167 + 0.126\) \(= 0.126\) |
We can similarly add petal area to the model with a + sign. After accounting for petal area, the estimated contribution of ASD to our prediction again decreases, although not by as much as it did after accounting for petal color:
lm(prop_hybrid ~ petal_area_mm + asd_mm, data = rils)
Call:
lm(formula = prop_hybrid ~ petal_area_mm + asd_mm, data = rils)
Coefficients:
(Intercept) petal_area_mm asd_mm
-0.132984 0.002957 0.100355 Again, we can use these coefficients to find the expected proportion of hybrid seed, \(\hat{Y}_i\), for a plant with a given petal area and ASD. To do so, start with the intercept, then add the contribution of each explanatory variable:
\[\hat{Y}_i = \text{intercept} + (\text{petal-area coefficient} \times \text{petal area}) + (\text{ASD coefficient} \times \text{ASD})\]
In our case:
\[ \hat{Y}_i = -0.133 + 0.003 \times \text{petal area} + 0.1 \times \text{ASD} \]
For example, we can find the predicted proportion of hybrid seed for RILs with a petal area of 50 or 100 mm^2 and an ASD of 1 or 2 mm as follows:
| ASD = 1 add \(0.1 \times 1 = 0.1\) |
ASD = 2 add \(0.1 \times 2 = 0.2\) |
|
|---|---|---|
| Petal area = 50 add \(0.003 \times 50 = 0.15\) |
\(\hat{Y} = -0.133 + 0.15 + 0.1\) \(= 0.117\) |
\(\hat{Y} = -0.133 + 0.15 + 0.2\) \(= 0.217\) |
| Petal area = 100 add \(0.003 \times 100 = 0.3\) |
\(\hat{Y} = -0.133 + 0.3 + 0.1\) \(= 0.267\) |
\(\hat{Y} = -0.133 + 0.3 + 0.2\) \(= 0.367\) |
There is nothing special about stopping at two explanatory variables. We can continue adding explanatory variables to our model with + signs.
For example, we can account for petal area and petal color when estimating the association between ASD and the proportion of hybrid seed:
lm(prop_hybrid ~ petal_area_mm + asd_mm + petal_color, data = rils)
Call:
lm(formula = prop_hybrid ~ petal_area_mm + asd_mm + petal_color,
data = rils)
Coefficients:
(Intercept) petal_area_mm asd_mm petal_colorwhite
-0.010664 0.003408 0.028254 -0.174958 Now the ASD coefficient describes the expected change in proportion of hybrid seed associated with a one-mm increase in ASD among RILs with the same petal area and petal color. Importantly, after adjusting for petal area and petal color, proportion of hybrid seed only modestly increases with ASD, going from a slope of about 0.13 in a model with ASD alone to a slope of about 0.03 in a model with the other two explanatory variables. This suggests that much of the increase in proportion of hybrid seed with increasing ASD can be explained by petal area and petal color.
As before, each coefficient tells us how much to add to our prediction. The only difference is that we now have three explanatory variables to account for.
\[ \hat{Y}_i = \text{intercept} + (\text{petal-area coefficient} \times \text{petal area}) + (\text{ASD coefficient} \times \text{ASD}) + \text{petal-color term} \]
In our case:
\[ \hat{Y}_i = -0.011 + 0.003 \times \text{petal area}+ 0.028 \times \text{ASD} - 0.175 \times \begin{cases} 0, & \text{if petal color is pink} \\ 1, & \text{if petal color is white} \end{cases} \]
The petal-color term is zero for a pink-petaled RIL and -0.175 for a white-petaled RIL. For example, the expected proportion of hybrid seed for a white-petaled RIL with a petal area of 50 mm^2 and an ASD of 1 mm is:
\[ \begin{aligned} \hat{Y} &= -0.011 + (0.003 \times 50) + (0.028 \times 1) + -0.175\\ &= -0.008. \end{aligned} \]
So our predictions now look like this:
| petal color | petal area | ASD | intercept | + petal area | + petal color | + ASD | predicted hybrid seed |
|---|---|---|---|---|---|---|---|
| pink | 50 | 1 | -0.011 | 0.150 | 0.000 | 0.028 | 0.167 |
| pink | 50 | 2 | -0.011 | 0.150 | 0.000 | 0.056 | 0.195 |
| pink | 100 | 1 | -0.011 | 0.300 | 0.000 | 0.028 | 0.317 |
| pink | 100 | 2 | -0.011 | 0.300 | 0.000 | 0.056 | 0.345 |
| white | 50 | 1 | -0.011 | 0.150 | -0.175 | 0.028 | -0.008 |
| white | 50 | 2 | -0.011 | 0.150 | -0.175 | 0.056 | 0.020 |
| white | 100 | 1 | -0.011 | 0.300 | -0.175 | 0.028 | 0.142 |
| white | 100 | 2 | -0.011 | 0.300 | -0.175 | 0.056 | 0.170 |
We can continue in the same way and add location to the model:
lm(prop_hybrid ~ location + petal_area_mm + petal_color + asd_mm, data = rils)
Call:
lm(formula = prop_hybrid ~ location + petal_area_mm + petal_color +
asd_mm, data = rils)
Coefficients:
(Intercept) locationLB locationSR locationUS
0.009550 0.073067 0.017549 -0.114336
petal_area_mm petal_colorwhite asd_mm
0.003218 -0.173831 0.030596 This model estimates the association between ASD and proportion of hybrid seed after accounting for location, petal area, and petal color. Again, a prediction begins with the intercept, and then we add the contribution of each explanatory variable:
\[ \hat{Y}_i = \text{intercept} + \text{location term} + (\text{petal-area coefficient} \times \text{petal area}) + \text{petal-color term} + (\text{ASD coefficient} \times \text{ASD}). \]
Rather than work through this somewhat burdensome math by hand, I have developed an app below that allows you to build a linear model, set trait values, and find \(\hat{Y_i}\):
Do not trust these predictions for variable combinations outside your data.
Prediction is not the same as causation. This app (or any multiple regression) lets us change values of petal area, petal color, ASD, and location and see how the predicted proportion of hybrid seed changes under a linear model. This is useful for understanding how multiple regression makes predictions.
But moving the sliders in the app (or changing values in a linear model) is not the same thing as doing an experiment.
In an experimental study, we manipulate an explanatory variable and compare outcomes across treatments. For example, if we could experimentally increase ASD while holding other relevant traits constant, then a difference in hybrid seed could support a causal claim about ASD.
In an observational study, we measure traits as they naturally occur. In that case, regression can tell us whether ASD is associated with hybrid seed after accounting for other measured variables, such as petal area, petal color, and location. But it cannot, by itself, prove that ASD causes a change in hybrid seed, because RILs that differ in ASD may also differ in other ways we did not measure or include in the model.
So, in this section, we should interpret regression coefficients as conditional associations: expected differences in predicted hybrid seed among RILs with the same values of the other variables in the model.
This is an observational study: we measured traits of RILs as they occurred. Multiple regression can ask whether ASD is associated with hybrid seed after accounting for other measured variables. But that is not the same as experimentally changing ASD while holding other relevant factors constant. So here we interpret coefficients as conditional associations, not necessarily causal effects.
The app can make predictions for combinations of variables that may be rare, or even absent, in the original data. For example, if very large petal area is mostly seen in one petal color or one location, then a prediction for that petal area in a different color or location may be extrapolating beyond the data that informed the model. Multiple regression can make a numerical prediction for such cases, but that does not mean the prediction is equally trustworthy.