notebooks/a. BEST.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Estimation Supersedes the T-Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Problem\n",
    "\n",
    "Several statistical inference procedures involve the comparison of two groups. We may be interested in whether one group is larger than another, or simply different from the other. We require a statistical model for this because true differences are usually accompanied by measurement or stochastic noise that prevent us from drawing conclusions simply from differences calculated from the observed data. \n",
    "\n",
    "The *de facto* standard for statistically comparing two (or more) samples is to use a statistical test. This involves expressing a null hypothesis, which typically claims that there is no difference between the groups, and using a chosen test statistic to determine whether the distribution of the observed data is plausible under the hypothesis. This rejection occurs when the calculated test statistic is higher than some pre-specified threshold value.\n",
    "\n",
    "Unfortunately, it is not easy to conduct hypothesis tests correctly, and their results are very easy to misinterpret. Setting up a statistical test involves several subjective choices (*e.g.* statistical test to use, null hypothesis to test, significance level) by the user that are rarely justified based on the problem or decision at hand, but rather, are usually based on traditional choices that are entirely arbitrary (Johnson 1999). The evidence that it provides to the user is indirect, incomplete, and typically overstates the evidence against the null hypothesis (Goodman 1999). \n",
    "\n",
    "A more informative and effective approach for comparing groups is one based on **estimation** rather than **testing**, and is driven by Bayesian probability rather than frequentist. That is, rather than testing whether two groups are different, we instead pursue an estimate of how different they are, which is fundamentally more informative. Moreover, we include an estimate of uncertainty associated with that difference which includes uncertainty due to our lack of knowledge of the model parameters (epistemic uncertainty) and uncertainty due to the inherent stochasticity of the system (aleatory uncertainty)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Drug trial evaluation\n",
    "\n",
    "To illustrate how this Bayesian estimation approach works in practice, we will use a fictitious example from Kruschke (2012) concerning the evaluation of a clinical trial for drug evaluation. The trial aims to evaluate the efficacy of a \"smart drug\" that is supposed to increase intelligence by comparing IQ scores of individuals in a treatment arm (those receiving the drug) to those in a control arm (those recieving a placebo). There are 47 individuals and 42 individuals in the treatment and control arms, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "sns.set(color_codes=True)\n",
    "\n",
    "np.random.seed(20090425)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "drug = (101,100,102,104,102,97,105,105,98,101,100,123,105,103,100,95,102,106,\n",
    "        109,102,82,102,100,102,102,101,102,102,103,103,97,97,103,101,97,104,\n",
    "        96,103,124,101,101,100,101,101,104,100,101)\n",
    "placebo = (99,101,100,101,102,100,97,101,104,101,102,102,100,105,88,101,100,\n",
    "           104,100,100,100,101,102,103,97,101,101,100,101,99,101,100,100,\n",
    "           101,100,99,101,100,102,99,100,99)\n",
    "\n",
    "y1 = np.array(drug)\n",
    "y2 = np.array(placebo)\n",
    "y = pd.DataFrame(dict(value=np.r_[y1, y2], group=np.r_[['drug']*len(drug), ['placebo']*len(placebo)]))\n",
    "\n",
    "y.hist('value', by='group');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in a Bayesian approach to inference is to specify the full probability model that corresponds to the problem. For this example, Kruschke chooses a Student-t distribution to describe the distributions of the scores in each group. This choice adds robustness to the analysis, as a T distribution is less sensitive to outlier observations, relative to a normal distribution. The three-parameter Student-t distribution allows for the specification of a mean $\\mu$, a precision (inverse-variance) $\\lambda$ and a degrees-of-freedom parameter $\\nu$:\n",
    "\n",
    "$$ f(x|\\mu,\\lambda,\\nu) =\n",
    "           \\frac{\\Gamma(\\frac{\\nu + 1}{2})}{\\Gamma(\\frac{\\nu}{2})}\n",
    "           \\left(\\frac{\\lambda}{\\pi\\nu}\\right)^{\\frac{1}{2}}\n",
    "           \\left[1+\\frac{\\lambda(x-\\mu)^2}{\\nu}\\right]^{-\\frac{\\nu+1}{2}}$$\n",
    "           \n",
    "the degrees-of-freedom parameter essentially specifies the \"normality\" of the data, since larger values of $\\nu$ make the distribution converge to a normal distribution, while small values (close to zero) result in heavier tails.\n",
    "\n",
    "Thus, the likelihood functions of our model are specified as follows:\n",
    "\n",
    "$$y^{(treat)}_i \\sim T(\\nu, \\mu_1, \\sigma_1)$$\n",
    "\n",
    "$$y^{(placebo)}_i \\sim T(\\nu, \\mu_2, \\sigma_2)$$\n",
    "\n",
    "As a simplifying assumption, we will assume that the degree of normality $\\nu$ is the same for both groups. We will, of course, have separate parameters for the means $\\mu_k, k=1,2$ and standard deviations $\\sigma_k$.\n",
    "\n",
    "Since the means are real-valued, we will apply normal priors on them, and arbitrarily set the hyperparameters to the pooled empirical mean of the data and twice the pooled empirical standard deviation, which applies very diffuse information to these quantities (and importantly, does not favor one or the other *a priori*).\n",
    "\n",
    "$$\\mu_k \\sim N(\\bar{x}, 2s)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "μ_m = y.value.mean()\n",
    "μ_s = y.value.std() * 2\n",
    "\n",
    "with pm.Model() as model:\n",
    "    \n",
    "    group1_mean = pm.Normal('group1_mean', μ_m, sd=μ_s)\n",
    "    group2_mean = pm.Normal('group2_mean', μ_m, sd=μ_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The group standard deviations will be given a uniform prior over a plausible range of values for the variability of the outcome variable, IQ.\n",
    "\n",
    "In Kruschke's original model, he uses a very wide uniform prior for the group standard deviations, from the pooled empirical standard deviation divided by 1000 to the pooled standard deviation multiplied by 1000. This is a poor choice of prior, because very basic prior knowledge about measures of human coginition dictate that the variation cannot ever be as high as this upper bound. IQ is a standardized measure, and hence this constrains how variable a given population's IQ values can be. When you place such a wide uniform prior on these values, you are essentially giving a lot of prior weight on inadmissable values. In this example, there is little practical difference, but in general it is best to apply as much prior information that you have available to the parameterization of prior distributions. \n",
    "\n",
    "We will instead set the group standard deviations to have a $\\text{Uniform}(1,10)$ prior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "σ_low = 1\n",
    "σ_high = 10\n",
    "\n",
    "with model:\n",
    "    \n",
    "    group1_std = pm.Uniform('group1_std', lower=σ_low, upper=σ_high)\n",
    "    group2_std = pm.Uniform('group2_std', lower=σ_low, upper=σ_high)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We follow Kruschke by making the prior for $\\mu$ exponentially distributed with a mean of 30; this allocates high prior probability over the regions of the parameter that describe the range from normal to heavy-tailed data under the Student-T distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with model:\n",
    "    \n",
    "    ν = pm.Exponential('ν_minus_one', 1/29.) + 1\n",
    "\n",
    "sns.distplot(np.random.exponential(30, size=10000), kde=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since PyMC3 parameterizes the Student-T in terms of precision, rather than standard deviation, we must transform the standard deviations before specifying our likelihoods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with model:\n",
    "    \n",
    "    λ1 = group1_std**-2\n",
    "    λ2 = group2_std**-2\n",
    "\n",
    "    group1 = pm.StudentT('drug', nu=ν, mu=group1_mean, lam=λ1, observed=y1)\n",
    "    group2 = pm.StudentT('placebo', nu=ν, mu=group2_mean, lam=λ2, observed=y2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having fully specified our probabilistic model, we can turn our attention to calculating the comparisons of interest in order to evaluate the effect of the drug. To this end, we can specify deterministic nodes in our model for the difference between the group means and the difference between the group standard deviations. Wrapping them in named `Deterministic` objects signals to PyMC that we wish to record the sampled values as part of the output.\n",
    "\n",
    "As a joint measure of the groups, we will also estimate the \"effect size\", which is the difference in means scaled by the pooled estimates of standard deviation. This quantity can be harder to interpret, since it is no longer in the same units as our data, but the quantity is a function of all four estimated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with model:\n",
    "\n",
    "    diff_of_means = pm.Deterministic('difference of means', group1_mean - group2_mean)\n",
    "    diff_of_stds = pm.Deterministic('difference of stds', group1_std - group2_std)\n",
    "    effect_size = pm.Deterministic('effect size', \n",
    "                                   diff_of_means / pm.sqrt((group1_std**2 + group2_std**2) / 2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can fit the model and evaluate its output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with model:\n",
    "    trace = pm.sample(2000,  njobs=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the stochastic parameters of the model. PyMC's `plot_posterior` function replicates the informative histograms portrayed in Kruschke (2012). These summarize the posterior distributions of the parameters, and present a 95% credible interval and the posterior mean. The plots below are constructed with the final 1000 samples from each of the 2 chains, pooled together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.plot_posterior(trace[1000:], \n",
    "                  varnames=['group1_mean', 'group2_mean', 'group1_std', 'group2_std', 'ν_minus_one'],\n",
    "                  color='#87ceeb');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the group differences, we can conclude that there are meaningful differences between the two groups for all three measures. For these comparisons, it is useful to use zero as a reference value (`ref_val`); providing this reference value yields cumulative probabilities for the posterior distribution on either side of the value. Thus, for the difference in means, 98.9% of the posterior probability is greater than zero, which suggests the group means are credibly different. The effect size and differences in standard deviation are similarly positive.\n",
    "\n",
    "These estimates suggest that the \"smart drug\" increased both the expected scores, but also the  variability in scores across the sample. So, this does not rule out the possibility that some recipients may be adversely affected by the drug at the same time others benefit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.plot_posterior(trace[1000:], \n",
    "                  varnames=['difference of means', 'difference of stds', 'effect size'],\n",
    "                  ref_val=0,\n",
    "                  color='#87ceeb');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When `forestplot` is called on a trace with more than one chain, it also plots the potential scale reduction parameter, which is used to reveal evidence for lack of convergence; values near one, as we have here, suggest that the model has converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.forestplot(trace[1000:], varnames=['difference of means', 'difference of stds', 'effect size'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.plots.summary(trace[1000:], \n",
    "                 varnames=['difference of means', 'difference of stds', 'effect size'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1.\tGoodman SN. Toward evidence-based medical statistics. 1: The P value fallacy. Annals of Internal Medicine. 1999;130(12):995-1004. doi:10.7326/0003-4819-130-12-199906150-00008.\n",
    "2.\tJohnson D. The insignificance of statistical significance testing. Journal of Wildlife Management. 1999;63(3):763-772.\n",
    "3.\tKruschke JK. Bayesian estimation supersedes the t test. J Exp Psychol Gen. 2013;142(2):573-603. doi:10.1037/a0029146."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}