# Posterior distribution of proportional difference of two binomial variables

Can somebody point me in the right direction for a treatment of the following problem? I imagine this should be a fairly common problem in medical statistics...

Given two binomial random variables $X_1\sim Bin(n_1,\pi_1)$ and $X_2\sim Bin(n_2,\pi_2)$ I am looking for the posterior distribution of $\pi_1-\pi_2$.

In case this is relevant for simplifying approximations, in my problem both $\pi$ are approximately 2%, their difference about 0.1%. Both $n$ are of the order of $1000$.

• You need to tell us what prior distributions you are placing on $\pi_1$ and $\pi_2$
– Dan
Commented Aug 27, 2014 at 21:56
• If you place beta priors on $\pi_1$ and $\pi_2$ then the posterior distribution of $\pi_1$ and $\pi_2$ will be the difference of two beta distributions which have a closed form expression.
– Dan
Commented Aug 27, 2014 at 21:57
• @Dan I really have no idea. I am a physicist and did an experiment. In the first situation I am getting $n_1\times \pi_1$ occurences for $n_1$ trials, in the second situation I am getting $n_2\times \pi_2$ occurences for $n_2$ trials. I am looking for the probability of the two situations differing by a given amount. - Does this make any sense at all?
– Max
Commented Aug 27, 2014 at 22:03
• Yes it makes sense, but you still need to choose a distribution to place on the parameters you are estimating. This is always a tough and subjective decision. The default, especially since you have a large amount of data, is to place and "non-informative prior" which in this case is a Beta(1,1) which conveniently is the Uniform(0,1) prior.
– Dan
Commented Aug 27, 2014 at 22:06
• @Dan Re-reading your comments, somehow I am still not understanding you quite right: Assuming $\pi_1\sim$Uniform(0,1) and $\pi_2\sim$Uniform(0,1) the posterior distribution $P(\pi_1-\pi_2)\sim$Uniform(0,1)-Uniform(0,1). This uniform difference seems to rise linearly from (-1,0) to (0,1) and fall from there linearly to (1,0). But this does not incorporate the likelyhood of my data at all...
– Max
Commented Aug 27, 2014 at 22:29

Since you say that you have a large amount of data, it would be "fairly logical" (or at least easy to defend) the use of a diffuse prior, which in this case is as simple as placing a Beta(1,1) which conveniently is the conjugate prior for the binomial likelihood.

In general, if we have that $X\sim \text{Bin}(n,\pi)$ then placing a $\text{Beta}(\alpha,\beta)$ prior on $\pi$ yields the following:

\begin{align*} p(\pi|X)&=\frac{p(X|\pi)p(\pi)}{\int p(X|\pi)p(\pi)d\pi}\\ \,\\ &\propto p(X|\pi)p(\pi)\\ \,\\ &\propto \prod_{i=1}^N\text{Bin}(n,\pi)\text{Beta}(\alpha,\beta)\\ \,\\ &\propto \prod_{i=1}^N\pi^{X_i}(1-\pi)^{n_i-X_i}\pi^\alpha(1-\pi)^\beta\\ \,\\ &\propto \pi^{\alpha+\sum_{i=1}^NX_i}(1-\pi)^{\beta+\sum_{i=1}^Nn_i-\sum_{i=1}^NX_i} \end{align*} which is proportional to another beta distribution and so the posterior distribution for $\pi$ is $\pi|X\sim\text{Beta}(\alpha+\sum_{i=1}^NX_i, \beta+\sum_{i=1}^Nn_i-\sum_{i=1}^NX_i)$. You can verify these results here: http://en.wikipedia.org/wiki/Conjugate_prior.

Now to make the above have an uninformative prior then you lets $\pi\sim\text{Beta}(\alpha=1,\beta=1)$ which is really saying $\pi\sim\text{Uniform}(0,1)$.

Now, returning to your question, we have that independently the posterior distributions for $\pi_1$ and $\pi_2$ are

$$\pi_1|X_1\sim\text{Beta}\left(1+X_1,\, 1+n_1-X_{1}\right)$$

and $$\pi_2|X_2\sim\text{Beta}\left(1+X_2,\, 1+n_2-X_{2}\right)$$

and so now, $\pi_1-\pi_2|X_1,X_2$ should just be another distribution with updated parameters (I don't know the distributional parameters of the top of my head but I am sure its easy enough to google for difference (or sum) of beta random variables).

However, if you don't care about a closed form solution you can sample from the posterior distribution of $\pi_1$ and $\pi_2$ independently and then literally subtract the posterior sample from each other to obtain posterior samples of $\pi_1-\pi_2$. Algorithmically, what I mean by this is the following:

Step 1) Sample $\pi_1^*$ from $\pi_1|X_1$

