4
$\begingroup$

I am trying to recreate (in R) a frequentist hypothesis testing in Bayesian from, by calculating Bayes factors of the null (H0) and alternative (H1) models.

The model is simply a simple linear regression that tries to detect a trend in global temp. data from 1995 to 2009 (here). Therefore, H0 is no trend (i.e. slope = 0), or similary, the H0 model is a linear model with only the intercept.

So I calculated the lm() of both models to arrive at negative log likelihood values that are significantly different. The p-value for the H1 lm() model is 0.0877.

I also calculated this in a Bayesian way by using MCMCpack, and I get negative log likelihood values that are super duper uber different. Log likelihood values of 13.7 and 4.3 are about a 10000 fold difference in their likelihood ratios (where >100 is considered to be "decisive").

The means and sds of the estimates are very similar, so why am I getting such different likelihood values? (particularly for the Bayesian H0 model) I feel like there is a gap in my understanding on marginal likelihoods, but I can't pinpoint the problem.

Thanks

library(MCMCpack)

## data: http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3gl.txt

head(hadcru, 2)
##  Year      1      2      3      4      5      6      7      8      9     10
## 1 1850 -0.691 -0.357 -0.816 -0.586 -0.385 -0.311 -0.237 -0.340 -0.510 -0.504
## 2 1851 -0.345 -0.394 -0.503 -0.480 -0.391 -0.264 -0.279 -0.175 -0.211 -0.123
##       11     12    Avg
## 1 -0.259 -0.318 -0.443
## 2 -0.141 -0.151 -0.288

hadcru.lm <- lm(Avg ~ 1 + Year, data = subset(hadcru, (Year <= 2009 & Year >= 1995)))
hadcru.lm.zero <- lm(Avg ~ 1, data = subset(hadcru, (Year <= 2009 & Year >= 1995)))

hadcru.mcmc <- MCMCregress(Avg ~ 1 + Year, data = subset(hadcru, (Year <= 2009 & Year >= 1995)), thin = 100, mcmc = 100000, b0 = c(-20, 0), B0 = c(.00001, .00001), marginal = "Laplace")
hadcru.mcmc.zero <- MCMCregress(Avg ~ 1, data = subset(hadcru, (Year <= 2009 & Year >= 1995)), thin = 100, mcmc = 100000, b0 = c(0), B0 = c(.00001), marginal = "Laplace")

-logLik(hadcru.lm)
## 'log Lik.' -14.55338 (df=3)
-logLik(hadcru.lm.zero)
## 'log Lik.' -12.80723 (df=2)

attr(hadcru.mcmc, "logmarglike")
##           [,1]
## [1,] -13.65188
attr(hadcru.mcmc.zero, "logmarglike")
##           [,1]
## [1,] -4.310564

alt text

|cite|improve this question
$\endgroup$
4
$\begingroup$

When you're computing Bayes factors, the priors matter. The influence of the priors can persist even if you have a large amount of data. When you're doing posterior inference, the effect of the prior goes away as you collect more data, but not so with Bayes factors.

Also, you'll get faster convergence if your null and alternative priors have disjoint support. Details here.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ your paper will take time to digest. it looks like the -logLik changes fairly significantly if I changed the number of digits in the prior precision. MCMCpack doesn't let me specify other priors (aside from uniform), so I'm going to try another package. btw, what is a "disjoint support?" $\endgroup$ – apeescape Aug 12 '10 at 1:55
  • $\begingroup$ On second though, my paper may not be so important here. But the point about the priors is important. Look at the differences in the priors at the empirical mean of the data. If that difference is large, that could explain your problem. $\endgroup$ – John D. Cook Aug 12 '10 at 2:10
  • $\begingroup$ As for the paper, disjoint support means there is no region to which both hypotheses assign positive probability. i.e. the alternative really is an alternative. $\endgroup$ – John D. Cook Aug 12 '10 at 2:12
  • 2
    $\begingroup$ it looks this problem has been widely known since 1939, called Jeffrey-Lindley paradox. Gelman's book (p185 & 250) also says uninformative priors w/ Bayes Factors are a no-no. Basically, they say that formulating point null hypotheses are stupid, and tests should be constructed from the posteriors. artsci.uc.edu/collegedepts/economics/research/docs/Wppdf/… $\endgroup$ – apeescape Aug 12 '10 at 4:50
  • 1
    $\begingroup$ Srikant: The problem is that you have two posteriors becoming more concentrated and you're looking at their relative fit. Either prior alone would become irrelevant in the limit, but the priors continue to impact the ratio of posterior model probabilities. $\endgroup$ – John D. Cook Aug 12 '10 at 13:56
