Correct usage/understanding of Bayes Factor when comparing two proportions I'm just starting to learn R and explore Bayesian statistics, but I keep getting tripped on using Bayes Factor and (honestly), I'd love a little confirmation if my process is correct in interpreting difference between two sample proportions. I'd definitely appreciate any feedback.
Below is a scenario I'm looking at, and I'd love to know if the process below in R is right or if I'm applying anything incorrectly.
Scenario:
Web users are split into two groups to test how well a website design performs.  Since the designs are new, I'll use an uninformed prior (1,1). After running the split test for a few days, I have the following results. Results can only be successes/failures, so looking to use a binomial-beta.
Data:
Sample 1:


*

*n = 1000

*successes = 20

*failures = 980


Sample 2:


*

*n = 2000

*successes = 30

*failures = 1970


I want to use the following code to compare the resulting distributions and then determine the Bayes Factor.
theta1=rbeta(10000,20+1,980+1) #taking 10,000 random draws from distribution with sample 1
theta2=rbeta(10000,30+1,1970+1) #taking 10,000 random draws from distribution with sample 2
theta_dif = mean(theta1>theta2) #Find the posterior probability that someone from sample 1 will convert more than someone in sample 2
bf = theta_dif/(1-theta_dif) #This is the step I'm most unsure of is taking the probability that sample 1 is better over the complement (which is the probability that No. 2 is better) the right way to get a Bayes Factor? 

Am I correct in these steps? In this case, the result comes out to be 5.01 which would be "Substantial". The 5.01 represents the roughly 85% / 15% from theta_dif (posterior probability that sample1 is better) divided by its compliment (theta2). If I do (85/15)/(15/85) the number is about 33, but is that correct Bayes Factor? 
 A: The Bayes factor for testing $H_0:\ p\le q$ when $X\sim\mathcal{N}(n,p)$ and $Y\sim\mathcal{N}(m,p)$ is$$\mathfrak{B}_{01}(x,y)=\dfrac{\int_{p\le q} f_X(x|p)f_Y(y|q)\pi(p,q)\text{d}p\text{d}q}{\int_{p\ge q} f_X(x|p)f_Y(y|q)\pi(p,q)\text{d}p\text{d}q}=\dfrac{\int_{p\le q} p^x(1-p)^{n-x}q^y(1-q)^{m-y}\text{d}p\text{d}q}{\int_{p\ge q} p^x(1-p)^{n-x}q^y(1-q)^{m-y}\text{d}p\text{d}q}$$
which happens to be
$$\mathfrak{B}_{01}(x,y)=\dfrac{\mathbb{P}(p\le q|x,y)}{\mathbb{P}(p\ge q|x,y)}$$
indeed. Hence, the R code in the question is correct.

Note that in the question "uninformed (1,1)" should be either
  "Uniform(0,1)" or "Beta(1,1)".

A: Indeed you are now computing the BF of theta1 > theta2 against its complement, theta1 ≤ theta2. Note that if it may be that theta1 and theta2 are 'exactly' equal (this depends on the context) you would need to perform an precise test instead of a one-sided test. The choice of the prior (i.e., where do you expect theta1 and theta2 to lie if they would be unequal) would then affect the results. For the BF for the precise test you could evaluate the posterior and prior at the nul value (also known as the Savage-Dickey density ratio), i.e.,
theta1_post = rbeta(100000,20+1,980+1)
theta2_post = rbeta(100000,30+1,1970+1)
theta_dif_post = theta1_post-theta2_post
dens_post <- density(theta_dif_post, from=0, to=0, n=1)$y
theta1_prior = rbeta(100000,1,1)
theta2_prior = rbeta(100000,1,1)
theta_dif_prior = theta1_prior-theta2_prior
dens_prior <- density(theta_dif_prior, from=0, to=0, n=1)$y
BFprecise <- dens_post / dens_prior

For your data this results in a lot of evidence for equality of theta's (a BF of around 48). Note (again) that the choice of the prior matters more than for one-sided testing. The above computation is done when using uniform priors, which implies that before observing the data you assume that any value for the theta's in the interval (0,1) is equally likely a priori. Further note that the Savage-Dickey density ratio should be used with some care but here I think it is alright.
A: Thank you so much for making this post! I was looking to solve the same problem, and yours is the only one I found that satisfies both (1) compare two samples (2) The null hypothesis isn't just p1 = p2, but p1 < p2. Your approach makes sense to me. Below is a breakdown.
Setup.

*

*X: The data. X1, X2 ~ Bin(n1, p2), Bin(n2, p2), where n1= 1000, n2=2000, and the sampled data are: successes = [20,30] respectively.

*H0: p1 > p2, H1: p1 <= p2

*Assume uniform priors for p1, p2.

Author's procedure.

*

*Want: Compute Bayes factor $\frac{P(x | H1)}{P(x | H0)}$. Since we assume uniform priors, then P(H0) = P(H1). The bayes factor is the same as computing the posterior odds ratio.

*Want: Compute posterior odds ratio $\frac{P(H1 | x)}{P(H0 | x)}$. We note that the denominator is just 1 - numerator. So, let's just focus on the numerator.

*Want: Compute $P(H1 | x)$. We expand this out to be: $P(p1 > p2 | x1, x2)$

*The author approximates $P(p1 > p2 | x1, x2)$ by sampling the individual posteriors:

*

*Sample p1 from the individual posterior $P(p1 | x1)$, which is a Beta distribution (See 3.1 from these notes on conjugate priors), similarly p2 from its posterior

*Estimate $P(p1 > p2 | x1, x2)$ to be: out of the (p1,p2) pairs sampled from the individual posteriors, what fraction of them has p1 > p2?



Proof, and gaps.
So far, the entire procedure makes sense to me. The only thing we're missing is the proof that your sampling method is valid. I lack the background knowledge to do this, so can only offer an intuitive explanation: "P(p1 > p2 | x1, x2) = integral p1>p2 joint_pdf(p1,p2 | x1,x2) = integral p1>p2 f(p1|x2) f(p2|x2)". This last term, corresponds to your sampling approach. -- A request for someone who has learned about MCMC to verify the intuition please :)
