comparing distributions - bayesian decision analysis

I am attempting to use Bayesian analysis to compare distributions to help with decision analysis - when to treat a patient based on a blood measurement X.

Here you can see 1000 samples from two posterior distributions (student t) that I am trying to compare. I have computed the difference between means, and p(mu_delta > 0 | data) = 0.335 --> not entirely convincing.

That said, the blue distribution has a much longer tail on the right (higher level of X) than the red distribution. The blue distribution is patients who shouldn't be treated. The aim is to reduce the number of patients who receive treatment who shouldn't. Is there a way I can select a threshold to say any marker > VALUE(X) should not receive treatment since there is a higher probability they will not respond? Or is it only fair to compare means?

I understand that this is a sample distribution which may not behave in the same way with real data but that is of course something I can test.

If it makes a difference - this work was done using SciPy (python) and I have other measurements I wold also like to include (multivariate) using Stan.

EDIT:I have since realised that I am also sampling values below 0 which cannot actually happen. Is there a way to correct for this or should I use rejection/ MCMC sampling to reject all points below?

• The answers below are quite good, so I just wanted to comment that using rejection sampling after the fact is not a good idea. It implies that whatever priors you were using do not integrate to 1 (or another constant) over the admissible parameter space (in this case, the values of the parameter that imply the histogram is positive). And if your priors do not integrate to 1, then the posterior is not coherent even before you start throwing away draws from it. With Stan, you should be able to constrain things appropriately from the start and specify a prior that integrates to 1 over that space. Apr 11, 2018 at 16:31
• @BenGoodrich: sorry, I do not understand the comment! Even when proposing from a prior (which then necessarily integrates to one), rejecting values outside the range of possible values is in general valid. Jun 1, 2018 at 8:44
• @JB1: I have trouble with the question in that I do not understand the difference between the two posteriors: is it against a different prior? a different dataset? In both cases the natural tool would be a Bayes factor. Jun 1, 2018 at 8:47

The right way to go about this is to compute the expected benefit of giving treatment and the expected benefit of not giving treatment, and then take whichever action has the greater expected benefit. A good introduction is "Making Hard Decisions" by Robert Clemen.

In more detail here it is. Let U(blue, T) = benefit for giving treatment when patient is blue, U(blue, not T) = benefit for not giving treatment when patient is blue, U(red, T) = benefit for giving treatment when patient is red, U(red, not T) = benefit for not giving treatment when patient is red. Any of the "benefit" values U(blue, T), etc, can be negative to indicate it's actually a loss. If you can assign reasonable values to these, you can go ahead.

The expected benefit for T is EU(T) = U(blue, T) p(blue | X) + U(red, T) p(red | X), and the expected benefit for not T is EU(not T) = U(blue, not T) p(blue | X) + U(red, not T) p(red | X). So the recommendation is to give treatment if EU(T) > EU(not T) and not otherwise.

You'll want to be careful with plus and minus signs throughout; figure out a sign convention and stick to it. For example, so far I'm assuming losses are negative and benefits are positive, and you want maximize benefit. You can work it the other way too, you just have to be consistent.

You can get some insight by simplifying the inequality a little. You can derive something like p(red|X)/p(blue|X) is less or greater than some ratio of benefits. Be careful about the signs of the benefit differences when you do that, since it changes the sense of the comparison.

The histograms you showed are p(X|red) and p(X|blue) -- from these you can get p(red|X) and p(blue|X) via Bayes' rule. You'll need to estimate the background rates p(red) and p(blue) in order to do that.

The most obvious take-away from this whole analysis is that if you have p(red|X)/p(blue|X) less or greater than some ratio of benefits, you can rework that into something involving p(X|red)/p(X|blue), and that means you can choose a place on the overlapping histograms plot and say, apply treatment if the ratio of densities is more or less than a threshold. That threshold ratio will be a function of the benefits and the background rates.

