Rstanarm - for my own learning, trying to match credible interval to lm() confidence interval

Hi all,

I’m a relative beginner with regards to Bayesian inference. I’ve worked my way through Statistical Rethinking (a fantastic book), and I see it’s been translated into brms, which is great, but I’m really enjoying using rstanarm and am learning to use it as my next step. I like to tinker with things to see how they work, and rstanarm makes that so easy because it’s so fast. I can change things and re-run code in seconds. No compiling needed! Seems like rstanarm doesn’t get as much love as it should, but I’m certainly digging it.

One of the things I’m trying to do is figure out how to set everything up so that I can create credible intervals that coincide with confidence intervals. This is strictly as a learning tool, since I’m a firm believer in using informative priors if possible, and regularizing priors if not. My understanding is that getting Cred-I and Conf-I to coincide is possible for a Gaussian linear model using flat priors. Below is my first attempt. I’ve spent quite a few hours going down paths and learning a lot, but I’m stuck at this point. Any help would be appreciated.

The fundamental question is, why does the credible interval not match the confidence interval from lm()? Many more details within the code as comments.


# First, create the simple 5-point dataset with mean 0 and sd of 2

outcome <- c(-2,-2,0,2,2)
dat <- as.data.frame(outcome)

# lm() --------------------------------------------------------------------

m.lm <- lm(outcome ~ 1, data=dat)
summary(m.lm)
confint(m.lm, level = 0.95)

2/sqrt(5)
qt(p=0.025, df=4)
qt(p=0.025, df=4) * 2/sqrt(5)

# The SEM is 0.8944 for the lm() analysis. Matches my manual calculation of 0.8944272.
# The 95% confidence interval for the mean is -2.483328 2.483328

# rstanarm ----------------------------------------------------------------

# Ok, let's try this with stan_glm and flat priors

library(rstanarm)

m.stan_glm <- stan_glm(
  formula = outcome ~ 1 ,
  family = gaussian ,
  prior_intercept = NULL ,
  prior_aux = NULL ,
  chains = 4 ,
  cores = 4 ,
  iter = 2e5 ,
  seed = 42 ,
  data = dat
)

print(m.stan_glm$stanfit, digits=8)

# sd of the intercept is 1.772181, far greater than 0.8944 from lm()
# The 95% interval is also far wider than for lm()
# -3.297692 to  +3.303978
# -2.483328 to  +2.483328 for lm()
# Is this because of sigma not being a point (without uncertainty)?  There's no
# uncertainty around sigma in the lm() analysis, right?  I'll try to fix that.

# What if I use a super-tight prior at 2 for sigma? As far as I can tell from my experimentation,
# the prior of a normal for sigma cuts off the left side at the mean, even if it's in postive
# territory. It's forced to be a half-normal, so I'll start it 1 standard deviation (sd=1e-6) below 2.

m.stan_glm2 <- stan_glm(
  formula = outcome ~ 1 ,
  family = gaussian ,
  prior_intercept = NULL ,
  prior_aux = normal(2-1e-6, 1e-6) ,
  chains = 4 ,
  cores = 4 ,
  iter = 2e5 ,
  seed = 42 ,
  data = dat
)

prior_summary(m.stan_glm2)
print(m.stan_glm2$stanfit, digits=8)

# sigma is now 2, and the Intercept sd is 0.892 (close to 0.8944)
# However, the width of the 95% intervals are far more _narrow_ than the lm() intervals
# -1.747336 to  +1.755980
# -2.483328 to  +2.483328 for lm()

posterior <- as.array(m.stan_glm2)
posterior_vector_intercept <- as.numeric(posterior[,1:4,1])
posterior_vector_sigma <- as.numeric(posterior[,1:4,2])
hist(posterior_vector_intercept)
hist(posterior_vector_sigma)
( posterior_intercept <- mean(posterior_vector_intercept) )
( posterior_sigma <- mean(posterior_vector_sigma) )
( posterior_intercept_sd <- sd(posterior_vector_intercept) )

xlim_min <- -5
xlim_max <- 5
hist(posterior_vector_intercept, breaks=200, freq=FALSE, xlim=c(xlim_min,xlim_max))
curve(dnorm(x, mean=posterior_intercept, sd=posterior_intercept_sd),
      from=xlim_min, to=xlim_max, col="red", add=TRUE)

# The histogram and the curve are on top of each other. So, the histogram is Gaussian.
# If the posterior_intercept_sd is 0.892, which is close to the lm() SEM of 0.8944,
# why are the uncertaity limits so different?
# -1.747    to  +1.756
# -2.483328 to  +2.483328 for lm()

min(posterior_vector_intercept)
#[1] -4.004504
max(posterior_vector_intercept)
#[1] 4.113033
# No severe outliers to mess with things.

# eof

Operating system: Windows 10
R version 4.0.2
rstanarm version 2.21.1
rstan version 2.21.2 (not used when using rstanarm?)

I suspect there are a couple of problems here:

  1. c(-2,-2,0,2,2) isn’t a very likely sample from N(0, 2). There’s much more weight in the tails than you’d expect. lm() just does a least squares fit and doesn’t worry about that.

  2. With only 5 data points the posterior distribution for a variance (or standard deviation) parameter will be heavily influenced by the prior.

See what happens if you set outcome like this:

outcome <- rnorm(1000, sd = 2)

I get

(-0.019, 0.108) from lm()
(-0.018, 0.110) from stan_glm()

Kent