Confidence intervals for regression parameters: Bayesian vs. classical

Given two arrays x and y, both of length n, I fit a model y = a + b*x and want to calculate a 95% confidence interval for the slope. This is (b - delta, b + delta) where b is found in the usual way and

delta = qt(0.975,df=n-2)*se.slope


and se.slope is the standard error in the slope. One way to get the standard error of the slope from R is summary(lm(y~x))$coef[2,2]. Now suppose I write the likelihood of the slope given x and y, multiply this by a "flat" prior and use a MCMC technique to draw a sample m from the posterior distribution. Define lims = quantile(m,c(0.025,0.975))  My question: is (lims[[2]]-lims[[1]])/2 approximately equal to delta as defined above? Addendum Below is a simple JAGS model where these two seem to be different. model { for (i in 1:N) { y[i] ~ dnorm(mu[i], tau) mu[i] <- a + b * x[i] } a ~ dnorm(0, .00001) b ~ dnorm(0, .00001) tau <- pow(sigma, -2) sigma ~ dunif(0, 100) }  I run the following in R: N <- 10 x <- 1:10 y <- c(30.5,40.6,20.5,59.1,52.5, 96.0,121.4,78.9,112.1,128.4) lin <- lm(y~x) #Calculate delta for a 95% confidence interval on the slope delta.lm <- qt(0.975,df=N-2)*summary(lin)$coef[2,2]

library('rjags')
jags <- jags.model('example.bug', data = list('x' = x,'y' = y,'N' = N),
update(jags, 1000)
params <- jags.samples(jags,c('a', 'b', 'sigma'),7500)
lims <- quantile(params$b,c(0.025,0.975)) delta.bayes <- (lims[[2]]-lims[[1]])/2 cat("Classical confidence region: +/-",round(delta.lm, digits=4),"\n") cat("Bayesian confidence region: +/-",round(delta.bayes,digits=4),"\n")  And get: Classical confidence region: +/- 4.6939 Bayesian confidence region: +/- 5.1605 Rerunning this multiple times, the Bayesian confidence region is consistently wider than the classical one. So is this due to the priors I've chosen? - add comment 4 Answers The 'problem' is in the prior on sigma. Try a less informative setting tau ~ dgamma(1.0E-3,1.0E-3) sigma <- pow(tau, -1/2)  in your jags file. Then update a bunch update(10000)  grab the parameters, and summarise your quantity of interest. It should line up reasonably well with the classic version. Clarification: The updating is just to make sure you get where you're going whatever choice of prior you decide on, although chains for models like this one with diffuse priors and random starting values do take longer to converge. In real problems you'd check convergence before summarising anything, but convergence is not the main issue in your example I don't think. - Thanks, this worked. – Ringold Jan 3 '11 at 23:37 @Ringold, what worked? The prior on sigma or the update? Or both? Have you tested them separately? – Tomas Jun 7 '12 at 14:16 should be sigma <- pow(tau, -1/2) or sigma <- 1/sqrt(tau) – Tomas Jun 20 '12 at 17:11 @Tomas, quite right. Typo fixed. – conjugateprior Jun 22 '12 at 14:09 Although frankly that might be the source of the difference since it's in the original code... – conjugateprior Jun 22 '12 at 14:25 add comment If you sample from the posterior of b | y and calculate lims (as you define) it should be same as (b - delta, b + delta). Specifically, if you calculate the posterior distribution of b | y under a flat prior, it is same as the classical sampling distribution of b. For more details refer to: Gelman et al. (2003). Bayesian Data Analysis. CRC Press. Section 3.6 Edit: Ringold, the behavior observed by you is consistent with the Bayesian idea. The Bayesian Credible Interval (CI) is generally wider than the classical ones. And the reason is, as you correctly guessed, the hyperpriors taken into account the variability because of the unknown parameters. For simple scenarios like these (NOT IN GENERAL): Baysian CI > Empirical Bayesian CI > Classical CI ; > == wider - I added some code using JAGS where I seem to be getting a different answer. Where is my mistake? Is this happening because of the priors? – Ringold Jan 3 '11 at 12:20 Now I'm confused. First you said that the posterior distribution of b|y under a flat prior is the same as the classical sampling distribution of b. Then you said the Bayesian CI is wider than the classical one. How could it be wider if the distributions are the same? – Ringold Jan 3 '11 at 15:28 Sorry - I should have told what @CP suggested in his comments. Theoretically, b|y under a flat prior and classical CI are the same, but you can't do that practically in JAGS unless you use a very very diffuse prior like CP suggested and use a lot of MCMC iterations. – suncoolsu Jan 3 '11 at 16:17 I have merged your accounts so you can edit your questions and add comments. Yet, please register your account by clicking here: stats.stackexchange.com/users/login ; you can use your Gmail OpenID to do it in a few seconds and you'd not loose your account here any more. – mbq Jan 3 '11 at 20:18 Thanks, I have registered. And many thanks to those who answered this question. I will try the gamma prior. – Ringold Jan 3 '11 at 20:37 add comment For linear Gaussian models it is better to use the bayesm package. It implements the semi-conjugate family of priors, and the Jeffreys prior is a limit case of this family. See my example below. These are classical simulations, there's no need to use MCMC. I do not remember whether the credibility intervals about the regression parameters are exactly the same as the usual least squares confidence intervals, but in any case they are very close. > # required package > library(bayesm) > # data > age <- c(35,45,55,65,75) > tension <- c(114,124,143,158,166) > y <- tension > # model matrix > X <- model.matrix(tension~age) > # prior parameters > Theta0 <- c(0,0) > A0 <- 0.0001*diag(2) > nu0 <- 0 > sigam0sq <- 0 > # number of simulations > n.sims <- 5000 > # run posterior simulations > Data <- list(y=y,X=X) > Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq) > Mcmc <- list(R=n.sims) > bayesian.reg <- runireg(Data, Prior, Mcmc) > beta.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw > sigmasq.sims <- bayesian.reg$sigmasqdraw
> apply(beta.sims, 1, quantile, probs = c(0.025, 0.975))
[,1] [,2]
2.5% 53.33948 1.170794
97.5% 77.23371 1.585798
> # to be compared with:
> frequentist.reg <- lm(tension~age)

-
Another reason for using simulations: you also get posterior simulations for any function $f(\beta_0, \beta_1, \ldots, \beta_p, \sigma)$ of the parameters. For instance I use this possibility to get a credibility interval about the probability $\Pr(Y > 10 \mid x)$ for some value $x$ of the covariate. I don't know how to get a frequentist confidence interval about this probability. –  Stéphane Laurent Jun 5 '12 at 15:01