• Thanks for the response. The U(T) concept is a nice one and not one I am familiar with although is very intuitive (my brain was going in circles i think). Im not sure if it was clear from the description, but assuming U(T) is a percentage, U(red, T) = 1 and U(blue, T) = 0. This means that EU(T) is just the p(red | X), and of course the counter. the example here is with a single parameter. I shall extend to multiparametric but of course the same principle holds. Thanks again :)
– JB1
Apr 9, 2018 at 10:58
• The "benefit" function is abstractly called "utility", which is the motivation for choosing U instead of B for the name of the function. It's hardly any more complicated to work with general values so I encourage you to go ahead and work through the problem without assuming 0/1 values. I am guessing it you will get some useful insight about the problem if you look at ranges of possible values after working through it, instead of substituting 0/1 up front. Apr 9, 2018 at 14:44
• Thanks @robert-dodier. What I mean by 0/1 values is that these values are known labels, and I have split the population into the 0/1, sampled and visualised. This is maybe what you already understood but im just checking! By this I mean, the red distribution has been sampled from the posterior of patients who DO respond to treatment, and blue do not. I will work through the problem to get a ratio and see how things look. Thanks for you r help!
– JB1
Apr 9, 2018 at 19:18
1. if you're getting values below zero for the measurements (I assume you're referring to them and not to your chart, as the difference in means seems like something which can perfectly well be negative), you probably want to change your parameterization from a normal/student-t to something like a lognormal, which is always positive and otherwise will look more or less like a normal with fatter tails and seems sensible for a medical application where you might have sick patients with occasional extreme values. I don't know if this is very important to your application, but if nothing else, it'll make the results a little more sensible and your MCMC will probably run a little faster. (If the lognormal or exponential or other positive distributions don't work for you, I think Stan does let you add some constraints to ensure >=0 values. No idea about SciPy.)
2. if you want a decision analysis, you need to do more work and define a loss function: what is the cost/gain of treating a patient when they should be treated? what is the cost of not treating a patient when they should be treated? what is the cost of not treating a patient when they should not be treated? what is the cost of treating a patient when they should not be treated? When you define those, you can simulate various treatment thresholds in various scenarios, weighted by their posterior probability, and see what thresholds reduce cost the most.

Indeed, sometimes a decision problem is a foregone conclusion: it's entirely possible that you'll learn that no patients should be treated based on this biomarker, or, vice-versa, that regardless of biomarker status all patients should be treated. (If it's really expensive, or really cheap, say.) Your description is inadequate for any kind of answer here; if you really want a decision analysis, you may need to take a look at some textbooks and papers to get an idea of how to actually go about it for your example. You can find a bunch here: https://www.reddit.com/r/DecisionTheory/

• "It's entirely possible that you'll learn that no patients should be treated based on this biomarker." This is an important point to bear in mind. Apr 8, 2018 at 23:09
• Thanks for your response! with regards to prior - youre right a log normal makes much more sense, I will switch to this. I can see that Stan also accepts contraints. I will move to stan now for MCMC.
– JB1
Apr 9, 2018 at 8:27
• With regards to decision analysis and loss functions I think maybe I should put this on hold (slightly). My main aim is to first check if the distributions between the two populations are different using a Bayes approach, and if so, use this information practically. I understand that based on a biomarker we may learn no patients should be treated (ethically one could argue all patients should be treated?) however it would seem its possible to remove patients at the higher end of blue (>8.0)? I will trawl the subreddit for a while :)
– JB1
Apr 9, 2018 at 9:31
• @JB1 "My main aim is to first check if the distributions between the two populations are different using a Bayes approach, and if so, use this information practically." -- this isn't really consistent. The decision theoretic stuff about expected loss is exactly the formalization of the practical aspect. It doesn't require quantifying how different the distributions are in order to get to a decision. ... "it would seem its possible to remove patients at the higher end of blue (>8.0)?" -- agreed, it is possible, and the stuff I sketched about the ratio p(X|red)/p(X|blue) is the way to get there. Apr 9, 2018 at 14:51