# Estimating Linear Intervals

When using a linear model, we can generate prediction and confidence intervals without much effort. By way of example, we can use workboots to approximate linear model intervals. Let’s start by building a baseline model. In this example, we’ll predict a home’s sale price based on the first floor’s square footage with data from the Ames housing dataset.

library(tidymodels)

# setup our data
data("ames")
ames_mod <- ames %>% select(First_Flr_SF, Sale_Price)

# baseline plot
ames_mod %>%
ggplot(aes(x = First_Flr_SF, y = Sale_Price)) +
geom_point(alpha = 0.25) +
scale_x_log10(labels = scales::comma_format()) +
scale_y_log10(labels = scales::label_number(scale_cut = scales::cut_short_scale()))

We can use a linear model to predict the log transform of Sale_Price based on the log transform of First_Flr_SFand plot our predictions against a holdout set with a prediction interval.

# log transform
ames_mod <-
ames_mod %>%
mutate(across(everything(), log10))

# split into train/test data
set.seed(918)
ames_split <- initial_split(ames_mod)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)

# train a linear model
set.seed(314)
ames_lm <- lm(Sale_Price ~ First_Flr_SF, data = ames_train)

# predict on new data with a prediction interval
ames_lm_pred_int <-
ames_lm %>%
predict(ames_test, interval = "predict") %>%
as_tibble()

ames_lm_pred_int %>%

# rescale predictions to match the original dataset's scale
bind_cols(ames_test) %>%
mutate(across(everything(), ~10^.x)) %>%

# plot!
ggplot(aes(x = First_Flr_SF)) +
geom_point(aes(y = Sale_Price),
alpha = 0.25) +
geom_line(aes(y = fit),
size = 1) +
geom_ribbon(aes(ymin = lwr,
ymax = upr),
alpha = 0.25) +
scale_x_log10(labels = scales::comma_format()) +
scale_y_log10(labels = scales::label_number(scale_cut = scales::cut_short_scale()))

We can use workboots to approximate the linear model’s prediction interval by passing a workflow built on a linear model to predict_boots().

library(workboots)

# setup a workflow with a linear model
ames_wf <-
workflow() %>%
add_recipe(recipe(Sale_Price ~ First_Flr_SF, data = ames_train)) %>%

# generate bootstrap predictions on ames test
set.seed(713)
ames_boot_pred_int <-
ames_wf %>%
predict_boots(
n = 2000,
training_data = ames_train,
new_data = ames_test
)

By overlaying the intervals on top of one another, we can see that the prediction interval generated by predict_boots() (in blue) is a good approximation of the theoretical interval from lm().

ames_boot_pred_int %>%
summarise_predictions() %>%

# rescale predictions to match original dataset's scale
bind_cols(ames_lm_pred_int) %>%
bind_cols(ames_test) %>%
mutate(across(.pred:Sale_Price, ~10^.x)) %>%

# plot!
ggplot(aes(x = First_Flr_SF)) +
geom_point(aes(y = Sale_Price),
alpha = 0.25) +
scale_x_log10(labels = scales::comma_format()) +
scale_y_log10(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +

# add prediction interval created by lm()
geom_line(aes(y = fit),
size = 1) +
geom_ribbon(aes(ymin = lwr,
ymax = upr),
alpha = 0.25) +

# add prediction interval created by workboots
geom_point(aes(y = .pred),
color = "blue",
alpha = 0.25) +
geom_errorbar(aes(ymin = .pred_lower,
ymax = .pred_upper),
color = "blue",
alpha = 0.25,
width = 0.0125)

Both lm() and summarise_predictions() use a 95% prediction interval by default but we can generate other intervals by passing different values to the parameter interval_width:

ames_boot_pred_int %>%

# generate a 95% prediction interval
summarise_predictions(interval_width = 0.95) %>%
rename(.pred_lower_95 = .pred_lower,
.pred_upper_95 = .pred_upper) %>%
select(-.pred) %>%

# generate 80% prediction interval
summarise_predictions(interval_width = 0.80) %>%
rename(.pred_lower_80 = .pred_lower,
.pred_upper_80 = .pred_upper) %>%

# rescale predictions to match original dataset's scale
bind_cols(ames_test) %>%
mutate(across(.pred_lower_95:Sale_Price, ~10^.x)) %>%

# plot!
ggplot(aes(x = First_Flr_SF)) +
geom_point(aes(y = Sale_Price),
alpha = 0.25) +
geom_line(aes(y = .pred),
size = 1,
color = "blue") +
geom_ribbon(aes(ymin = .pred_lower_95,
ymax = .pred_upper_95),
alpha = 0.25,
fill = "blue") +
geom_ribbon(aes(ymin = .pred_lower_80,
ymax = .pred_upper_80),
alpha = 0.25,
fill = "blue") +
scale_x_log10(labels = scales::comma_format()) +
scale_y_log10(labels = scales::label_number(scale_cut = scales::cut_short_scale()))

Alternatively, we can estimate the confidence interval around each prediction by passing the argument "confidence" to the interval parameter of predict_boots().

# generate linear model confidence interval for reference
ames_lm_conf_int <-
ames_lm %>%
predict(ames_test, interval = "confidence") %>%
as_tibble()

# generate bootstrap predictions on test set
set.seed(867)
ames_boot_conf_int <-
ames_wf %>%
predict_boots(
n = 2000,
training_data = ames_train,
new_data = ames_test,
interval = "confidence"
)

Again, by overlaying the intervals on the same plot, we can see that the confidence interval generated by predict_boots() is a good approximation of the theoretical interval.

ames_boot_conf_int %>%
summarise_predictions() %>%

# rescale predictions to match original dataset's scale
bind_cols(ames_lm_conf_int) %>%
bind_cols(ames_test) %>%
mutate(across(.pred:Sale_Price, ~10^.x)) %>%

# plot!
ggplot(aes(x = First_Flr_SF)) +
geom_point(aes(y = Sale_Price),
alpha = 0.25) +
scale_x_log10(labels = scales::comma_format()) +
scale_y_log10(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +

# add prediction interval created by lm()
geom_line(aes(y = fit),
size = 1) +
geom_ribbon(aes(ymin = lwr,
ymax = upr),
alpha = 0.25) +

# add prediction interval created by workboots
geom_point(aes(y = .pred),
color = "blue",
alpha = 0.25) +
geom_ribbon(aes(ymin = .pred_lower,
ymax = .pred_upper),
fill = "blue",
alpha = 0.25)