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$.
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.
R code to make the point explicit (hopefully).
$\endgroup$
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.