# Extending a Hierarchical Beta-Binomial Model to account for higher-level groups

I searched all over but was unable to find an answer to this question. Please forgive me if I missed something obvious.

In order to analyze an experiment, I recently implemented a Hierarchical Beta-Binomial Model like the one detailed in this blog post. This model fit my needs because I had 5 variants and wanted to measure the proportion of "successes" (i.e., clicks on a particular link) in each one.

(In brief, this model works by assuming that each variant's success rate can be modeled as a Binomial Random Variable and that the rate for each Binomial Random Variable is drawn from a common Beta distribution.)

Now, this is awesome for my simple experiment; however, let's say I'm also interested in analyzing this experiment for two separate populations: Wizards and Warlocks.

I could theoretically split my experiment into ten variants -- WizA, WarA, WizB, WarB, etc. -- but I had a few questions:

• Would that even be valid?

• Is there a better way? I suppose I could implement a version of a multilevel regression model, but I was wondering if there was a simple extension either of the conceptual model or of the code in that blog post that would elegantly handle this experimental use case: must I immediately and necessarily look into something more complex?

(Also: apologies for my naivete here -- I took basic stats in college but never learned much about these Bayesian approaches. I'm hoping to refine my understanding since this seems much more elegant than running a bunch of t-tests and applying manual Bonferroni corrections, which is all I ever learned how to do.)

• To clarify, is the idea to find the best variant (highest proportion of clicks) of the now 10 buckets? – ilanman Aug 16 '16 at 17:18
• Ah, good question: the goal is to find the best variant for Wizards and the best variant for Warlocks, even if those variants are different. Building two completely separate hierarchical beta-binomial models for each population seems weird since -- at least conceptually -- it seems like you might want to make your model aware of "all" the information possible. – drodd Aug 17 '16 at 2:07

The most obvious option is a multilevel logistic regression with $$\text{logit}(\pi_i) = u_i + \beta \times 1\{ i \text{ is warlock} \}$$ that was already mentioned as an option, but I thought I'd spell out the model in something other than R code. In this you can e.g. assume that $u_1 \sim N(\mu, \sigma^2)$ or that $\text{inverse-logit}(u_i) \sim \text{Beta}(a,b)$ with the latter option being attractive, because you can write the likelihood also in closed form as a beta-binomial model (see Kuss, O. "Statistical methods for meta‐analyses including information from studies without any events—add nothing to nothing and succeed nevertheless." Statistics in medicine 34.7 (2015): 1097-1116.).