How to compare two bayesian linear models?

Question

Let say, we have two bayesian linear models learned from some data (not necessarily same but from the same data distribution), how can we compare or what is the notion of similarity between the two models?

For example, consider two bayesian linear models (second order) as given below: $$f(x) = \sum_{i=1}^n \alpha_i x_i + \sum_{i=1}^n\sum_{j=1}^n \alpha_{ij}x_i x_j \quad\quad (1)$$ $$f(x) = \sum_{i=1}^n \beta_i x_i + \sum_{i=1}^n\sum_{j=1}^n \beta_{ij}x_i x_j\quad \quad (2)$$

where each $\alpha_i$ and $\beta_i$ is distributed according to a posterior learned over data. My initial thought was to compare the KL-divergence of each $\alpha_i$ with corresponding $\beta_i$ but I am not sure whether it is a good idea.

To describe more context, once the models are learned, I would optimize the model over the input space $x$ to find the best point with the trained model as the objective function. This can be thought of as some form of thompson sampling. My overall goal is to compare how similar the two optimization problems will be over these two models.

Could you please explain some mechanism or insight to find how similar they are? I am also looking to extend the same idea to two gaussian processes trained on different sampled data sets but from the same data distribution.

Answer 1

This is more a comment than an answer, but I hope it may be of some help.

A difficulty, more frequent than one might think, in answering a probability & statistical question is that the question itself is actually ill-defined. And unfortunately our language doesn't help us in recognizing this problem: the question we ask is grammatically correct and seem at first sight to have a meaning – but it actually doesn't. The situation often becomes even worse when the question is dressed in technical terms.

"Which model is best?" – well, what do we mean by "best"? "How to compare these two models?" – well, what do we exactly mean with "compare"? (and what do we mean by "model"?)

We can introduce metrics of various kinds to make this kind of questions more precise. Then the problem is the rational justification of that particular metric. In the end there's often a lot of subjectivity swept under the carpet.

In your specific case, what I wonder first of all is why the whole dataset has not been used to obtain a unique distribution for the $\alpha$s (or $\beta$s), rather than dividing it into two copies of the same parametric model. The probability distribution for the model parameters would be sharper.

A second point is that if we are uncertain between two models $M_1$ and $M_2$, probability theory tells us that we should use a convex combination of the two, weighted by their probabilities given the data $D$: $$\mathrm{p}[x \mid (M_1 \lor M_2) \land D] = \mathrm{p}(x \mid M_1 \land D) \ \mathrm{p}(M_1 \mid D) + \mathrm{p}(x \mid M_2 \land D) \ \mathrm{p}(M_2 \mid D)$$ this is called "model averaging". See for example

• Draper: Assessment and Propagation of Model Uncertainty (also here)

for an early reference. I believe that in your case this is equivalent to using just one copy of your model, with all data determining a distribution for its $\alpha$s.

A third, final point is that if you want to choose a single $x$ value then you have a decision-theory problem, so you should choose and justify an appropriate utility function – which makes precise what you mean with "best" $x$ – and make your choice by maximization of expected utility. See for example chapters 13–14 of

• Jaynes: Probability Theory: The Logic of Science (also here and here)

Answer 2

I would argue that a standard Bayesian solution would make them directly incomparable, but I could be wrong. Consider two posteriors $$\pi(\theta|X_1)$$ and $$\pi(\theta|X_2).$$ Each posterior is optimal, given its observations and prior. It differs from a Frequentist method which optimal on average.

To give you a thought experiment as to why you cannot compare them directly, consider a school where you are going to split the students into two samples by tossing a fair coin. Your concern is estimating average height. Ignore the preposterousness of the setup for a moment. Assume you are in group one if the coin comes up heads.

Sample one’s estimators were $\bar{x}_1=36 in$ and $s^2_1=9in^2$. Sample two’s estimators were $\bar{x}_2=36.1$ and $s^2_2=8.9in^2$. What about flipping the coin, the intervention, caused students to be taller or the sample variability narrower when the coin came up tails? If nothing, then if you repeated this, say 1000 times, you could discuss the sensitivity of the parameter estimates to the sample space, but that is about it. You could construct the sampling distribution of the posterior means, medians, and modes or just one of them as they should be somewhat predictable. You can compare them indirectly as part of a distribution, but not directly because the difference is randomness.

A Frequentist process, of course, is built around the sample space rather than the parameter space and well understood analytic properties exist. While the Bayesian solution would be analytic under a conjugate prior, you will lack things such as the guaranteed minimum level against false positives.

Bayes factors will not help you here as there can be technical issues with Bayes factors you probably want to avoid. What I would do, if it were me, is that I would marginalize out the parameters leaving only a model parameter (Model 1, Model 2, etc.). That would give you an understanding of the impact of the sampling distribution on the trustworthiness of the posterior.

I have not looked, but as this is just a linear model, I would guess there is a paper on the Frequentist properties of Bayesian conjugate linear models. If it isn’t conjugate, there won’t be an analytic solution. Do note that your prior must be a proper prior, you will not be able to use a flat prior, because once your independent variables exceed 3, then a flat prior will prevent the posterior from integrating to unity and you may get spurious results. Still, it could have such a wide prior variance that it is effectively a flat line such as a standard deviation of one million centimeters over $\mu$ to measure the average size of a pea. Surely the size of the pea will be withing $\pm{3,000,000}cm$ of your prior estimate of 1 cm.

Because the difference is a coin toss, there can be no direct comparison. You can discuss the impact of the sample space on the sampling distribution of the parameters and some estimator of the likelihood function such as its supremum.