Confidence intervals for regression parameters: Bayesian vs. classical

Question

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),
                   n.chains = 4,n.adapt = 100)
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?

Answer 1

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.

Answer 2

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

Answer 3

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)

Answer 4

Given that simple linear regression is analytically identical between classical and Bayesian analysis with Jeffrey's prior, both of which are analytic, it seems a bit odd to resort to a numerical method such as MCMC to do the Bayesian analysis. MCMC is just a numerical integration tool, which allows Bayesian methods to be used in more complicated problems which are analytically intractable, just the same as Newton-Rhapson or Fisher Scoring are numerical methods for solving classical problems which are intractable.

The posterior distribution p(b|y) using the Jeffrey's prior p(a,b,s) proportional to 1/s (where s is the standard deviation of the error) is a student t distribution with location b_ols, scale se_b_ols ("ols" for "ordinary least squares" estimate), and n-2 degrees of freedom. But the sampling distribution of b_ols is also a student t with location b, scale se_b_ols, and n-2 degrees of freedom. Thus they are identical except that b and b_ols have been swapped, so when it comes to creating the interval, the "est +- bound" of the confidence interval gets reversed to a "est -+ bound" in the credible interval.

So the confidence interval and credible interval are analytically identical, and it matters not which method is used (provided there is no additional prior information) - so take the method which is computationally cheaper (e.g. the one with fewer matrix inversions). What your result with MCMC shows is that the particular approximation used with MCMC gives a credible interval which is too wide compared to the exact analytic credible interval. This is probably a good thing (although we would want the approximation to be better) that the approximate Bayesian solution appears more conservative than the exact Bayesian solution.