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}$).


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:

  1. 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}$)
  2. 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 with dpar(0.5,1) priors but not with dgamma(2,20) priors. Thanks 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 calling 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)

A visualization of the 4 datasets included in the Excel file

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
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)

    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)

    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)
    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:


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)


"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))


model in Try1.bug
data in Try1.r
compile, nchains(2)
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.

  • $\begingroup$ You have skipped crucial aspects in your description. After you uniformly draw an $x$, what happens next? Something must, for otherwise you have absolutely no information about any of the $p_x$! Then, what do you observe? Do you observe $x$ and, say, an independent draw from $o_x$? Or do you just observe the draw from $o_x$ and do not record $x$? To estimate all these density distributions, surely you must repeat the process many, many times--how many? Or, by "estimate" do you assume you know all the $p_x$ and you would like to calculate the distribution of the outcome of your experiment? $\endgroup$
    – whuber
    Aug 28, 2013 at 17:59
  • $\begingroup$ Sorry, perhaps the phrase "$x$ is uniformly drawn from 1...N" was misleading. Given a large number of existing observations about various $o_x$, I want to determine the distribution of $p_x$. It's somewhat analogous to sampling 50 people from a population, observing their heights, and then asking for the distribution of the height for a person "randomly drawn from the population". I just used uniform since the popuation in this case was enumerated, and I see it's confusing. (Continuing with the analogy, in this scenario, we don't get to measure the people, just observe a few binary responses) $\endgroup$ Aug 28, 2013 at 20:20
  • $\begingroup$ I've added an example to the end of the question that hopefully makes it clear what I'm trying to do. $\endgroup$ Aug 28, 2013 at 20:53
  • $\begingroup$ Thanks. The question remains confusing, though. In the example you are observing an ordered pair of values: times rolled and successes. Exactly what "distribution" are you inquiring about, then? Precisely what property of the big bag of dice would it correspond to? $\endgroup$
    – whuber
    Aug 28, 2013 at 20:54
  • 1
    $\begingroup$ The pair of values (rolls & successes) conveys information about a latent variable, p. The distribution I want is the distribution of p's in the bag. Of course, if we don't assume a class of distributions that it belongs to, there are arbitrarily many possibilities, but if we say that p is for example beta-distributed, then it comes down to just selecting the best fitting parameters for this class of distributions. If it helps, now consider instead of dice with p=1/4,1/6,1/8, that each object has a p from Beta(2,2)... or Beta (4,4), etc. $\endgroup$ Aug 29, 2013 at 3:03

1 Answer 1


Since the time for the bounty expired and I received no answers, I'll post the answer I was able to come up with, though my limited experience with Bayesian inference suggests that this should be taken with a healthy dose of skepticism.

I) Setup

  1. I downloaded & installed JAGS 3.4 and R 3.0.1
  2. I installed the packages rjags and coda by initiating R and using install.packages(pkgname)

II) Model & Data - Used the model & data files already detailed in the question. To answer question #1, I added one additional observation onto the data with all four variables as 0.

III) Answering questions

  1. I ran JAGS on the model/data (open the command line, go the the directory with your files, and type > jags-terminal Command.cmd . It ran and output a few files
  2. In R, I used the following commands:
    • library("rjags") to load the installed package (and its required package coda)
    • setwd() to get to the directory where the outputfiles were
    • results=read.coda("STEMchain1.txt","STEMindex.txt")
  3. To answer the first question:
    • As a PDF plot, "plot(results[,3*N])"
    • As quantiles, "quantile(results[,3*N],c(0.025,0.25,0.5,0.75,0.975))"
    • Where N is the number of observations, and the last observation corresponds to the position of the "all 0's" observation. (1 to n is for variable p_A, n+1 to 2n is for p_B, and 2n+1 to 3n is for delta)
  4. To answer the second question, same as above, but change 3*N -> 2*N+y

I'm not sure this is the correct way to get the answer, or whether a more complex model would produce better results, particularly in the case of correlation, but hopefully eventually someone more experienced will chime in...


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.