Fit data to model (fit to constant term)

set.seed(264)
fit <- stan_glm(y ~ 1, data = newcomb, refresh = 0)

Simulate from the predictive distribution.

set.seed(970)

sims <- as_tibble(fit)

n_sims <- nrow(sims)
n_newcomb <- nrow(newcomb)

y_rep_tidy <- 
  sims %>% 
  mutate(rep = row_number()) %>% 
  group_by(rep) %>% 
  reframe(y = rnorm(n_newcomb, mean = `(Intercept)`, sd = sigma))   

y_rep_tidy is a tidy tibble with 4000 * 66 rows.

Visual comparison of actual and replicated datasets

Plot histograms for 20 sample replicates.