comparing distributions - bayesian decision analysis

Question

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?

Answer 1

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.

Answer 2

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/