Why Hierarchical? : I've tried researching this problem, and from what I understand, this is a "hierarchical" problem, because you are making observations about observations from a population, rather than making direct observations from that population. Reference: http://www.econ.umn.edu/~bajari/iosp07/rossi1.pdf
Why Bayesian? : Also, I've tagged it as Bayesian because an asymptotic/frequentist solution may exist for an "experimental design" wherein each "cell" is assigned sufficient observations, but for practical purposes real-world/non-experimental data sets (or at least mine) are sparsely populated. There's a lot of aggregate data, but individual cells may be blank or have just a few observations.
The model in abstract:
Let U be a population of units ${u_1, u_2, u_3 ... u_N}$ to each of which we can apply a treatment, $T$, of either $A$ or $B$, and from each of which we observe i.i.d. observations of either 1 or 0 a.k.a. successes and failures. Let $p_{iT}$ for $i \in \{1...N\}$ be the probability that an observation from object $i$ under treatment $T$ results in a success. Note that $p_{iA}$ may be correlated with $p_{iB}$.
To make the analysis feasible, we (a) assume that the distributions $p_A$ and $p_B$ are each an instance of a specific family of distributions, such as the beta distribution, (b) and select some prior distributions for hyperparameters.
Example of the model
I have a really big bag of Magic 8 Balls. Each 8-ball, when shaken, can reveal "Yes", or "No". Also, I can shake them with the ball right-side-up, or upside-down (assuming our Magic 8 Balls work upside-down...). The orientation of the ball may completely change the probability of resulting in a "Yes" or a "No" (in other words, initially you have no belief that $p_{iA}$ is correlated with $p_{iB}$).
Questions:
Someone has randomly sampled a number, $n$, of units from the population, and for each unit has taken and recorded an arbitrary number of observations under treatment $A$ and an arbitrary number of observations under treatment $B$. (In practice, in our data set, most units will have observations only under one treatment)
Given this data, I need to answer the following questions:
- If I take a new unit $u_x$ randomly from the population, how can I calculate (analytically or stochastically) the joint posterior distribution of $p_{xA}$ and $p_{xB}$? (Primarily so that we can determine the expected difference in proportions, $\Delta=p_{xA}-p_{xB}$)
- For a specific unit $u_y$, $y \in \{1,2,3...,n\}$, with observations of $s_y$ successes and $f_y$ failures, how can I calculate (analytically or stochastically) the joint posterior distribution for $p_{yA}$ and $p_{yB}$, again to build a distribution $\Delta_y$ of the difference in proportions $p_{yA}-p_{yB}$
Bonus question: Although we really do expect $p_A$ and $p_B$ to be very correlated, we are not explicitly modeling that. In the likely case of a stochastic solution, I believe this would cause some samplers, including Gibbs, the be less effective at exploring the posterior distribution. Is this the case, and if so, should we use a different sampler, somehow model the correlation as a separate variable and transform the $p$ distributions to make them uncorrelated, or just run the sampler longer?
Answer criteria
I'm looking for an answer that:
Has code, using preferably Python/PyMC, or barring that, JAGS, that I am able to run
Is able to handle an input of thousands of units
Given enough units and samples, is able to output distributions for $p_A$, $p_B$, and $\Delta$ as an answer to Question 1, that can be shown to match the underlying population distributions (to be checked against excel sheets provided in the "challenge datasets" section)
Given enough units and samples, is able to output the right distributions for $p_A$, $p_B$, and $\Delta$ (I'll use the excel sheets provided in the "challenge datasets" section to check) as an answer to Question 2, and provides some rationale as to why these distributions are correct
If the answer is similar to the last JAGS model I posted, explain why it works withThanks to Martyn Plummer on the JAGS forum for catching the error in my JAGS model. In trying to set a prior of Gamma(Shape=2,Scale=20), I was callingdpar(0.5,1)
priors but not withdgamma(2,20)
priors.dgamma(2,20)
which actually set a prior of Gamma(Shape=2,InverseScale=20) = Gamma(Shape=2,Scale=0.05).
Challenge datasets
I've generated a few sample datasets in Excel, with a few different possible scenarios, changing the tightness of the p distributions, the correlation between them, and making it easy to change other inputs. https://docs.google.com/file/d/0B_rPBjs4Cp0zLVBybU1nVnd0ZFU/edit?usp=sharing (~8Mb)
My attempted/partial solution(s) to date
1) I downloaded and installed Python 2.7 and PyMC 2.2. Initially, I got an incorrect model to run, but when I tried to reformulate the model, the extension freezes. By adding/removing code, I've determined the code that triggers the freeze is mc.Binomial(...) , though this function did work in the first model, so I assume there is something wrong with how I've specified the model.
import pymc as mc
import numpy as np
import scipy.stats as stats
from __future__ import division
cases=[0,0]
for case in range(2):
if case==0:
# Taken from the sample datasets excel sheet, Focused Correlated p's
c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightCorr.tsv", unpack=True)
if case==1:
# Taken from the sample datasets excel sheet, Focused Uncorrelated p's
c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightUncorr.tsv", unpack=True)
scale=20.0
alpha_A=mc.Uniform("alpha_A", 1,scale)
beta_A=mc.Uniform("beta_A", 1,scale)
alpha_B=mc.Uniform("alpha_B", 1,scale)
beta_B=mc.Uniform("beta_B", 1,scale)
p_A=mc.Beta("p_A",alpha=alpha_A,beta=beta_A)
p_B=mc.Beta("p_B",alpha=alpha_B,beta=beta_B)
@mc.deterministic
def delta(p_A=p_A,p_B=p_B):
return p_A-p_B
obs_n_A=mc.DiscreteUniform("obs_n_A",lower=0,upper=20,observed=True, value=n_A_arr)
obs_n_B=mc.DiscreteUniform("obs_n_B",lower=0,upper=20,observed=True, value=n_B_arr)
obs_c_A=mc.Binomial("obs_c_A",n=obs_n_A,p=p_A, observed=True, value=c_A_arr)
obs_c_B=mc.Binomial("obs_c_B",n=obs_n_B,p=p_B, observed=True, value=c_B_arr)
model = mc.Model([alpha_A,beta_A,alpha_B,beta_B,p_A,p_B,delta,obs_n_A,obs_n_B,obs_c_A,obs_c_B])
cases[case] = mc.MCMC(model)
cases[case].sample(24000, 12000, 2)
lift_samples = cases[case].trace('delta')[:]
ax = plt.subplot(211+case)
figsize(12.5,5)
plt.title("Posterior distributions of lift from 0 to T")
plt.hist(lift_samples, histtype='stepfilled', bins=30, alpha=0.8,
label="posterior of lift", color="#7A68A6", normed=True)
plt.vlines(0, 0, 4, color="k", linestyles="--", lw=1)
plt.xlim([-1, 1])
2) I downloaded and installed JAGS 3.4. After getting a correction on my priors from the JAGS forum, I now have this model, which successfully runs:
Model
var alpha_A, beta_A, alpha_B, beta_B, p_A[N], p_B[N], delta[N], n_A[N], n_B[N], c_A[N], c_B[N];
model {
for (i in 1:N) {
c_A[i] ~ dbin(p_A[i],n_A[i])
c_B[i] ~ dbin(p_B[i],n_B[i])
p_A[i] ~ dbeta(alpha_A,beta_A)
p_B[i] ~ dbeta(alpha_B,beta_B)
delta[i] <- p_A[i]-p_B[i]
}
alpha_A ~ dgamma(1,0.05)
alpha_B ~ dgamma(1,0.05)
beta_A ~ dgamma(1,0.05)
beta_B ~ dgamma(1,0.05)
}
Data
"N" <- 60
"c_A" <- structure(c(0,6,0,3,0,8,0,4,0,6,1,5,0,5,0,7,0,3,0,7,0,4,0,5,0,4,0,5,0,4,0,2,0,4,0,5,0,8,2,7,0,6,0,3,0,3,0,8,0,4,0,4,2,6,0,7,0,3,0,1))
"c_B" <- structure(c(5,0,2,2,2,0,2,0,2,0,0,0,5,0,4,0,3,1,2,0,2,0,2,0,0,0,3,0,6,0,4,1,5,0,2,0,6,0,1,0,2,0,4,0,4,1,1,0,3,0,5,0,0,0,5,0,2,0,7,1))
"n_A" <- structure(c(0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3))
"n_B" <- structure(c(9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3))
Control
model in Try1.bug
data in Try1.r
compile, nchains(2)
initialize
update 400
monitor set p_A, thin(3)
monitor set p_B, thin(3)
monitor set delta, thin(3)
update 1000
coda *, stem(Try1)
The actual application for anyone who would rather pick apart the model :)
On the web, typical A/B testing considers the impact to conversion rate from a single page or unit of content, with possible variations. Typical solutions include a classical significance test against the null hypotheses or two equal proportions, or more recently analytical Bayesian solutions leveraging the beta distribution as a conjugate prior.
Instead of this single-unit-of-content approach, which incidentally would require a lot of visitors to each unit of I am interested in testing, we want to compare variations in a process that generates multiple units of content (not an unusual scenario really...). So on the whole, the units/pages produced by process A or B have a lot of visits/data, but each individual unit may only have a few observations.