Obtaining a Bayes Factor for the difference between two proportions (R code provided) Below (in R code), I'm showing the Bayesian estimation of the Difference between two proportions resulted from two binomially distributed groups (groups) of scores (y).
My question is how can I now compute a Bayes Factor to understand which hypothesis (i.e., $H_1$ or $H_0$) my data supports?
set.seed(0)
n1 = 100 ; n2 = 100 ; p1 = .65 ; p2 = .5 ; trials = 15

     y = as.vector(unlist(mapply(FUN = rbinom, n = c(n1, n2), size = c(trials, trials), prob = c(p1, p2))))
groups = factor( rep(1:2, times = c(n1, n2)) )

alpha1 <- 500 ; beta1 <- 500 ; alpha2 <- 300 ; beta2 <- 500 # Beta priors parametrs for each group of scores

post.p1 = rbeta(1e6, alpha1 + sum(y[groups == 1]), beta1 + n1*trials - sum(y[groups == 1]))
post.p2 = rbeta(1e6, alpha2 + sum(y[groups == 2]), beta2 + n2*trials - sum(y[groups == 2]))

delta.p = post.p1 - post.p2

hist(delta.p, freq = FALSE, col = 4)  ;  (CI = quantile(delta.p, c(0.025, 0.975)))

 A: I am assuming from your question, code, naming conventions, and comment that you are working with a beta prior on a binomial likelihood. Let's start by defining the Bayes Factor (BF) which is 
$BF_{12} = \frac{Pr(D|H_1)}{Pr(D|H_2)}$. Here I am using the notation from the Kass and Raftery http://www.stat.cmu.edu/~kass/papers/bayesfactors.pdf paper a worthwhile read if you'd like to learn more about Bayes Factors. Now the question is how to calculate $\Pr(D|H_i)$? To do this you need your likelihood and prior. Assuming binomial likelihood and beta prior you calculate the integral, good notation here helps. Here is the integral: 
$\Pr(D|H_j)= \int_0^1 \Pr(D|p,H_j)\pi(p|H_j)$. In this case the 2 separate hypthesis are a bit different in approach. The first one assumes the two $p_i$s are the same therefore the priors must have the same values of the parameters $\alpha$ and $\beta$. So the Typical Bayesian calculus of prior times likelihood ad integrating out the $p$ parameter yields
$m1 ={n_1+n_2\choose x_1+x_2}\frac{ \Gamma(\alpha+\beta) }{ \Gamma(\alpha) \Gamma(\beta) }\frac{ \Gamma(x_1+x_2+\alpha)\Gamma(n_1+n_2-x_1-x_2+\beta)}{\Gamma(n_1+n_2+\alpha+\beta )}$. 
Note that this result is obtained by incorporating the null hypothesis of equality ($p_1=p_2$) into the likelihood and prior. 
For the second hypothesis you have no constraint on $p_i, i=1,2$. Although not stated in the OP I will assume independence of $p_i$s. Therefore, a similar integral of the product of 2 likelihoods (but this time with different $p_i$s and 2 priors gives
$m2 ={n_1\choose x_1}{n_2\choose x_2}\frac{ \Gamma(\alpha_1+\beta_1) }{ \Gamma(\alpha_1) \Gamma(\beta_1) }\frac{ \Gamma(\alpha_2+\beta_2) }{ \Gamma(\alpha_2) \Gamma(\beta_2) }\frac{ \Gamma(x_1+\alpha_1)\Gamma(n_1-x_1+\beta_1)}{\Gamma(n_1+\alpha_1+\beta_1 )} \frac{ \Gamma(x_2+\alpha_2)\Gamma(n_2-x_2+\beta_2)}{\Gamma(n_2+\alpha_2+\beta_2 )}$
Now these two values ($m_i,i=1,2$) can be easily calculated for the specific values of $n_i$, $\alpha_i$, $\beta_i$, for $i=1,2$ in the OP. This can be done in R language for example via the corresponding choose() and Gamma() functions to calculate $m_i$ and then divide $m_1/m_2$. 
