6
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ One method that comes to mind is the Bayes factor, though it is generally not recommended. $\endgroup$ May 16, 2020 at 23:57
  • $\begingroup$ One does not simply reject a model. The question assumes it makes sense for Bayesians to reproduce frequentist artifacts. At best you will end up with a posterior distribution over models. The most Bayesian thing to do is hang on to that full posterior distribution, not apply some arbitrary procedure to it so that one hypothesis is 'rejected' and the other (wink, wink) 'not rejected'. $\endgroup$ May 18, 2020 at 9:36
  • 1
    $\begingroup$ Can you clarify the sentence "I impose in the model the restriction that $\mu_1=c\mu_2$ [...] the mean of the to samples can't be more than $c$ times apart". This seems to imply $\mu_1<c\mu_2$. Also, is $c$ known outright or do you place a prior on it? $\endgroup$ May 18, 2020 at 14:25
  • $\begingroup$ @DemetriPananos It means that int my model the restriction $\mu_1 = c\mu_2$ holds. In the real world anything can hold. $\endgroup$
    – user272422
    May 18, 2020 at 15:34
  • $\begingroup$ @conjectures Go ahead then and write down your answer of how you should solve this toy problem. $\endgroup$
    – user272422
    May 18, 2020 at 15:36

4 Answers 4

7
+50
$\begingroup$

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.

$\endgroup$
11
  • $\begingroup$ Can you numerically apply your proposal in my simple two gaussians example? $\endgroup$
    – user272422
    May 18, 2020 at 0:56
  • $\begingroup$ I'm asking because both your and Tim's answer are very generic. It should be easy to demonstrate the two approaches in this example: the posterior predictive check versus having a hyper prior or something (is that what you are proposing, a hyper prior, one that selects the constrained model $M_1$ with X% chance versus the unconstrained model $M_2$ with 1-X% chance?). $\endgroup$
    – user272422
    May 18, 2020 at 0:59
  • 2
    $\begingroup$ In a basic sense, every assumption restricts things, so you just need to list your assumptions and then consider if you want to test these, and if so how. One of the difficulties of the Bayesian framework (as opposed to classical hypothesis testing) is that you have to specify your alternative as a model, rather than just leaving it as non-compliance with the null. Nevertheless, if you can identify the assumptions you want to test, and the alternative, you can usually set up an alternative model. $\endgroup$
    – Ben
    May 18, 2020 at 3:12
  • 3
    $\begingroup$ @BayesianNewbie: You are certainly free to award the bounty to whomever you wish, on whatever basis you wish. I will just note that it is considered somewhat "bad form" to move the goalposts with edits requiring substantially more work after a question has already been asked and answered. (Since you are new to the site, you probably didn't know that, so please don't take this as a criticism.) In any case, numerical implementation of the model is a much larger project, that is outside the scope of my answer. Cheers, Ben. $\endgroup$
    – Ben
    May 19, 2020 at 2:43
  • 1
    $\begingroup$ Ben I gave the bounty to you since you were going to get it anyway, but this answer is too generic... I'm going to ask several more specific questions next to complement this. $\endgroup$
    – user272422
    May 24, 2020 at 1:00
4
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ May 18, 2020 at 16:52
  • $\begingroup$ Please post any new comments in the dedicated chat room linked above. All comments posted here will be deleted from now on. $\endgroup$ May 18, 2020 at 17:04
2
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Will this be equivalent to posterior odds for the hypotheses $H_0: c \in ROPE$ versus $H_1: c \not\in ROPE$ with prior for $p(H_0)=\int_{ROPE} p(c) dc$ etc $\endgroup$
    – innisfree
    May 18, 2020 at 8:38
  • 4
    $\begingroup$ @innisfree you’re doing a lot of commenting, which is great, but it seems you have a preferred way of answering this question. I think you should post it to compare among the existing answers. $\endgroup$ May 18, 2020 at 14:13
  • 1
    $\begingroup$ Adam, it doesn't ROPE doesn't matter that much here. You are just changing the constraint to $|\mu_1 - \mu_2| < c$. But the logic is the same. $\endgroup$
    – user272422
    May 18, 2020 at 15:45
  • $\begingroup$ I more of less agree with the answer of @ben $\endgroup$
    – innisfree
    May 18, 2020 at 20:57
  • $\begingroup$ @BayesianNewbie yeah I think you're right, if I understand you correctly, the only real difference between mine & Demetri's answer is that you quantify what proportion of the total density falls inside the ROPE around c, as opposed to just concluding that c falls outside of the 95% HDI. I don't think there's going to be much difference for symmetric, two-sided intervals, but I think that I can imagine a situation where c would fall inside the HDI but at the same time have < alpha proportion of density in the ROPE (marginal posterior with long tails). $\endgroup$
    – Adam B.
    May 18, 2020 at 23:04
1
$\begingroup$

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

enter image description here

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.

$\endgroup$
5
  • $\begingroup$ You can’t use credible regions to test a point null. The fact that $1$ lies outside a credible region tells nothing about the relative plausibility of $c=1$ to $c\neq 1$ or an error rate associated with any decision procedure associated with $c$. $\endgroup$
    – innisfree
    May 17, 2020 at 6:33
  • $\begingroup$ On top those those two problems, the procedure depends on an ordering rule for the credible region, and worse still, isn’t necessarily invariant under reparameterisations of the problem. $\endgroup$
    – innisfree
    May 17, 2020 at 6:37
  • 1
    $\begingroup$ Thanks, I don't think innisfree is correct, your answer is sensible, but I believe your modeling choice is not correct for this example. Here you should model both means. $\endgroup$
    – user272422
    May 17, 2020 at 19:28
  • $\begingroup$ @innisfree I don't understand your objection to this answer. Can you elaborate on how you would go about solving this question in my example on a separate answer? $\endgroup$
    – user272422
    May 17, 2020 at 19:37
  • 1
    $\begingroup$ Why couldn't you just estimate how much posterior density falls in region of practical equivalence (ROPE) around 1? Point null c = 1 is always going to be wrong. $\endgroup$
    – Adam B.
    May 18, 2020 at 4:35

Your Answer

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