library(dplyr)
library(readr)
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 scenario: You ran a model with interactions, estimated parameters, and even evaluated significance. Now what? The results are hard to think about and explain. To understand the biology, I show you how to translate this model into meaningful predictions, estimates, and measures of explained variation.
Learning goals: By the end of this section, you should be able to:
We modelled the proportion of hybrid seeds on a RIL as a function of petal area, petal color, and their interaction. If you recall there was some weirdness - the model returned a large and positive coefficient for petal_colorwhite, even though we found no region of our data in which the model predicted that white-petaled RILs would have a higher proportion of hybrid seeds than pink-petaled RILs (Figure 1 A).
We realized that this was because the coefficient of 0.2373 for petal_colorwhite describes the difference between its expected value at the intercept (i.e. petal_area_mm = 0, Figure 1 B) and that of the pink-petled RIls. This prediction with no biological application as no plants have petals of size zero \(mm^2\).
| name | value |
|---|---|
| (Intercept) | -0.2356 |
| petal_colorwhite | 0.2373 |
| petal_area_mm | 0.0075 |
| petal_colorwhite:petal_area_mm | -0.0067 |
So, if we can’t get simple biological meaning from raw model coefficients, how can we make sense of our models? My best advice is to generate meaningful estimates of predicted means + uncertainty in them, for a few illustrative values of the explanatory variables.
Let’s start by doing this ourselves, by simply plugging in petal areas of 45 and 85 \(mm^2\) for white- and pink- petaled RILs into the model with interactions. I will also provide predictions from a linear model with no interactions for reference. Click the section below if you want to see this math:
model_interaction <- lm(prop_hybrid ~ petal_color * petal_area_mm,data = rils)Interaction model:
Call:
lm(formula = prop_hybrid ~ petal_color * petal_area_mm, data = rils)
Coefficients:
(Intercept) petal_colorwhite
-0.235614 0.237309
petal_area_mm petal_colorwhite:petal_area_mm
0.007477 -0.006729 Interaction predictions:
# Color: pink; petal area: 45 mm^2
-.2356 + 45 * .007477 [1] 0.100865# Color: white; petal area: 45 mm^2
-.2356 + 0.237309 + 45 * (.007477-0.006729) [1] 0.035369 # Color: pink; petal area: 85 mm^2
-.2356 + 85 * .007477 [1] 0.399945# Color: white; petal area: 85 mm^2
-.2356 + 0.237309 + 85 * (.007477-0.006729)[1] 0.065289Additive model:
model_additive <- lm(prop_hybrid ~ petal_color + petal_area_mm, data = rils)
Call:
lm(formula = prop_hybrid ~ petal_color + petal_area_mm, data = rils)
Coefficients:
(Intercept) petal_colorwhite petal_area_mm
0.007045 -0.180748 0.003563 # Color: pink; petal area: 45 mm^2
0.007045 + 45 * 0.003563 [1] 0.16738# Color: white; petal area: 45 mm^2
0.007045 - 0.180748 + 45 * 0.003563 [1] -0.013368# Color: pink; petal area: 85 mm^2
0.007045 + 85 * 0.003563 [1] 0.3099# Color: white; petal area: 85 mm^2
0.007045 - 0.180748 + 85 * 0.003563 [1] 0.129152| Petal area | Petal color | Interaction model | No-interaction model |
|---|---|---|---|
| 45 | pink | 0.101 | 0.167 |
| 45 | white | 0.035 | -0.013 |
| 85 | pink | 0.400 | 0.310 |
| 85 | white | 0.065 | 0.129 |
In the interaction model, the predicted proportion of hybrid seed increases strongly with petal area among pink-petaled RILs but changes little among white-petaled RILs (Figure 2 A). So the predicted difference between colors is much larger at \(85 mm^2\) than at 45 \(mm^2\). Because the additive model predicts the same increase with petal area for both petal colors, the predicted difference between pink- and white-petaled RILs is constant across petal areas (Figure 2 B).
\(R^2\), the proportion of variation explained is a standard summary of “effect size.” In the section on multiple regression, we saw that we could calculate partial \(R^2\) as:
\[R^2_\text{partial} = 1-\frac{SS_\text{error (full model)}}{SS_\text{error (model without focal predictor)}}\]
For an interaction model we do the same thing, first for a main effects model:
library(honestlm)
lm(prop_hybrid ~ petal_color + petal_area_mm, data = rils)|>
partial_r2()# A tibble: 2 × 4
term df partial_r2 f2
<chr> <dbl> <dbl> <dbl>
1 petal_color 1 0.160 0.190
2 petal_area_mm 1 0.0625 0.0667We see that petal color and petal area “explain” 16% and 6% of the variance in the proportion of hybrid seeds on a RIL, respectively, after adjusting for the other variable.
And then for the interaction:
library(honestlm)
lm(prop_hybrid ~ petal_color * petal_area_mm, data = rils)|>
partial_r2()# A tibble: 1 × 4
term df partial_r2 f2
<chr> <dbl> <dbl> <dbl>
1 petal_color:petal_area_mm 1 0.0579 0.0614After adjusting for these main effects, the interaction between petal area and petal color explains 5% of the remaining variance in the proportion of hybrid seeds. Thus this interaction explains nearly as much variation in the proportion of hybrid seeds as petal area itself!!!
Of course, no quick statistical descriptions is complete without a measure of uncertainty. Again we can accomplish this with the emmeans or modelbased packages. Here I use estimate_means() function in the modelbased package, because I find it easier to code.
So to find uncertainty about the point estimates of the proportion of hybrid seeds for pink- and white-petaled RILs with petal aread of 45 and 85 \(mm^2\), we type:
library(modelbased)
estimate_means(
lm(prop_hybrid ~ petal_color * petal_area_mm, data = rils),
by = list(
petal_color = c("pink", "white"),
petal_area_mm = c(45, 85)))| petal_color | petal_area_mm | Mean | SE | CI_low | CI_high | t | df |
|---|---|---|---|---|---|---|---|
| pink | 45 | 0.101 | 0.023 | 0.056 | 0.145 | 4.439 | 398 |
| pink | 85 | 0.400 | 0.028 | 0.345 | 0.455 | 14.373 | 398 |
| white | 45 | 0.035 | 0.021 | -0.006 | 0.076 | 1.696 | 398 |
| white | 85 | 0.065 | 0.025 | 0.017 | 0.113 | 2.661 | 398 |
Here, each row is a model-predicted mean at one color–area combination, the confidence interval describes uncertainty around that predicted mean.
We can also report uncertainty in the slope for each group with the estimate_slopes() function in the modelbased package:
library(modelbased)
estimate_slopes(
lm(prop_hybrid ~ petal_color * petal_area_mm, data= rils),
slope = "petal_area_mm",
by = list(petal_color = c("pink", "white")))| petal_color | Slope | SE | CI_low | CI_high | t | df | p |
|---|---|---|---|---|---|---|---|
| pink | 0.0075 | 0.0010 | 0.0054 | 0.0095 | 7.2041 | 398 | 0.000 |
| white | 0.0007 | 0.0009 | -0.0010 | 0.0025 | 0.8498 | 398 | 0.396 |
The two confidence intervals describe uncertainty in each slope. Note that they do not directly test the difference between slopes, as we do this by testing the interaction.
Let’s try to synthesize these results into a brief “results” paragraph:
We tested the contribution of petal color, petal area and their interaction to the proportion of hybrid seeds on a RIL. Using Type II sums of squares, we found significant associations between the proportion of hybrid seeds and petal color (\(F_{1,398} = 80.2\), p < 2.2e-16), petal area (\(F_{1,398} = 28.1\), p < 2 e-7), and their interaction (\(F_{1,398} = 24.5\), p < 2 e-6). This interaction reflects the sharp increase in the proportion of hybrid seed with petal area in pink-, but not white-petaled RILs (compare slopes of 0.0075 – 95% CI: {0.0054 : 0.0095} and 0.0007 – 95% CI: {-0.0010 : 0.0025}, respectively).
.