check-workboots-heteroskedasticity.R

library(tidyverse)
library(tidymodels)
library(workboots)

data <- tibble(x = abs(rnorm(1000)), 
               y = x + rnorm(1000, sd = 0.7) * x / (max(x))
)

rec <- recipe(y ~ ., data = data)

mod <- parsnip::linear_reg() %>% 
  set_engine("lm") %>% 
  set_mode("regression")

workflow <- workflows::workflow() %>% 
  add_recipe(rec) %>% 
  add_model(mod) 

set.seed(345)
pred_int <-
  workflow %>%
  predict_boots(
    n = 500,
    training_data = data,
    new_data = data
  )

# summarise predictions with a 95% prediction interval
bind_cols(data, 
          pred_int %>% summarise_predictions()
) %>% 
  ggplot(aes(x = x, y = y))+
  geom_line(aes(y = .pred, colour = ".pred"))+
  geom_ribbon(aes(ymin = .pred_lower, ymax = .pred_upper), alpha = 0.2)+
  geom_point()+
  labs(title = "workboots also does not capture heteroskedasticity in intervals")