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: You built a multiple regression model but now want to get onto the dirty business of null hypothesis significance testing.
Learning goals: By the end of this section, you should be able to:
car::Anova() to test model terms in a multiple regression.anova(),to find p-values for terms in a multiple regression, and explain why these are wrong.summary() and tidy() as term-level tests, and explain why these can be misleading.The point of multiple regression is to model a response variable as a function of several explanatory variables in one coherent equation, rather than trying to reason our way through a pile of separate models. With multiple regression, we can ask about the association between one explanatory variable and the response while accounting for the others in the model.
But how do we calculate p-values for each term in our model? The answer comes down to how we break up the sums of squares for the overall model into the specific terms in the model. Type II sums of squares – a common approach, and the one I advocate here – has a logic similar to the partial \(R^2\) calculation from the previous section. For each term in the model, we compare the sums of squares from the full model to a reduced model that contains all the other terms but leaves out the focal term.
I think about this as asking the question “how much extra variation in the response is captured when a given term is added to the model?”:
You may also see the same idea written using residual sums of squares:
\[SS_\text{term} = SS_\text{residual(model without term)}-SS_\text{residual(full model)}\]
In the (relatively simple) linear models in this and the next chapter, both equations give the same answer.
\[SS_\text{term} =SS_\text{model(full model)}-SS_\text{model(model without term)}\]
summary() for NHSTThe p-values from summary() (as well as broom’s tidy() function) test the significance of model coefficients away from zero. These are not the tests we want and these p-values will often mislead you.
Let’s walk through these calculations with broom’s augment() function. Recall that augment() takes a linear model and provides detailed info for each observation. Most importantly for this effort, it provides \(\hat{Y_i}\) in the column, .fitted. For this example we will focus on asd_mm
library(broom)
full_model <- lm(prop_hybrid ~ asd_mm + petal_color + petal_area_mm + location, data = rils )
reduced_model <- lm(prop_hybrid ~ petal_color + petal_area_mm + location, data = rils )
ss_model_full <- augment(full_model) |>
summarize(ss_model_full = sum((.fitted - mean(prop_hybrid))^2))|>
pull()
ss_model_reduced <- augment(reduced_model) |>
summarize(ss_model_reduced = sum((.fitted - mean(prop_hybrid))^2))|>
pull()
ss <- tibble(ss_model_full , ss_model_reduced ) |>
mutate(ss_asd = ss_model_full - ss_model_reduced)| ss_model_full | ss_model_reduced | ss_asd |
|---|---|---|
| 6.376 | 6.334 | 0.043 |
ss_error_full <- augment(full_model) |>
summarize(ss_error_full = sum((.fitted - prop_hybrid)^2))|>
pull()
ss_asd <- pull(ss,ss_asd )
df_asd <- 1
ms_asd <- ss_asd / df_asd
ss_error <- ss_error_full
df_error <- df.residual(full_model)
ms_error <- ss_error / df_error
f_val <- ms_asd / ms_error
p_val <- pf(f_val, df_asd, df_error,lower.tail = FALSE)
tibble(ss_asd, df_asd, ms_asd, ss_error, df_error, ms_error, f_val, p_val)| ss_asd | df_asd | ms_asd | ss_error | df_error | ms_error | f_val | p_val |
|---|---|---|---|---|---|---|---|
| 0.043 | 1 | 0.043 | 15.302 | 395 | 0.039 | 1.105 | 0.294 |
So, with a p-value of 0.294 we do not have strong evidence that anther-stigma distance is associated with proportion of hybrid seeds after adjusting for petal area, petal color, and location. This is despite a very low p-value (\(p \approx 0.00005\)) for the simple linear regression: lm(prop_hybrid ~ asd_mm, data = rils).
Because sums of squares for each term come after adjusting for the other terms, the full model’s model sum of squares will usually not equal the sum of the sums of squares for each term in the model.
Anova() function can use Type II Sums of SquaresThe section above was meant to explain how to calculate Type II sums of squares from a multiple regression and then get a p-value. But it would be a tremendous pain to do this for all terms in our model. Instead you can use the Anova() with a capital A function from the car package to do this for you.
library(car)
lm(prop_hybrid ~ asd_mm + petal_color + petal_area_mm + location, data = rils ) |>
Anova(type = "II")Anova Table (Type II tests)
Response: prop_hybrid
Sum Sq Df F value Pr(>F)
asd_mm 0.0428 1 1.1053 0.2937
petal_color 2.7777 1 71.7046 4.970e-16 ***
petal_area_mm 0.8840 1 22.8200 2.512e-06 ***
location 1.9435 3 16.7233 3.016e-10 ***
Residuals 15.3016 395
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We now see that after adjusting for other variables, all but anther-stigma distance remain strongly associated with the proportion of hybrid seeds.
anova() for multiple regressionBase R’s anova() function uses Type I sums of squares. In multiple regression, this means the p-value for a term can depend on where that term appears in the formula. This is usually not what we want (see below).
In sequential sums of squares, sums of squares for the first term are calculated as if it were the only term in the model. Then sums of squares for the second term are calculated adjusting for the first term. The sums of squares for the third term adjust for the first two terms, etc etc.
This means that if we accidentally use Type I sums of squares (which R does when we use the anova(), the p-value for a term can literally depend on the order it is entered into the lm() function.
So, for example, with Type I sums of squares we can wrongly conclude that anther-stigma distance is associated with proportion of hybrid seeds (if it is entered first)
lm(prop_hybrid ~ asd_mm + petal_color + petal_area_mm + location, data = rils ) |>
anova()Analysis of Variance Table
Response: prop_hybrid
Df Sum Sq Mean Sq F value Pr(>F)
asd_mm 1 0.8630 0.86302 22.278 3.280e-06 ***
petal_color 1 2.5753 2.57530 66.480 4.738e-15 ***
petal_area_mm 1 0.9945 0.99453 25.673 6.223e-07 ***
location 3 1.9435 0.64783 16.723 3.016e-10 ***
Residuals 395 15.3016 0.03874
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Or not (if it is entered later)
lm(prop_hybrid ~ petal_color + petal_area_mm + location + asd_mm, data = rils ) |>
anova()Analysis of Variance Table
Response: prop_hybrid
Df Sum Sq Mean Sq F value Pr(>F)
petal_color 1 3.2439 3.2439 83.7389 < 2.2e-16 ***
petal_area_mm 1 1.1524 1.1524 29.7495 8.688e-08 ***
location 3 1.9372 0.6457 16.6691 3.238e-10 ***
asd_mm 1 0.0428 0.0428 1.1053 0.2937
Residuals 395 15.3016 0.0387
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In my view Type I sums of squares are never “right,” but there are two cases in which they are “not wrong”.
In a “balanced design” in which we experimentally ensure that there are no associations between our explanatory variables, Type I and Type II sums of squares provide the same answer.
If you only care about the p-value of one term in your model, you can enter it last and then use Type I sums of squares, as the last term in a multiple regression will have the same p-value regardless of the type of sums of squares.
In practice, I think it’s best to always use Type II sums of squares for additive multiple regression models as there is no case where Type I is better, just some cases where it doesn’t really matter.
From this analysis we now see that petal area, petal color, and location are significantly associated with proportion of hybrid seeds after accounting for all other variables in the model. But for categorical variables like location, we do not know which classes differ from one another. As we saw in our introduction to the ANOVA, we answer this question with a post-hoc test.
For multiple regression, a post-hoc test asks which values of a categorical variable differ from each other after adjusting for the other terms in the model. We can do this in R with the contrast() function in the emmeans package:
library(emmeans)
location_means <- emmeans(full_model, ~ location)
emmeans(full_model, ~ location) |>
contrast(method = "pairwise", adjust = "tukey") contrast estimate SE df t.ratio p.value
GC - LB -0.0731 0.0290 395 -2.517 0.0589
GC - SR -0.0175 0.0280 395 -0.626 0.9236
GC - US 0.1143 0.0276 395 4.144 0.0002
LB - SR 0.0555 0.0282 395 1.969 0.2016
LB - US 0.1874 0.0277 395 6.758 <0.0001
SR - US 0.1319 0.0267 395 4.937 <0.0001
Results are averaged over the levels of: petal_color
P value adjustment: tukey method for comparing a family of 4 estimates From this post-hoc test, we conclude that, after adjusting for the other explanatory variables in the model, plants at location US (Upper Sawmill) have a significantly smaller proportion of hybrid seeds than plants at each other location. No other pairs of locations significantly differ from one another.
xxxx https://019b2da8-edfb-a262-61be-7973c056d9ae.share.connect.posit.cloud/intro.html