16
$\begingroup$

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?

$\endgroup$

4 Answers 4

11
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ @Ringold, what worked? The prior on sigma or the update? Or both? Have you tested them separately? $\endgroup$
    – Tomas
    Commented Jun 7, 2012 at 14:16
  • $\begingroup$ Although frankly that might be the source of the difference since it's in the original code... $\endgroup$ Commented Jun 22, 2012 at 14:25
6
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – Ringold
    Commented Jan 3, 2011 at 12:20
  • $\begingroup$ 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? $\endgroup$
    – Ringold
    Commented Jan 3, 2011 at 15:28
  • $\begingroup$ 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. $\endgroup$
    – suncoolsu
    Commented Jan 3, 2011 at 16:17
  • $\begingroup$ 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. $\endgroup$
    – user88
    Commented Jan 3, 2011 at 20:18
  • $\begingroup$ Thanks, I have registered. And many thanks to those who answered this question. I will try the gamma prior. $\endgroup$
    – Ringold
    Commented Jan 3, 2011 at 20:37
6
$\begingroup$

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)
$\endgroup$
4
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Not so odd really. One reason for using a numerical method to find an answer to a problem that can be solved analytically is to ensure that one is using the software correctly. $\endgroup$
    – Ringold
    Commented Jan 23, 2011 at 14:49
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jun 5, 2012 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.