Comparing two Bayesian models under disjoint prior supports using MCMC

Question

I have a Bayesian model involving three parameters $\theta_1$,$\theta_2$ and $\theta_3$. Experts think that $\theta_1 > \theta_2 > \theta_3$. So I would like to test the submodel $M_0$ corresonding to $\theta_1 > \theta_2 > \theta_3$ versus $ M_1:\{(\theta_1,\theta_2,\theta_3) \not\in M_0 \}$ regarding to the observations. I see two main options:

1) run the full model and check the percentage of samples for which $\theta_1 > \theta_2 > \theta_3$. Conclude to significance if its greater than $95\%$.

2) compute the bayes factor for the two models.

Are these two methods valid to you ? Are they equivalent in a certain sense ? One of my point is that it seems to me that there is no standard method to compute Bayes factor from MCMC sample without much trouble and point 1) looks as a good opportunity to me. Nevertheless any hints for bayes factor from MCMC samples (I used R and jags) are welcomed.

Answer 1

There are a few ways to evaluate this. First, you could run your model with priors that make no assumptions as to the relationship between $\theta_1$, $\theta_2$, and $\theta_3$, then examine the posterior distribution.

You could also use model selection methods to compare the above model to one with a prior that explicitly assumes $\theta_1 > \theta_2 > \theta_3$ (this can be done using the sort() function in JAGS). While you suggest Bayes Factors, there is quite a variety of methods that would work well. WAIC is becoming a widely used method for Bayesian model selection, and it can be computed much more easily in JAGS than Bayes Factors.

I would generally recommend going with the first option, but the second may be more valuable if you have a small amount of data.

Answer 2

So let call hypothesis $H_1$ is $\{(\theta_1,\theta_2,\theta_3) \mbox{ such that } \theta_1<\theta_2<\theta_3\}$, $H_2$ is it complementary. The Bayes factor writes: $$ K = \frac{P(x|H_1)}{P(x|H_2)} =\frac{P(H_1|x)}{P(H_2|x)} \frac{P(H_2)}{P(H_1)} $$

Then supposing that I have a model to infer $P(\theta_1,\theta_2,\theta_3 | x)$ under the complete parameter space $H_1 \cup H_2 $ then $P(H_1|x)$ can simply be computed as: $$ P(H_1|x)=\int_m^M P(\theta_1,\theta_2,\theta_3|x) 1_{\theta_1>\theta_2>\theta_3}(\theta_1,\theta_2,\theta_3) d\theta_1 d\theta_2 d\theta_3 $$ by simply counting the posterior MCMC samples satisfying the condition and $P(H_2|x)=1-P(H_1|x)$.

So it remains to extract the associated $\frac{P(H_2)}{P(H_1)}$ by simply considering the subpart of $H_1$ over the overall parameter space: $$ P(H_1)=\int p(\theta_1,\theta_2,\theta_3) \cdot 1_{\theta_1>\theta_2>\theta_3}(\theta_1,\theta_2,\theta_3) d\theta_1 d\theta_2 d\theta_3 $$ which can be computed analytically for the prior $p(\theta_1,\theta_2,\theta_3)$ associated to my posterior and again $ P(H_2)=1-p(H_1)$.

So the two options are very close each other.In particular if $P(H_1|x)>0.95$ then the BF would be greater than (0.95/0.05)*((5/6)/(1/6))=95