When using a linear model, we can generate a prediction interval without much effort. By way of example, we can use workboots to approximate a linear model’s prediction interval. 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_y_log10() +
scale_x_log10()
We can use a linear model to predict the log transform of Sale_Price based on the log transform of First_Flr_SF and 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)
mod <- lm(Sale_Price ~ First_Flr_SF, data = ames_train)
# predict on new data with a prediction interval
ames_preds <-
mod %>%
predict(ames_test, interval = "predict") %>%
as_tibble()
ames_preds %>%
# re-scale 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_y_log10() +
scale_x_log10()
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(tidymodels)
library(workboots)
# setup a workflow with a linear model
ames_wf <-
workflow() %>%
add_recipe(recipe(Sale_Price ~ First_Flr_SF, data = ames_train)) %>%
add_model(linear_reg())
# generate bootstrap predictions on ames_test
set.seed(713)
ames_preds_boot <-
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_preds_boot %>%
summarise_predictions() %>%
bind_cols(ames_preds) %>%
bind_cols(ames_test) %>%
mutate(across(c(.pred_lower:Sale_Price), ~10^.x)) %>%
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) +
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) +
scale_x_log10() +
scale_y_log10()
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 conf:
ames_preds_boot %>%
# generate 95% prediction interval
summarise_predictions(conf = 0.95) %>%
rename(.pred_lower_95 = .pred_lower,
.pred_upper_95 = .pred_upper) %>%
select(-.pred) %>%
# generate 80% prediction interval
summarise_predictions(conf = 0.80) %>%
rename(.pred_lower_80 = .pred_lower,
.pred_upper_80 = .pred_upper) %>%
bind_cols(ames_test) %>%
mutate(across(c(.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() +
scale_y_log10()