4
$\begingroup$

I'm note sure I follow the R-code as I have only used R once or twice, but it looks to me as if you are comparing the marginal likelihood of a model with only an intercept and no slope (hadcru.mcmc.zero) and the marginal likelihood of a model with a slope and an intercept (hadcru.mcmc). However, while hadcru.mcmc.zero seems to be the correct model for H0, hadcru.mcmc does not seem to me to correctly represent H1 as there is nothing as far as I can see that constrains the slope to be positive. Is the something in the prior for the slope that makes it strictly positive (I don't know enough about MCMC in R to know)? If not, that may be where your problem lies as the marginal likelihood would then have a component representing the likelihood of the data for all of then egative values of the slop permitted under the prior (and 0) as well as the positive.

It is debatable whether the H0 for this question should be that the slope is exactly zero, nobody would believe that to be plausible a-priori. Perhaps a test using the Bayes factor for a model where the slope is strictly positive (H1) against a model where it is zero or negative (H0).

HTH (and I am not just confusing things)

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I think you have the interpretation of the R code correct (to be sure, hadcru.mcmc.zero has no slope parameter, it is just a mean model). I don't understand why the prior should be constrained just to the positive values though (I also wouldn't know how to code or write it down mathematically either; I want a diffuse prior but ignore negative values?). Shouldn't the H1 model include all values of the slope? Maybe I need to test whether H0: \theta <= 0 and H1: \theta > 0, in that case, I wouldn't know how to code it either :) $\endgroup$ – apeescape Aug 12 '10 at 17:07
  • 2
    $\begingroup$ The Bayes factor is the ratio of the marginal likelihood assuming H1 to the marginal likelihood assuming H0. In this case, the alternative hypothesis H1 is that the slope is positive, so your prior fo H1 must exclude negative values for the slope parameter as H1 says they are implausible a-priori. One way of achieving this might be to use an exponential prior, which is strictly positive, for the slope and a Gaussian prior for the intercept. You could then implement a "negative slope" null hypothesis by flipping the data upside-down and using the same model. HTH $\endgroup$ – Dikran Marsupial Aug 13 '10 at 7:23
1
$\begingroup$

I do not know the packages you are using or their internal working but perhaps the choice of priors matter? Perhaps, you should consider using different prior structures to see how sensitive the mcmc marginal likelihoods are to your choice of priors.

In particular, I suspect that the mcmc and the traditional likelihoods are likely to converge better as the priors become more diffuse. Note that in mcmc the marginal likelihoods are computed by integrating out the likelihood function with respect to the priors. Thus, I have a feeling that the 'diffuseness' of the priors may matter (could be wrong on this issue but worth checking out).

|cite|improve this answer
$\endgroup$
  • $\begingroup$ the priors are distributed normal N(b0, B0^-1). note above, that I use a precision (B0) of .00001, which is pretty diffuse. $\endgroup$ – apeescape Aug 12 '10 at 0:59
  • 1
    $\begingroup$ Actually, I think one other factor that matters is sample size. I am not sure if you can expect identical values with a sample size of just 15 odd points. With only 15 points there is not enough information in the likelihood to outweigh the prior. In other words, the likelihood is not 'concentrated enough' which suggests that when you integrate out the parameters you may get a different mcmc marginal likelihood. Does that make any sense? You can do a small simulation to check the above. $\endgroup$ – user28 Aug 12 '10 at 1:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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