incubator/generalized-abcde-tests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import janitor\n",
    "\n",
    "from scipy.stats import ttest_ind\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Bayesian A/B\\[/C/D/E...\\] tests\n",
    "\n",
    "## Abstract\n",
    "\n",
    "\"If you use Bayesian methods, there should be nothing stopping you from doing comparisons of more than two groups.\" With this pedagogical article, I hope to make the case that mastering a language that allows us to construct arbitrary models enables us to go beyond simple case/control (A/B) testing, giving us a tool to flexibly analyze data originating from complex experimental designs. \n",
    "\n",
    "## Introduction\n",
    "\n",
    "### Observations on A/B Testing\n",
    "\n",
    "Bayesian A/B testing has been gaining popularity over the years. I believe that there are two conceptual and two practical reasons for this.\n",
    "\n",
    "#### Conceptual Reasons\n",
    "\n",
    "Firstly, A/B tests are conceptually simple. I have a baseline (group \"A\"), and compare it to a perturbed setting (group \"B\"), and do an analysis comparing some summary statistic of both groups.\n",
    "\n",
    "Secondly, established methods exist for doing A/B analyses, the most common being the t-test. Under certain assumptions, say, equal variance, with a t-distributed likelihood function, having carefully determined the degrees of freedom... At the risk of revealing a spoiler of this article, how often do we actually encounter these assumptions? Can we go further and do A/B\\[/C/D/E\\] testing? What about breaking these rigid assumptions?\n",
    "\n",
    "#### Practical Reasons\n",
    "\n",
    "Firstly, we have probabilistic programming languages that allow us to express the classical t-test in a Bayesian fashion. \n",
    "\n",
    "Secondly, we have a wealth of freely available teaching material online, enabling a new generation of data scientists to conduct Bayesian A/B tests.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### On Probabilistic Programming Languages (PPLs)\n",
    "\n",
    "The first practical reason touched on PPLs. PPLs, such as PyMC3 (#cite), Stan (#cite) and Edward (#cite), allow us to express a statistical model in terms of (a) probability distributions describing our beliefs about a particular parameter, essentially assigning \"credibility points\" to parameter values and/or ranges, and (b) mathematical functions (a.k.a. transformations) that link what statisticians call \"independent\" or \"predictor\" variables to one or more \"target(s)\" or \"output(s)\". (#diagram) \n",
    "\n",
    "By use of simulation methods (essentially Markov Chain Monte Carlo (MCMC) sampling), PPLs provide a level of abstraction away from the difficult mathematics underlying Bayesian statistical inference. Without PPLs, Bayesian statistical analysis would often involve intractable integrals; without PPLs, in order to make these integrals tractable, we would usually have to limit ourselves to models that take advantage of distribution conjugacy, such as the Beta-Bernoulli/Binomial model, or assume Gaussian-distributed parameters and likelihoods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Going Beyond A/B Testing\n",
    "\n",
    "In retrospect, the idea that one should be able go beyond A/B testing should be obvious if you know the mechanics behind doing a Bayesian A/B test. Therefore, nothing I write below should be \"new knowledge\". However, I do not think this point is sufficiently popularized at conferences. Here's just a sampling of titles from conferences, based on a keyword search of \"A/B testing conferences\" (given the personalization that happens on search engines, your replication mileage may vary):\n",
    "\n",
    "- https://odsc.com/training/portfolio/percent-change-estimation-ab-testing\n",
    "- https://odsc.com/training/portfolio/mastering-ab-testing-from-design-to-analysis\n",
    "- https://github.com/datapythonista/odsc-bayesian-ab-testing-talk/blob/master/Bayesian%20AB%20testing.ipynb\n",
    "- https://www.youtube.com/watch?v=8DMU9h4dwi8\n",
    "- https://www.youtube.com/watch?v=PSqtcNZDj4A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An Example of Bayesian A/B testing\n",
    "\n",
    "This section is intended for readers who are not familiar with the mechanics of a Bayesian A/B test. Here, we will use a modified example from the Bayesian Estimation Supersedes the t-Test paper (#cite).\n",
    "\n",
    "Let us look at a fictional scenario, where we have data measuring the effect of a drug on IQ. The dose has been set; here we are essentially doing a gold-standard, two-arm, randomized controlled trial (RCT), where one group has been given a placebo, and the other has been given a drug, and participants are randomized between the groups.\n",
    "\n",
    "To get ourselves warmed up, let us look at the data that comes from this trial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drug = [\n",
    "    99.0,\n",
    "    110.0,\n",
    "    107.0,\n",
    "    104.0,\n",
    "    103.0,\n",
    "    105.0,\n",
    "    105.0,\n",
    "    110.0,\n",
    "    99.0,\n",
    "    109.0,\n",
    "    100.0,\n",
    "    102.0,\n",
    "    104.0,\n",
    "    104.0,\n",
    "    100.0,\n",
    "    104.0,\n",
    "    101.0,\n",
    "    104.0,\n",
    "    101.0,\n",
    "    100.0,\n",
    "    109.0,\n",
    "    104.0,\n",
    "    105.0,\n",
    "    112.0,\n",
    "    97.0,\n",
    "    106.0,\n",
    "    103.0,\n",
    "    101.0,\n",
    "    101.0,\n",
    "    104.0,\n",
    "    96.0,\n",
    "    102.0,\n",
    "    101.0,\n",
    "    100.0,\n",
    "    92.0,\n",
    "    108.0,\n",
    "    97.0,\n",
    "    106.0,\n",
    "    96.0,\n",
    "    90.0,\n",
    "    109.0,\n",
    "    108.0,\n",
    "    105.0,\n",
    "    104.0,\n",
    "    110.0,\n",
    "    92.0,\n",
    "    100.0,\n",
    "]\n",
    "\n",
    "placebo = [\n",
    "    95.0,\n",
    "    105.0,\n",
    "    103.0,\n",
    "    99.0,\n",
    "    104.0,\n",
    "    98.0,\n",
    "    103.0,\n",
    "    104.0,\n",
    "    102.0,\n",
    "    91.0,\n",
    "    97.0,\n",
    "    101.0,\n",
    "    100.0,\n",
    "    113.0,\n",
    "    98.0,\n",
    "    102.0,\n",
    "    100.0,\n",
    "    105.0,\n",
    "    97.0,\n",
    "    94.0,\n",
    "    104.0,\n",
    "    92.0,\n",
    "    98.0,\n",
    "    105.0,\n",
    "    106.0,\n",
    "    101.0,\n",
    "    106.0,\n",
    "    105.0,\n",
    "    101.0,\n",
    "    105.0,\n",
    "    102.0,\n",
    "    95.0,\n",
    "    91.0,\n",
    "    99.0,\n",
    "    96.0,\n",
    "    102.0,\n",
    "    94.0,\n",
    "    93.0,\n",
    "    99.0,\n",
    "    99.0,\n",
    "    113.0,\n",
    "    96.0,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It usually is instructive to compare the two distributions of data. To avoid binning bias, we will plot the empirical cumulative distributions of both groups on the same plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ECDF(data):\n",
    "    x = np.sort(data)\n",
    "    y = np.cumsum(x) / np.sum(x)\n",
    "\n",
    "    return x, y\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 4))\n",
    "ax.plot(*ECDF(drug), label=\"drug\")\n",
    "ax.plot(*ECDF(placebo), label=\"placebo\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequentist t-test\n",
    "\n",
    "If we were to do a two-sample t-test, working with the relaxed assumption that variances are not necessarily equal between the two groups, we would get the following result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ttest_ind(drug, placebo, equal_var=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the significance level of $\\alpha = 0.05$, we obtained a p-value of 0.02, which means that at a \"significance level\" of 0.05, we reject the null hypothesis that drug treatment did not increase the IQ of participants.\n",
    "\n",
    "If that sounded confusing, fret not because you are not alone: this author usually needs to articulate that sentence a few times to properly comprehend that we have only rejected the null hypothesis, but not accepted the \"alternate\" hypothesis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian estimation\n",
    "\n",
    "According to Kruschke's paper, we can move beyond the t-test, and flexibly create a model that explicitly encodes our assumptions about the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_vals = np.concatenate([drug, placebo])\n",
    "labels = [\"drug\"] * len(drug) + [\"placebo\"] * len(placebo)\n",
    "\n",
    "data = pd.DataFrame([y_vals, labels]).T\n",
    "data.columns = [\"IQ\", \"treatment\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this model, we explicitly state that we are modelling the data with t-distributions, with a shared degrees of freedom nuisance parameter and separately modelled means and variances. (This also encodes the not-necessarily-equal variances assumption into the model.) Finally, we also assume that the means of both likelihoods are normally distributed, that without seeing the data, we assume a very flat prior centered on the average human IQ, and that the variances are Cauchy-distributed, but can only take positive values (i.e. positive support), as expressed by the HalfCauchy distributions used below.\n",
    "\n",
    "Using PyMC3, a PPL that allows statistical models to be expressed using Python syntax, the aforementioned model can be implemented as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as kruschke_model:\n",
    "\n",
    "    mu_drug = pm.Normal(\"mu_drug\", mu=100, sd=100 ** 2)\n",
    "    mu_placebo = pm.Normal(\"mu_placebo\", mu=100, sd=100 ** 2)\n",
    "    sigma_drug = pm.HalfCauchy(\"sigma_drug\", beta=100)\n",
    "    sigma_placebo = pm.HalfCauchy(\"sigma_placebo\", beta=100)\n",
    "    nu = pm.Exponential(\"nu\", lam=1 / 29) + 1\n",
    "\n",
    "    drug_like = pm.StudentT(\"drug\", nu=nu, mu=mu_drug, sd=sigma_drug, observed=drug)\n",
    "    placebo_like = pm.StudentT(\n",
    "        \"placebo\", nu=nu, mu=mu_placebo, sd=sigma_placebo, observed=placebo\n",
    "    )\n",
    "    diff_means = pm.Deterministic(\"diff_means\", mu_drug - mu_placebo)\n",
    "    pooled_sd = pm.Deterministic(\n",
    "        \"pooled_sd\", np.sqrt(np.power(sigma_drug, 2) + np.power(sigma_placebo, 2) / 2)\n",
    "    )\n",
    "    effect_size = pm.Deterministic(\"effect_size\", diff_means / pooled_sd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using PyMC3, inference is performed by using two lines of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with kruschke_model:\n",
    "    kruschke_trace = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As is good practice in Bayesian statistical modelling, before arriving at conclusions, we use visual diagnostics to check that sampling has converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.traceplot(kruschke_trace, varnames=[\"mu_drug\", \"mu_placebo\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The visual diagnostic for a good sampling run is a \"hairy caterpillar\" type of MCMC trace, as shown on the right hand side panels.\n",
    "\n",
    "Non-hairy-caterpillar MCMC sampling runs is usually an indication of mathematical issues in the model, in which case a reparameterization or respecification of priors might be in order. Fixing these issues, however, is not a central focus of this article, hence the examples in this article will not show those issues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should compare the two posterior distributions on the same x-axis scale, in order to see the effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.forestplot(kruschke_trace, varnames=[\"mu_drug\", \"mu_placebo\"], vline=101.7)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.forestplot(kruschke_trace, varnames=[\"diff_means\"], vline=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the difference in means, we can see that the posterior distribution is centered at approximately 2.5 IQ units, with the 95% highest posterior density ranging from approximately 0.5 to 4.5 IQ units, and that the posterior does not overlap with a difference in 0 IQ units. Hence we might be tempted to conclude that the drug does indeed increase IQ!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While that statement is true under the data provided and the model that we have specified, looking at the two posterior means is generally not enough, though. What we are really interested in is whether there is a large enough effect size to justify the treatment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.forestplot(kruschke_trace, varnames=[\"effect_size\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By computing the posterior distribution over the effect size, we realize that the effect size is most probably of a \"medium\" size, but that it could also credibly be very large or very small. I would be most comfortable concluding that we need more experimentation to know what is going on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beyond two-group Bayesian estimation\n",
    "\n",
    "Bayesian estimation has the following advantages:\n",
    "\n",
    "- We have more assumptions that are explicitly stated in the model, and fewer assumptions that are left implicitly stated.\n",
    "- We only have to compare posterior distribution beliefs. Multiple hypothesis correction procedures are probably irrelevant.\n",
    "- Instead of making convoluted claims about 95% confidence intervals (#cite), we instead can state that there is a 95% probability that the true effect size is within the posterior highest density interval, *given the data*.\n",
    "\n",
    "The Bayesian estimation paper only dealt with two-group comparisons. That said, if we think carefully about Bayesian estimation logic, there should be nothing, in principle, that stops us from including a third, fourth, fifth (and more) groups into the modelling procedure. \n",
    "\n",
    "Probabilistic programming allows us to add that third, fourth... up to nth group easily.\n",
    "\n",
    "To do this, let's add in a 3rd and 4th group, which is drug + a non-pharmaceutical intervention group, and a non-pharmaceutical intervention-only group. Those familiar with experiment design will recognize this as a 2x2 matrix, with ±intervention and ±drug on each axis respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intervention_drug = np.random.randint(low=100, high=110, size=23)\n",
    "intervention_only = np.random.randint(low=110, high=120, size=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_new = (\n",
    "    data.append(\n",
    "        pd.DataFrame({\"IQ\": intervention_drug, \"treatment\": \"intervention+drug\",})\n",
    "    )\n",
    "    .append(pd.DataFrame({\"IQ\": intervention_only, \"treatment\": \"intervention\"}))\n",
    "    .label_encode(\"treatment\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With very compact syntax,  we can express a similar model. The key here is to use the `shape` keyword argument for each distribution object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as multigroup_model:\n",
    "\n",
    "    mu = pm.Normal(\"mu\", mu=100, sd=10, shape=(4,))\n",
    "    sd = pm.HalfCauchy(\"sd\", beta=10, shape=(4,))\n",
    "    nu = pm.Exponential(\"nu\", 1 / 29.0) + 1\n",
    "\n",
    "    like = pm.StudentT(\n",
    "        \"like\",\n",
    "        mu=mu[data_new[\"treatment_enc\"].values],\n",
    "        sd=sd[data_new[\"treatment_enc\"].values],\n",
    "        nu=nu,\n",
    "        observed=data_new[\"IQ\"].values,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can sample from the posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with multigroup_model:\n",
    "    trace_mg = pm.sample(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, the bare minimum that we should do is a sampling convergence check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.traceplot(trace_mg, varnames=[\"mu\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can examine the posterior distributions on the same scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = dict(sorted(zip(data_new[\"treatment_enc\"], data_new[\"treatment\"]))).values()\n",
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.forestplot(trace_mg, varnames=[\"mu\"], ylabels=labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For brevity here, we will ignore the computation of effect sizes. \n",
    "\n",
    "Comparing all four groups together, we can see very clearly that while the drug might increase IQ relative to a placebo, when compared to an intervention, its effect is less pronounced, and when combined with an intervention, the drug actually reduces the effect of an intervention! \n",
    "\n",
    "There are two points to be made here. \n",
    "\n",
    "Firstly, this is an example of an A/B/C/D test, where there are four treatments (one of them being the control group). There is no reason to stop at A/B testing.\n",
    "\n",
    "Secondly, our inferences can be made without requiring the usual ANOVA+post-hoc t-testing. All we have to do is to look at the computed posterior distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalizing\n",
    "\n",
    "### Beyond t-distributed data\n",
    "\n",
    "The t-distribution is the generalization of the Normal and Cauchy distributions, and hence is a flexible distribution for modelling purposes. Here, we assumed that our data were t-distributed. However, even this assumption can be changed for a different distribution, with different parameters, and hence different priors placed on those parameters. For example, if our measured variable were count data, then we could change the distribution to a discrete one, such as the poisson or negative binomial, or even zero-inflated versions of them.\n",
    "\n",
    "### Going hierarchical\n",
    "\n",
    "While not discussed here, hierarchical models are very easy to implement in PyMC3. Using the example in this notebook, the major change is that instead of specifying a single number for the `mu` and `sd` parameters across all groups, we specify a probability distribution that it is drawn from, with its own \"hyper-parameters\". More specifically:\n",
    "\n",
    "```python\n",
    "with pm.Model() as multigroup_model:\n",
    "    \n",
    "    # Group-level hyperparameters\n",
    "    group_mu = pm.Normal('group_mu', mu=100, sd=10)\n",
    "    group_sd = pm.Exponential('group_sd', lam=5)\n",
    "    \n",
    "    # Individual group parameters\n",
    "    mu = pm.Normal('mu', mu=group_mu, sd=group_sd, shape=(4,))\n",
    "    sd = pm.HalfCauchy('sd', beta=group_beta, shape=(4,))\n",
    "    nu = pm.Exponential('nu', lam=1/29.) + 1\n",
    "    like = pm.StudentT('like', mu=mu[groups], sd=sd[groups], nu=nu, observed=y_data)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "My hope is that this notebook article highlights for you the fact that we need not stop at two group, case/control (A/B) comparisons. With probabilistic programming, we can flexibly and easily specify multi-group comparisons.\n",
    "\n",
    "For further reading, I would suggest checking out the [Bayesian Analysis Recipes](https://github.com/ericmjl/bayesian-analysis-recipes) repository that I have been building over the years. In the [README](https://github.com/ericmjl/bayesian-analysis-recipes/blob/master/README.md), there are further resources for learning available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bayesian",
   "language": "python",
   "name": "bayesian"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}