Bayesian inference and testable implications

Question

I have one question regarding testable implications of a model and Bayesian inference. My main doubt is how to exploit testable implications to reject a model.

Here is a simple example.

Suppose my model is that I have an iid sample from two gaussians with means $\mu_1$ and $\mu_2$ (and known variance) and I impose in the model the restriction that $\mu_1 = c \mu_2$ where $c$ is a known constant. Note this model imposes restrictions on the observed data---the mean of the two samples can't be more than $c$ times apart, barring, of course, sampling error.

Now suppose the data is discrepant with the posited model. Is there a principled way in Bayesian inference to reject this model given discrepant data?

PS: The generic comments saying that we could use posterior predictive checks, bayes factors and what not are ok, but not very useful in practice. Please also show how you would actually solve this toy problem.

Edit for the bounty: I will give the bounty to an answer that compares the three current suggestions (or more) both theoretically and more importantly with numerical examples of the toy problem. The three suggestions are: (i) posterior predictive checks; (ii) bayes factors; (iii) credible intervals (with or without ROPE).

For people potentially interested in answering: it would be helpful to actually perform a posterior predictive check and the hypothesis test or whatever you choose in your answer. The answer that does this and compares approaches will get the bounty.

Answer 1

There are only two "principled" ways you can get out of your posited model that operate within the framework of the Bayesian paradigm. Once is to initially set a broader class of models, and give some non-zero prior probability for the alternative models in that class (i.e., have a prior probability less than one for your posited model class). The other is to observe some evidence that has zero density under all distributions in the posited model class, which then allows you to update to any belief you want (see discussion here).

If you have assigned a prior probability of one to a class of models, and you never observe evidence that is inconsistent with those models, you can never "escape" that set of models within the Bayesian paradigm. Note that this is by design --- if you assign a prior probability of one to a set of models, you are saying that any alternative class of models has zero probability. In short, you are choosing to stick with your posited class of models no matter how strongly the evidence turns against them, so long as it is not inconsistent with those models. If you would like to have a principled "escape route" operating within the Bayesian paradigm, you will need to posit some broader class of alternative models and give it a non-zero prior probability. You could certainly give the alternative models a very low prior probability, so that they only become important a posteriori when the main model class starts to be (probabilistically) falsified by the data.

---

Implementation in your problem: In the problem you raise, it would be usual to handle this by framing the problem as a Bayesian hypothesis test, with hypotheses:

$$H_0: \mu_1 = c \mu_2 \quad \quad \quad H_A: \mu_1 \neq c \mu_2.$$

For example, under $H_0$ you could posit an overall model like this:

$$\begin{aligned} X_{11}, X_{12}, ... , X_{1n} | \mu_2,\sigma_1^2,\sigma_2^2 &\sim \text{N}(c \mu_2,\sigma_1^2), \\[6pt] X_{21}, X_{22}, ... , X_{2n} | \mu_2,\sigma_1^2,\sigma_2^2 &\sim \text{N}(\mu_2,\sigma_2^2), \\[6pt] \mu_2 &\sim \text{N}(0, \eta^2), \\[6pt] \sigma_1^2 &\sim \text{Ga}(\alpha, \beta), \\[6pt] \sigma_2^2 &\sim \text{Ga}(\alpha, \beta), \\[6pt] \end{aligned}$$

and under $H_A$ you could posit an overall model like this:

$$\begin{aligned} X_{11}, X_{12}, ... , X_{1n} | \mu_1,\mu_2,\sigma_1^2,\sigma_2^2 &\sim \text{N}(\mu_1,\sigma_1^2), \\[6pt] X_{21}, X_{22}, ... , X_{2n} | \mu_1,\mu_2,\sigma_1^2,\sigma_2^2 &\sim \text{N}(\mu_2,\sigma_2^2), \\[6pt] \mu_1 &\sim \text{N}(0, \eta^2), \\[6pt] \mu_2 &\sim \text{N}(0, \eta^2), \\[6pt] \sigma_1^2 &\sim \text{Ga}(\alpha, \beta), \\[6pt] \sigma_2^2 &\sim \text{Ga}(\alpha, \beta). \\[6pt] \end{aligned}$$

You can obtain the Bayes' factor for the above hypothesis test and use this to see how you update prior probabilities for the hypotheses to posterior probabilities. If the data makes $H_0$ highly implausible, this will manifest in a lower posterior probability for $H_0$. Given some prior probability $\lambda = \mathbb{P}(H_0)$ for your posited subclass of models, you will be able to update this to a posterior probability.

