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: We can now build linear models with multiple explanatory variables and use them to estimate the expected value of a response for an individual, conditional on its values for those explanatory variables. But how do we compare which explanatory variables are most strongly associated with the response when they are measured on different scales or represented in different ways? In this section, we introduce several ways to compare effect sizes across predictors while keeping track of what each comparison does, and does not, mean.
Learning goals: By the end of this section, you should be able to:
We care not just to model the world, and not just to test the null hypothesis that two variables are unassociated. Rather we want to know if an explanatory variable matters for our response. Ultimately, deciding the strength that “matters” is up to the decision maker – for example, a government might decide how much increase in cancer risk is acceptable when putting in a new power plant or whatever, or in the context of our Clarkia data, an evolutionary theorist might know the maximum amount of gene flow two populations can experience before they merge into one species.
For both the examples above, the key measure of effect was simply the coefficient in the multiple regression model. This measure of effect is very handy – especially in the hands (or brains) of experts, who can make use of this information.
Our model coefficients are:
| term | coefficient |
|---|---|
| (Intercept) | -0.0107 |
| petal_area_mm | 0.0034 |
| petal_colorwhite | -0.1750 |
| asd_mm | 0.0283 |
lm(prop_hybrid ~ petal_area_mm + petal_color+ asd_mm, data = rils) |>
coef()But say I wanted to know which floral phenotype was most strongly associated with the proportion of hybrid seeds, after adjusting for the other variables. For example,
The coefficient for ASD, measured in \(mm\), is more than eight times greater than the coefficient for petal area, measured in \(mm^2\). Does this mean that ASD is nearly eight times more important?
Not necessarily. Petal area ranges from 34 to 112 \(mm^2\), while ASD ranges from 0.39 to 2.15 \(mm\).
As a first attempt, we could compare the predicted difference between the smallest and largest observed value of each trait. In that case, the model predicts a difference of 0.268 in the proportion of hybrid seeds between the RILs with the largest and smallest petal areas. It predicts a difference of 0.05 between the RILs with the largest and smallest anther–stigma distances. This “difference in extremes” is more informative than simply comparing two coefficients. However, because it focuses on extremes, this approach is not a perfect summary of effect size either.
Similarly, how do we compare the importance of petal color vs petal area?
We saw that our model predicts a difference of 0.268 in the proportion of hybrid seed between RILs with the largest and smallest petals. How do we compare this with the expected difference between pink- and white-flowered RILs?
In this subsection, I introduce methods to make such fair comparisons.
Effect sizes are great, but don’t let the name fool you.
| term | coefficient |
|---|---|
| (Intercept) | 0.378 |
| scale(petal_area_mm) | 0.221 |
| scale(asd_mm) | 0.044 |
| petal_colorwhite | -0.752 |
One way to compare the “effect” of continuous explanatory variables is to Z-transform them (as well as the response variable). After doing so, we can interpret linear model coefficients as the increase in the response, measured in standard deviations, for every standard deviation increase in a continuous explanatory variable.
lm(scale(prop_hybrid) ~ scale(petal_area_mm) + scale(asd_mm) + petal_color, data = rils)These “standardized regression coefficients” are useful because they can be directly compared between continuous explanatory variables. So we can now easily see the stronger standardized association between petal area and the proportion of hybrid seeds than between ASD and the proportion of hybrid seeds.
However, standardized coefficients have some drawbacks:
First, the coefficients for continuous and categorical explanatory variables represent different kinds of comparisons: a one-standard-deviation increase for a continuous variable versus a difference between categories.
Second, although standardization does not require normally distributed variables, one standard deviation may not represent a natural or typical change when an explanatory variable is strongly skewed or contains extreme values.
Finally, it is not always natural or useful to think in standard deviation units.
This doesn’t mean we should never use standardized regression coefficients. It simply means that we should know when they are and are not appropriate, and how to cautiously present and interpret them.
Above, I Z- transformed the continuous response and continuous predictors, but did not change the categorical predictor. This still changes the meaning (and value) of the coefficient(s) for the categorical predictor(s). Now the petal color coefficient describes the number of standard deviations that separate the proportion of hybrid seeds on pink- and white-flowered RILs.
Perhaps the most fair way to compare the “effect size” of a mix of categorical and continuous explanatory variables is the partial \(R^2\) (\(R^2_\text{partial}\)).
The partial \(R^2\) describes the proportion of variation in the response variable “explained” by a focal explanatory variable. This is equal to one minus the ratio of \(SS_\text{error}\) in a model with and without the focal explanatory variable.
\[R^2_\text{partial} = 1-\frac{SS_\text{error (full model)}}{SS_\text{error (model without focal predictor)}}\]
Because partial \(R^2\) adjusts for other variables in the model, \(\sum R^2\) does not generally add up to the \(R^2\) of the full model.
The drop1() function compares the full model to a series of reduced models, each missing one model term. It returns various summaries of each “reduced” model (as well as the full model in the first row, labeled <none>, which represents the full model). For our purposes, we will focus on the residual sums of squares (aka \(SS_\text{error}\) aka \(SS_\text{residual}\)), RSS. To show you that this works for categorical predictors with more than two levels, I will add location into the model.
full_model <- lm( prop_hybrid ~ petal_area_mm + asd_mm + petal_color +location, data = rils)
drop1(full_model)Single term deletions
Model:
prop_hybrid ~ petal_area_mm + asd_mm + petal_color + location
Df Sum of Sq RSS AIC
<none> 15.302 -1299.9
petal_area_mm 1 0.88400 16.186 -1279.4
asd_mm 1 0.04282 15.344 -1300.8
petal_color 1 2.77770 18.079 -1234.9
location 3 1.94348 17.245 -1257.9We can find the partial \(R^2\) as:
drop1(full_model)|>
dplyr::select(RSS )|>
as_tibble(rownames = "term")|>
mutate(partial_r2 = 1 - min(RSS) / RSS)# A tibble: 5 × 3
term RSS partial_r2
<chr> <dbl> <dbl>
1 <none> 15.3 0
2 petal_area_mm 16.2 0.0546
3 asd_mm 15.3 0.00279
4 petal_color 18.1 0.154
5 location 17.2 0.113 This result makes it clear that, after adjusting for other variables in the model, petal color and planting location explain most of the variation in the proportion of hybrid seeds on a RIL (15% and 11%, respectively), with petal area explaining another 5% of this variation. By contrast ASD explains less 0.3% of the variation in the proportion of hybrid seed.
partial_r2() in honestlmI added a partial_r2() function to my honestlm package to calculate \(R^2_\text{partial}\) from an lm object. The syntax is as follows:
library(honestlm)
partial_r2(full_model)# A tibble: 4 × 4
term df partial_r2 f2
<chr> <dbl> <dbl> <dbl>
1 petal_area_mm 1 0.0546 0.0578
2 asd_mm 1 0.00279 0.00280
3 petal_color 1 0.154 0.182
4 location 3 0.113 0.127 This function also returns Cohen’s \(f^2\) a related measure of effect size. Cohen’s \(f^2\) is just a transformation of partial \(r^2\), where \(f^2 = R^2_\text{partial} / (1 - R^2_\text{partial})\), so larger values mean a larger term-level contribution.
Just like a traditional \(R^2\), a partial \(R^2\) is only interpretable in the study we are conducting. We could not say, for example, that “petal area explains 5% of the variation in proportion hybrid seed set.” We could only say “In this study petal area explains 5% of the variation in proportion hybrid seed set.”