Step 2) Sample $\pi_2^*$ from $\pi_2|X_2$

Step 3) Obtain posterior samples of $\pi_1-\pi_2$ by calculating $\pi_1^*- \pi_2^*$

Not sure if this is helpful or not but here is some R code validating what I am proposing to you. The code is the following:

#Both pi's are approximately 2% and their difference is 0.1%
pi1 = .02
pi2 = .019

#Sample sizes of 1000
n1 = 1000
n2 = 1000

#X1 and X2 equal to their pi*n
x1 = pi1*n1
x2 = pi2*n2

#Sampling 10,000 random variables from the posterior distributions of pi1 and pi2
post1 = rbeta(10000,1+x1,1+n1-x1)
post2 = rbeta(10000,1+x2,1+n2-x2)

#Calulating posterior samples from pi1 - pi2
posterior = post1 - post2


Now we know that the true difference between $\pi_1-\pi_2=0.02-0.019=.001$ (in my computer code example) so now taking the mean (average) of the posterior samples I obtain for $\pi_1-\pi_2|X_1,X_2$ we obtain

> mean(posterior)
[1] 0.001002589


which is extremely close to the 0.1% that we expected to see.

• Dan, your first section on the conditional posteriors for $\pi_1$ and $\pi_2$ looks good (and on that basis, +1), but your subsequent statement "so now $\pi_1-\pi_2|X_1,X_2$ should just be another Beta distribution with updated parameters" - I don't see how that is justified. It's the distribution of the difference of two beta r.v.s (as you said in your previous comment) but that isn't itself beta, not even rescaled beta. Consider, for example, the distribution of the difference of two uniforms, which is triangular. Commented Aug 27, 2014 at 22:45
• Thanks for this awesome answer. This really clears things up. Given the issue with the sum of two beta distributions, I am happy to go with the sampling solution you suggested. Am I understanding right that I will need to carry out steps 1-3 several times to obtain a distribution for $\pi_1-\pi_2|X_1,X_2$ rather than a point estimate? I will thus implement nested sampling, e.g. M times each in steps 1 & 2 and N times the whole procedure giving a total of 2*M*N samples?
– Max
Commented Aug 27, 2014 at 22:47
• @Glen_b you are completely correct and I updated my answer to reflect that. Thanks!
– Dan
Commented Aug 27, 2014 at 22:52
• @Max I included some R code to make the point explicit (hopefully).
– Dan
Commented Aug 27, 2014 at 22:53
• @Max And no, you sample from steps 1 and 2 and then subtract 1 from 2 to obtain one sample from 3. You do this M times to obtain M samples.
– Dan
Commented Aug 27, 2014 at 22:55

Dan suggested using simulation, which is fine, and another alternative would be to do the convolution itself (one can find results for sums of beta random variables in various places, from which the present result can be obtained); alternatively it can be done numerically fairly easily. Dan asked me to post about the normal approximation I suggested in comments, so I will expand on that possibility here.

The process follows his derivation out to the posteriors for the individual $\pi$'s.

At that point we notice that (for say a uniform prior or a Jeffreys prior or some similar conjugate prior) we have a pair of $\text{beta}(a_i,b_i)$ with $a$'s near 20 and $b$'s near 1000.

They look somewhat like a normal, though very mildly skew, something like this:

The distribution of the difference of two such, where $a_1$ is close to $a_2$ will be much more symmetric and the tails will be much nicer. So a normal approximation is probably going to work well.

Simulation confirms this. Here's a histogram of a large sample from the difference between a beta(20,1000) and a beta(19,1000):

and a Q-Q plot:

(those two displays are from different samples, not that you can tell)

Given the posteriors for the $\pi$'s, we have the posterior mean and variance for them immediately (See here, which gives them anyway).

From there we just use independence, plus elementary properties of mean and variance to get the mean and variance of the difference (the mean is the difference of means, the variance is the sum of the variances).

Once we have the mean and variance of the posterior of the difference in proportion, we simply take the posterior distribution to be a normal with that same mean and variance and then use that approximation to do whatever calculations are required.

If the samples had been smaller, or the proportions were smaller/more different, I'd be inclined to suggest considering a scaled-shifted-beta approximation (matching the first 4 moments, say), though if we go to that effort, convolution might be as easy.