Answer 2

Prior predictive and posterior predictive checks may be helpful in here. In both cases you sample the predictions from the model (the "fake data"), in first case from the prior, in the second case from the posterior distribution, and then compare the distributions of the fake data, with the distribution of the observed data. Prior predictive checks are aimed to diagnosing the prior-data conflict, i.e. the model a priori does not make reasonable predictions that cover the possible range of the values observed in the data, it is ill-defined a priori. In posterior predictive checks you sample from the predictions after estimating the parameters (i.e. from posterior), so you check if the predictions that the model does fit the observed data. In both cases, there are many ways of doing this, depending on particular problem, ranging form eyeballing the histograms, density plots, scatter plots, summary statistics etc, up to defining more formal tests (data falls within the per-specified interval, hypothesis tests to compare the distributions, etc). This is a routine practice in Bayesian modeling.

If I understand you correctly, the model that you use as example assumes that your data $X$ comes from a mixture of two Gaussians, with unknown means $\mu_1, \mu_2$ and known variances $\sigma^2_1, \sigma^2_2$, and known constraint $c$, such that $\mu_2 = c\mu_1$. Simple way to test this model is to treat $c$ as free parameter, to be estimated. You know what $c$ should be, so you can come up with a strong, informative prior for it. In such case, it would surprise you if estimated $c$ differed from the true value. If I understand you correctly, that's the property of the model that you want to test. To test the validity of this assumption, you could take samples from the posterior distribution $\hat c_i$, and compare them to the true value of $c$, e.g. you would accept the model if at least in in $100\alpha\%$ cases, the predicted values for $c$ would be within the $\pm \varepsilon$ range from the truth

$$ \alpha \le 1/n \sum_{i=1}^n \mathbf{1}(|c - \hat c_i| < \varepsilon) $$

This is not exactly a posterior predictive check, since we may argue if $c$ is data, or not, but it follows the spirit of the kind of checks you would make to test model validity.

Accidentally, Michael Betancourt has just published a lengthy Towards A Principled Bayesian Workflow tutorial, where among other things, he discusses importance of prior and posterior checks discussed above.

Answer 3

I'm not a Bayesian expert and I'm happy to stand corrected, but to me the most straightforward & principled way to test this would be to define region of practical equivalence (ROPE) around c and then estimate how much posterior density falls inside this region.

For example, let's say that, based on theory and domain knowledge, you know that for all practical purposes, if c deviates from exactly 1 by less than 0.01 then it might as well be 1 (outside of simulation, c is never going to be exactly 1 anyway and so you will always reject point null hypothesis with enough data). Anyway, using the deviation of 0.01 you define a ROPE of 0.99 - 1.01. After that, you run your model, and estimate how much density falls inside the ROPE region. If the proportion of density $k$ that falls inside the rope is smaller than whatever you decide your alpha is, then you should feel comfortable rejecting your model, with $k$ confidence. See this vignette: https://easystats.github.io/bayestestR/articles/region_of_practical_equivalence.html

PS: You'll probably want a large tail effective sample (ESS) size for this kind of testing. This is because Monte Carlo samplers tend to explore the typical set & give increasingly less precise estimates towards the tails of the distribution, which is where your ROPE might be. So you'll want to run your sampler with a lot of iterations.

Answer 4

EDIT: innisfree is right. Bayes factors seem like a better approach than what I have provided here. I'm leaving it up for posterity, but it isn't the right approach.

Because this problem really relies on a single assertion (namely, that $c$ has some value), we can simply estimate the following model

$$ y \sim \mathcal{N}(b_0 + b_1x, \sigma)$$

and determine the posterior probability that either $b_0/(b_0+b_1)<c$ or $b_0/(b_0+b_1)>c$. Here is an example. Say we had a hypothesis that $c=1$ and we know that the variance is 4 and that the intercept (or the mean of one populatio

n) is 2. We can fit the following model in Stan

stan_model = '
data{
  int n;
  vector[n] x;
  vector[n]y;
}
parameters{
  real b;
}
model{
  b~normal(0,1);
  y~normal(2+b*x, 2);
}

'

This will allow is to freely estimate the parameter $b1$ assuming we know $b_0$ and $\sigma$. After fitting the model with a standard normal prior on $b_1$, here is a histogram of the posterior

The model provides a 95% posterior credible interval for $c$ covering (0.465, 0.686). We can be fairly certain that the value of $c$ is not 1.