This seems like it should be a widely-answered, but I've spent a lot of time searching and reading textbooks and have not been able to figure it out.

Let's say I give my exciting new treatment pill to a random sample of N people and a placebo pill to a random sample of M people. A of the N people die and B of the M people die. I'd like to know the probability that more people die taking the placebo than the treatment.

Currently I do:

> data <- matrix(c(A,B,N-A,M-B), nrow=2)
> chisq.test(data, correct=FALSE)

and look at the p-value (which I realize isn't actually the answer to my question, but that's frequentists for you). And if I want to know the interval there's a 95% chance the actual value is in, I just do:

> pbar = A/N; se = sqrt(pbar * (1-pbar)/N); e = qnorm(.975)*se; pbar + c(-e, e)

But I want the Bayesian equivalents. Now I know I can compute a Beta distribution with a uniform prior and then use Monte Carlo sampling to figure out the answers, but that seems hackish and inefficient. What's the more direct way to compute the answer?

  • $\begingroup$ Note that in this specific example the p-value is all you need, provided you ignore your prior information about the two death rates. The Bayesian Hypothesis basically gives you a different way to choose the significance level. $\endgroup$ – probabilityislogic Jan 5 '13 at 13:22
  • $\begingroup$ Given your question is about life/death type of information, I would have thought taking a representative sample is way more important than bayes vs frequentist differences. Random sampling does not gaurantee you this in the specific case. This is part of that prior information I was talking about, where you control for known causes/differences between the people who survive and the people who die (age is usually an obvious example). $\endgroup$ – probabilityislogic Jan 5 '13 at 14:01

Well you could calculate the bivariate beta distribution of the posteriors of the 2 parameters (independent assumption seems reasonable), then find the distribution of the difference and see what the area greater than 0 is (or less than depending on the direction of interest).

Personally I think it would be more efficient (at least for me compared to looking up and remembering how to compute the distribution of the difference) to generate a bunch of observations from each of the beta posteriors, pair them up, and see what proportion have the one value in each pair being greater. Yes it is a bit hackish, but you make that sound like a bad thing when I don't think it is.

| cite | improve this answer | |

To take the second question first - what is a 95% credibility interval (no longer a confidence interval, now we are in the Bayesian world) for the probability that a person receiving the treatment will die (this seems to be the "actual value" you are calculating a confidence interval for in your line of code starting > pbar = A/N;).

As you seem to be aware, Laplace favoured using a uniform prior distribution (with p taking a value in [0,1]) for this sort of problem with proportions. A uniform distribution over [0,1] is a beta(1,1) distribution, which has the advantage of being a conjugate prior distribution for p; ie once you modify the distribution based on your observation, the resulting posterior distribution will have the same form as the prior - a beta distribution of some sort. Any textbook will show that this means the posterior distribution for p will be beta (A+1, N-A+1). From here you can just use this distribution as the basis for calculating your credibility interval - there is no need to use monte carlo techniques so long as you have access to tables or software that can calculate the relevant values from the beta distribution.

This will generally give you similar results to your frequentist approach relying on the central limit theorem for the distribution of pbar to approximate normal, as in your code snippet. See for example:

> options(digits=3)
> A <- 100
> N <- 10000
> # Frequentist approach, using approximation to Normal:
> pbar = A/N; se = sqrt(pbar * (1-pbar)/N); e = qnorm(.975)*se; pbar + c(-e, e)
[1] 0.00805 0.01195
> # Bayesian approach, with uniform prior:
> qbeta(c(0.025, 0.975), A+1, N-A+1)
[1] 0.00823 0.01215

Now, things are a little more complicated if you're really interested in the difference (not in what your one line snipped had) ie you want the posterior distribution and hence a credibility interval for the difference in probability of dying between the treatment and control populations, but the principles are the same. As you suggest and Greg Snow says in his answer, if you can't be bothered calculating the exact posterior distribution of the difference between the propostions why not use simulation? Certainly there's less paperwork involved; and it has the advantage of being easily extendible to more complex comparisons of the estimated posterior distributions of independent parameters:

> A <- 100
> N <- 10000
> B <- 200
> M <- 10000
> sim <- rbeta(10000, B+1, M-B+1) - rbeta(10000, A+1, N-A+1)
> quantile(sim, c(0.025, 0.975))
  2.5%  97.5% 
0.0067 0.0134 
> plot(density(sim), bty="l", main="Posterior distribution of difference\nbetween the proportions")
> abline(v=quantile(sim, c(0.025, 0.975)), lty=2)

enter image description here

Your first question was about hypothesis testing. Typically choosing between two hypotheses is done by comparing the posterior and prior odds ratios of the two hypothesis ie by calculating the Bayes factor:

$\frac{P(H_0|z)/P(H_1|z)}{P(H_0)/P(H_1)} $

where z is the observed data.

From here it gets a little complicated. Exactly how the two hypotheses are formulated is important - in a case such as yours it will generally be necessary to have $H_1$ indicating the treatment proportion is more than a specified amount different than the original - because under the prior distribution of the continuous random variable p the chance of any exact equivalence of the two is zero (a problem with "sharp" hypotheses). I'd consult a text to see how this plays out, but hopefully this was enough to get you started.

| cite | improve this answer | |
  • $\begingroup$ Good Q&A on hypothesis testing here: stats.stackexchange.com/questions/31679/… $\endgroup$ – Peter Ellis Jan 5 '13 at 1:51
  • $\begingroup$ Note that if you read the OP's question carefully/pedantically, their hypothesis is actually discrete and so is their parameter space as the question asks about whether more people in the population would die from taking the placebo. $\endgroup$ – probabilityislogic Jan 5 '13 at 13:54
  • $\begingroup$ Sorry if I wasn't clear, bit I do want to calculate the answer exactly, not simulate it. $\endgroup$ – aaronsw Jan 5 '13 at 15:27
  • $\begingroup$ @aaronsw you've got too many questions for 1 question I think. I see your 3 questions as: what is a credibility region for a point estimate of the proportion of people with the treatment who die (which is what you have done the frequentist confidence interval for, and I've given that); what is a credibility region for the difference between the treatment and control proportions (which reduces to the distribution of one Beta distributed random variable minus another Beta distributed random variable - if you want the exact distribution it will be complex); and how to do the hypothesis test. $\endgroup$ – Peter Ellis Jan 5 '13 at 17:15
  • $\begingroup$ @probabilityislogic good point; but it also depends on how 'population' is being used. If it is meant to include people not yet born you might say that it is still infinite, or at least very large and unknown number? $\endgroup$ – Peter Ellis Jan 5 '13 at 17:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.