# Help me understand $p$-values in Bayesian glm

I am trying to run a Bayesian logit on the data here. I am using bayesglm() in the arm package in R. The coding is straightforward enough:

df = read.csv("http://dl.dropbox.com/u/1791181/bayesglm.csv", header=T)
library(arm)
model = bayesglm(PASS ~ SEX + HIGH, family=binomial(link="logit"), data=df)


summary(model) gives the following output:

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.10381    0.10240   1.014    0.311
SEXMale      0.02408    0.09363   0.257    0.797
HIGH        -0.27503    0.03562  -7.721 1.15e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2658.2  on 1999  degrees of freedom
Residual deviance: 2594.3  on 2000  degrees of freedom
AIC: 2600.3


Please walk me through this. I understand that this code uses a very weak prior (since I am not specifying the prior means) so the output is going to be practically the same if I used glm() instead of bayesglm(). But the output should still be in the Bayesian spirit, right? What are the $p$-values and $z$-values here? Aren't these frequentist inference tools? Are they interpreted differently here?

-
This is a comment and not an answer, but this is what would make some sense to me. You get estimates which are probably the values for which the posterior distribution is maximized. It might also be possible that they are just the means of the posterior? Worth checking out if you can. But no matter what the exact details are, once you have some estimates you can test them by the usual Estimate/Std. Error -> z-score procedure which works if the posterior is close enough to a normal (it goes to a normal under some conditions that usually hold). –  Erik Jul 2 '12 at 6:29
Erik... You are correct: the coefficients are indeed the means of the posterior densities. My question is about the p- and z-values. What do they represent here? –  user3671 Jul 2 '12 at 6:32
Ok. If you have a density that is approximately normally distributed you can test whether it's mean is zero by taking the z-score = mean / standard deviation of the distribution and comparing it with the standard normal distribution. Then you look how unlikely your value or a larger value would be under the standard normal distribution -> p-value. See z-score on Wikipedia for details. –  Erik Jul 2 '12 at 7:06
Well, yes. But why bother doing that in a Bayesian setting? In Bayesian inference, the point estimate is my best guess about the random parameter, so there's no need to test it. At the most, I can include a "credible interval" which is equivalent to a frequentist "confidence interval" but whose statistical interpretation is vastly different. This the confusing part about the summary() output. The spirit is Bayesian, but the output is frequentist? –  user3671 Jul 2 '12 at 7:13
One point is that the estimate you get will be different, since you used a prior. And while the point estimate is the "best guess" if you want to show in a Bayesian way that something has an effect you would try to show that the credible interval does not contain the zero. When you approximate the posterior by a normal with the same mean and sd (asymptotically correct) then the (1-p/2) credibility interval is the largest symmetrical credibility interval containing the zero, so your answer is basically the same. The p is p-value stated above. –  Erik Jul 2 '12 at 7:32

Great question! Although there are Bayesian p-values, and one of the authors of the arm package is an advocate, what you are seeing in your output is not a Bayesian p-value. Check the class of model

class(model)
"bayesglm" "glm"      "lm"


and you can see that class bayesglm inherits from glm. Furthermore, examination of the arm package shows no specific summary method for a bayesglm object. So when you do

summary(model)


you are actually doing

summary.glm(model)


and getting frequentist interpretation of the results. If you want a more Bayesian perspective the function in arm is display()

-
+1 Excellent answer! This is the trouble with R, there are so many highly intelligent statisticians who write horrible code that leaves these kinds of landmines lying around. –  Bogdanovist Jul 2 '12 at 22:05
It seems like a deliberate choice on the designers part, rather than an oversight. –  atiretoo Jul 2 '12 at 22:23
After reading the link I agree with the intent, but in that case summary() should have been re-implemented to simply call display() rather than giving nonsense results without warning. The person who asked this question got tripped up by a code that broke the user model for R that has been established by every other object they have ever used. That's terrible programming practice. –  Bogdanovist Jul 2 '12 at 22:37
Many thanks, atiretoo. This raises another question. What's the difference between display() and summary()? It seems to me that the output from the former is just the output from the latter, less two columns, and rounded to 2 digits. It would appear so, from Gelman's post you've linked above. –  user3671 Jul 3 '12 at 12:04
Yes, and from the discussion at Andrew Gelman's blog it sounds like they'll fix this in future versions of the arm package. –  atiretoo Jul 3 '12 at 19:54