Simple Linear Regression: ANOVA - Cats
Chapter 2
STA 4210
Looking at the cats data:
set.seed(4210)
cats |>
slice_sample(n = 10)## body heart
## 1 2.1 7.6
## 2 2.9 10.1
## 3 2.5 11.0
## 4 3.6 15.0
## 5 2.9 11.8
## 6 3.2 12.3
## 7 2.2 8.7
## 8 2.5 8.8
## 9 2.2 11.0
## 10 3.1 12.1Fitting the linear regression model
We’ll start by fitting a linear model for heart weight by body weight
using lm(hearty ~ body, cats)
cats_lm <- lm(heart ~ body, cats)
summary(cats_lm)##
## Call:
## lm(formula = heart ~ body, data = cats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5694 -0.9634 -0.0921 1.0426 5.1238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3567 0.6923 -0.515 0.607
## body 4.0341 0.2503 16.119 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.452 on 142 degrees of freedom
## Multiple R-squared: 0.6466, Adjusted R-squared: 0.6441
## F-statistic: 259.8 on 1 and 142 DF, p-value: < 2.2e-16Adding the predicted heart weight to the data
While there are multiple ways to adding \(\widehat{heart}\) to the
cats data set, the code chunk below uses
predict(model, data) to do it:
cats <-
cats |>
mutate(heart_hat = predict(cats_lm, newdata = cats))
cats |>
slice_sample(n = 10)## body heart heart_hat
## 1 3.0 13.0 11.745526
## 2 2.8 10.2 10.938713
## 3 2.3 9.6 8.921682
## 4 2.1 8.1 8.114869
## 5 2.9 10.1 11.342119
## 6 2.0 7.4 7.711463
## 7 2.2 11.0 8.518276
## 8 2.3 9.0 8.921682
## 9 3.4 14.4 13.359151
## 10 2.8 13.3 10.938713ANOVA: Sums of Squares
To perform an ANOVA test, we need three sums of squares from the regression model:
- \(SSTO\): The total sums of squares - The difference between \(Y\) and the average, \(\bar{Y}\)
\[SSTO = \sum_i^n(Y_i - \bar{Y})^2\]
SSTO <- sum((cats$heart - mean(cats$heart))^2)
round(SSTO, 2)## [1] 847.63- \(SSR\): The regression sums of squares - The difference between \(\hat{Y}_i\) and \(\bar{Y}\)
\[SSR = \sum_i^n(\hat{Y}_i - \bar{Y})^2\]
SSR <- sum((cats$heart_hat - mean(cats$heart))^2)
round(SSR, 2)## [1] 548.09- \(SSE\): The error sums of squares - The difference between \(Y_i\) and \(\hat{Y}_i\)
- AKA, the sum of the residuals squared, which is the quantity we minimized to estimate \(b_0\) and \(b_1\)
\[SSE = \sum_i^n(\hat{Y}_i - \bar{Y})^2\]
SSR <- sum((cats$heart - cats$heart_hat)^2)
round(SSR, 2)## [1] 299.53I’ll put them all together into one data frame using
summarize(...)
SS_table <-
cats |>
summarize(
Reg = sum((heart_hat - mean(heart))^2),
Error = sum((heart - heart_hat)^2),
Total = sum((heart - mean(heart))^2)
) |>
# Pivoting the data frame so it looks like an ANOVA table
pivot_longer(
cols = everything(),
names_to = 'source',
values_to = 'SS'
)
SS_table## # A tibble: 3 × 2
## source SS
## <chr> <dbl>
## 1 Reg 548.
## 2 Error 300.
## 3 Total 848.Full ANOVA table
Let’s add the degrees of freedom and MS columns to SS_table to make it into an ANOVA table:
ANOVA_table <-
SS_table |>
mutate( # SSR df # SSE df # SSTO df
df = c( 1, nrow(cats) - 2, nrow(cats) - 1),
MS = SS / df
)
ANOVA_table## # A tibble: 3 × 4
## source SS df MS
## <chr> <dbl> <dbl> <dbl>
## 1 Reg 548. 1 548.
## 2 Error 300. 142 2.11
## 3 Total 848. 143 5.93Note: The MS value in the Error column is our MSE we found the last few R examples!
Additional Note: While our hand made ANOVA table has \(MST\), we don’t typically calculate it since we never use it.
Conducting the ANOVA test
Once we have the complete ANOVA table, we can perform our ANOVA test.
The null and alternative hypotheses are:
\[H_0: \text{None of the predictors are linearly related to the response}\]
\[H_a: \text{At least one of the predictors is linearly related to the response}\]
or with context to our problem
\[H_0: \text{Body weight is linearly not related to heart weight in cats} \\ H_a: \text{Body weight is linearly related to heart weight in cats}\]
The test statistic is the ratio of the Mean Squared Regression vs Mean Squared Error:
\[F = MSR/MSE \sim F(1, n-2)\]
F_stat <- ANOVA_table$MS[1] / ANOVA_table$MS[2]
c(
'F stat' = F_stat,
'num df' = ANOVA_table$df[1],
'den df' = ANOVA_table$df[2]
) |>
round(2)## F stat num df den df
## 259.83 1.00 142.00From there, we can find our p-value using the stated \(F\)-distribution above:
\[P(F(1, n - 2) > F)\]
pf(F_stat, df1 = ANOVA_table$df[1], df2 = ANOVA_table$df[2], lower.tail = F)## [1] 6.969045e-34We have very, very strong evidence that body weight is linearly related to heart weight!
Performing the ANOVA test using built in functions
Since ANOVA tests with regression are very, very common, there are plenty of functions to do what we just did above without us having to manually calculate every \(SS\), df, and \(MS\) value.
Method 1: summary(model)
The summary function that we saw initially will perform our test for us:
summary(cats_lm)##
## Call:
## lm(formula = heart ~ body, data = cats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5694 -0.9634 -0.0921 1.0426 5.1238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3567 0.6923 -0.515 0.607
## body 4.0341 0.2503 16.119 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.452 on 142 degrees of freedom
## Multiple R-squared: 0.6466, Adjusted R-squared: 0.6441
## F-statistic: 259.8 on 1 and 142 DF, p-value: < 2.2e-16Where is it located in all that output?
If you look at the last line, you’ll see the test statistic, degrees of freedom, and p-value. All of them agree with what we found!
Method 2: glance(model) in broom
package
In a previous example, we used the tidy() function from
the broom package to find the p-value for the test of the
slope alone.
The same package has a function called glance() that
will return several fit statistics for the entire model in a one-row
data frame:
broom::glance(cats_lm)## # A tibble: 1 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.647 0.644 1.45 260. 6.97e-34 1 -257. 520. 529.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>That’s a lot of fit statistics! Let’s just look at the relevant ones:
broom::glance(cats_lm) |>
dplyr::select(statistic, df, df.residual, p.value)## # A tibble: 1 × 4
## statistic df df.residual p.value
## <dbl> <dbl> <int> <dbl>
## 1 260. 1 142 6.97e-34The above also agrees with what we found!
Method 3: anova(model)
The anova() function will perform an ANOVA test, but
instead of testing the entire model, it will be a test for each
predictor in the model. Which is more similar to the \(t\)-test that we did earlier, at least for
now…
anova(cats_lm)## Analysis of Variance Table
##
## Response: heart
## Df Sum Sq Mean Sq F value Pr(>F)
## body 1 548.09 548.09 259.83 < 2.2e-16 ***
## Residuals 142 299.53 2.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So what is the advantage of using anova() over
summary() or tidy(), since both are testing
each explanatory variable individually?
Once we start adding categorical explanatory variables,
anova() will be important, but for now, it’s not
beneficial.
SLR: F-test vs t-test
We get the same result conducting an F-test or a t-test. That’s not by coincidence!
The main four distributions (\(z\), \(t\), \(\chi^2\), \(F\)) are all connected!
It’s relevant here because the slope test is:
\[t = \frac{b_1}{\sqrt{MSE/S_{XX}}} \sim t(n-2)\]
If you take a \(t\)-distributed random variable and square it, you get an \(F\)-statistic with 1 numerator df and the denominator df are the same as the \(t\)-distributions!
\[t^2 \sim F(1; n-2)\]
If we take the test stat for the slope and square it, we should get the same as the \(F\) test stat:
broom::tidy(cats_lm) |>
slice(2) |>
dplyr::select(t_stat = statistic) |>
mutate(t_stat = t_stat^2)## # A tibble: 1 × 1
## t_stat
## <dbl>
## 1 260.which is the same as the \(F\) stat we’ve seen many times in this example!
Note: The F-test version is only equivalent to a two-tailed test. While we can do a left or right-tail test using the \(t\)-distribution, we can’t for the \(F\) version :(