Fake Data Loop

n_fake <- 1000
cover_68 <- rep(NA, n_fake)
cover_95 <- rep(NA, n_fake)
for (s in 1:n_fake){
  y <- a + b*x + rnorm(n,0,sigma)
  fake <- tibble(x,y)
  fit <- stan_glm(y ~ x, data = fake, refresh = 0)
  b_hat <- coef(fit)["x"]
  b_se <- se(fit)["x"]
  # Is true value in the corresponding interval?
  cover_68[s] <- abs(b - b_hat) < b_se
  cover_95[s] <- abs(b - b_hat) < 2*b_se
}

cat(paste("68% coverage: ", mean(cover_68), "\n"))
cat(paste("95% coverage: ", mean(cover_95), "\n"))

68% coverage:  0.667 
95% coverage:  0.941

Not far from nominal values, which gives us confidence in the fitting procedure!