Inference for Slope - Cats example
Chapter 2
STA 4210
Loading the data from the MASS package and cleaning the
names
Checking the data
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.1Scatter plot for body vs heart:
We should start with a plot of body vs heart weight:
gg_cats <-
ggplot(
data = cats,
mapping = aes(
x = body,
y = heart
)
) +
geom_point() +
theme_bw() +
labs(
title = 'Cats: Body weight vs heart weight',
subtitle = 'LOESS line added',
x = 'Body',
y = 'Heart'
) +
# Adding the units to the axes
scale_x_continuous(
labels = scales::label_number(suffix = ' kg')
) +
scale_y_continuous(
labels = scales::label_number(suffix = ' g')
)
gg_cats +
geom_smooth(
method = 'loess',
formula = y ~ x,
se = F
)
Fitting the linear model:
While we’ve been calculating the slope ‘by hand’, we typically use a
function to do it for us, like lm(y ~ x, data = ...) and we
can get the summary stats from summary(model)
cats_lm <- lm(heart ~ body, data = 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-16The broom package has some useful functions that are
popular as they return objects that are easier to work with
tidy(model)returns the coefficient table as a data frameglance(model)returns a 1 row data frame with the fit statistics of the model- We’ll get to that later
augment_columns(model, data)will add the predicted response, residuals, and other useful stats to the data frame
For now, we’ll just work with tidy()
cats_lm_table <- broom::tidy(cats_lm)
cats_lm_table## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.357 0.692 -0.515 6.07e- 1
## 2 body 4.03 0.250 16.1 6.97e-34# Finding the MSE:
sum(cats_lm$residuals^2) / (nrow(cats) - 2)## [1] 2.109388# Finding S_XX
sum((cats$body - mean(cats$body))^2)## [1] 33.67972Hypothesis test for the slope
The first step of any hypothesis test is to state the null and alternative hypotheses:
\[H_0: \beta_1 = 0 \\ H_a: \beta_1 \ne 0\]
We’re interested in seeing if there is any type of linear association between body and heart weight.
Next, we’d check the conditions, but that’s for a later chapter.
Next step: Test statistic
By hand
We’ll start by performing the hypothesis test for the slope manually, then use some of the above methods to find the test statistic and p-value
Test stat:
The test statistic for the slope is:
\[t = \frac{b_1-\beta_{1,0}}{\sqrt{\frac{MSE}{S_{XX}}}}\]
# Start by finding n, MSE, S_XX, and slope
n <- nrow(cats)
MSE_cats <- sum(cats_lm$residuals^2) / (n - 2)
S_XX <- sum((cats$body - mean(cats$body))^2)
slope_cats <- cats_lm$coefficients[2]
# Standard error next
slope_SE <- sqrt(MSE_cats/S_XX)
slope_ts <- slope_cats / slope_SE
slope_ts## body
## 16.11939We have a huge test statistic, which would imply we will have what is called in scientific terms, ‘teeny tiny’.
P-value
The test statistic follows a t-distribution with \(df = n - 2\)
\[t = \frac{b_1-\beta_{1,0}}{\sqrt{\frac{MSE}{S_{XX}}}} \sim t(n-2)\]
slope_p_val <- 2 * pt(abs(slope_ts), df = n-2, lower.tail = F)
slope_p_val## body
## 6.969045e-34So we have a p-value that has 33 zeros between the decimal and 6.
0.0000000000000000000000000000000006969
So we have very, very strong evidence that there is a linear association between body and heart weight, which is not surprising given the scatter plot we created in the beginning. But it is nice to have the methods reaffirm our subjective conclusions.
95% confidence interval ‘by hand’
Given the conclusion isn’t all that useful (there is a linear association), we want to follow it up with a confidence interval.
\[b_1 \pm t(1 - \alpha/2; n - 2) \sqrt{\frac{MSE}{S_{XX}}}\]
cats_slope_CI <-
c('lower' = slope_cats - qt(0.975, df = n - 2) * slope_SE,
'upper' = slope_cats + qt(0.975, df = n - 2) * slope_SE)
cats_slope_CI |>
round(2)## lower.body upper.body
## 3.54 4.53Slope inference using built-in functions
There are several functions we can use to calculate the test statistic, p-value, and confidence interval
summary(model)
Up first, the base function
summary(model).
cats_lm_sum <- summary(cats_lm)
cats_lm_sum##
## 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-16There’s a lot of output from the summary() function. If
you want some smaller pieces, we can extract them using the
$ operator. What can we put after the $?
names(object) will tell you all the names you can put
after the $
names(cats_lm_sum)## [1] "call" "terms" "residuals" "coefficients"
## [5] "aliased" "sigma" "df" "r.squared"
## [9] "adj.r.squared" "fstatistic" "cov.unscaled"# Let's look at coefficients
cats_lm_sum$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3566624 0.6922770 -0.5152019 6.072131e-01
## body 4.0340627 0.2502615 16.1193908 6.969045e-34The code above returns just the coefficients table, which is want we want if we’re just interested in the slope.
And hooray, it returns the same standard error, test statistic, and p-value that we found by hand.
What about a confidence interval?
There’s not a way to get it directly from the summary()
object. But there is a base function called
confint(model) that will create a confidence interval for
each parameter in the model
confint(cats_lm, level = 0.95) |> round(2)## 2.5 % 97.5 %
## (Intercept) -1.73 1.01
## body 3.54 4.53If you just want the interval for a specific parameter(s), you can
give parm the name of the row:
confint(cats_lm, parm = 'body', level = 0.95) |> round(2)## 2.5 % 97.5 %
## body 3.54 4.53Important: The confidence interval created by
confint() is different than the standard formula of \(stat \pm CV * SE\), but for larger samples
or if the data are approximately Normal, the two will be very, very
similar.
So what can we do if we want the test statistic, p-value, and confidence interval all in one place?
broom::tidy()
We can use the tidy(model) function in the
broom package, which returns just the table of the model
terms:
broom::tidy(cats_lm)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.357 0.692 -0.515 6.07e- 1
## 2 body 4.03 0.250 16.1 6.97e-34Test statistic and p-value, but no confidence interval :(
It doesn’t include the CI by default, but if you specify
conf.int = T, then it will create a 95% confidence
interval. If you want a different confidence level, you can include
conf.level = ... to pick a different one (needs to be
between 0 and 1, not a percentage).
broom::tidy(
cats_lm,
conf.int = T,
conf.level = 0.95
) |>
mutate(
across(
.cols = where(is.numeric),
.fns = ~ round(., 2)
)
)## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.36 0.69 -0.52 0.61 -1.73 1.01
## 2 body 4.03 0.25 16.1 0 3.54 4.53Same that we got the previous two times!
Personally, I like the broom functions for linear models
and will use them over summary() most of the time because
it’s easy to alter them using the dplyr functions like
mutate()