hypothesis_soln.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hypothesis Testing\n",
    "==================\n",
    "\n",
    "Copyright 2016 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "import numpy\n",
    "import scipy.stats\n",
    "\n",
    "import matplotlib.pyplot as pyplot\n",
    "\n",
    "import first\n",
    "\n",
    "# some nicer colors from http://colorbrewer2.org/\n",
    "COLOR1 = '#7fc97f'\n",
    "COLOR2 = '#beaed4'\n",
    "COLOR3 = '#fdc086'\n",
    "COLOR4 = '#ffff99'\n",
    "COLOR5 = '#386cb0'\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part One"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you observe an apparent difference between two groups and you want to check whether it might be due to chance.\n",
    "\n",
    "As an example, we'll look at differences between first babies and others.  The `first` module provides code to read data from the National Survey of Family Growth (NSFG)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "live, firsts, others = first.MakeFrames()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll look at a couple of variables, including pregnancy length and birth weight.  The effect size we'll consider is the difference in the means.\n",
    "\n",
    "Other examples might include a correlation between variables or a coefficient in a linear regression.  The number that quantifies the size of the effect is called the \"test statistic\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TestStatistic(data):\n",
    "    group1, group2 = data\n",
    "    test_stat = abs(group1.mean() - group2.mean())\n",
    "    return test_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the first example, I extract the pregnancy length for first babies and others.  The results are pandas Series objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group1 = firsts.prglngth\n",
    "group2 = others.prglngth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual difference in the means is 0.078 weeks, which is only 13 hours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "actual = TestStatistic((group1, group2))\n",
    "actual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null hypothesis is that there is no difference between the groups.  We can model that by forming a pooled sample that includes first babies and others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, m = len(group1), len(group2)\n",
    "pool = numpy.hstack((group1, group2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can simulate the null hypothesis by shuffling the pool and dividing it into two groups, using the same sizes as the actual sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RunModel():\n",
    "    numpy.random.shuffle(pool)\n",
    "    data = pool[:n], pool[n:]\n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of running the model is two NumPy arrays with the shuffled pregnancy lengths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RunModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we compute the same test statistic using the simulated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TestStatistic(RunModel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the model 1000 times and compute the test statistic, we can see how much the test statistic varies under the null hypothesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_stats = numpy.array([TestStatistic(RunModel()) for i in range(1000)])\n",
    "test_stats.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the sampling distribution of the test statistic under the null hypothesis, with the actual difference in means indicated by a gray line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyplot.axvline(actual, linewidth=3, color='0.8')\n",
    "pyplot.hist(test_stats, color=COLOR5)\n",
    "pyplot.xlabel('difference in means')\n",
    "pyplot.ylabel('count')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value is the probability that the test statistic under the null hypothesis exceeds the actual value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvalue = sum(test_stats >= actual) / len(test_stats)\n",
    "pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the result is about 15%, which means that even if there is no difference between the groups, it is plausible that we could see a sample difference as big as 0.078 weeks.\n",
    "\n",
    "We conclude that the apparent effect might be due to chance, so we are not confident that it would appear in the general population, or in another sample from the same population.\n",
    "\n",
    "STOP HERE\n",
    "---------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Part Two\n",
    "========\n",
    "\n",
    "We can take the pieces from the previous section and organize them in a class that represents the structure of a hypothesis test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HypothesisTest(object):\n",
    "    \"\"\"Represents a hypothesis test.\"\"\"\n",
    "\n",
    "    def __init__(self, data):\n",
    "        \"\"\"Initializes.\n",
    "\n",
    "        data: data in whatever form is relevant\n",
    "        \"\"\"\n",
    "        self.data = data\n",
    "        self.MakeModel()\n",
    "        self.actual = self.TestStatistic(data)\n",
    "        self.test_stats = None\n",
    "\n",
    "    def PValue(self, iters=1000):\n",
    "        \"\"\"Computes the distribution of the test statistic and p-value.\n",
    "\n",
    "        iters: number of iterations\n",
    "\n",
    "        returns: float p-value\n",
    "        \"\"\"\n",
    "        self.test_stats = numpy.array([self.TestStatistic(self.RunModel()) \n",
    "                                       for _ in range(iters)])\n",
    "\n",
    "        count = sum(self.test_stats >= self.actual)\n",
    "        return count / iters\n",
    "\n",
    "    def MaxTestStat(self):\n",
    "        \"\"\"Returns the largest test statistic seen during simulations.\n",
    "        \"\"\"\n",
    "        return max(self.test_stats)\n",
    "\n",
    "    def PlotHist(self, label=None):\n",
    "        \"\"\"Draws a Cdf with vertical lines at the observed test stat.\n",
    "        \"\"\"\n",
    "        pyplot.hist(ht.test_stats, color=COLOR4)\n",
    "        pyplot.axvline(self.actual, linewidth=3, color='0.8')\n",
    "        pyplot.xlabel('test statistic')\n",
    "        pyplot.ylabel('count')\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        raise UnimplementedMethodException()\n",
    "\n",
    "    def MakeModel(self):\n",
    "        \"\"\"Build a model of the null hypothesis.\n",
    "        \"\"\"\n",
    "        pass\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        raise UnimplementedMethodException()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`HypothesisTest` is an abstract parent class that encodes the template.  Child classes fill in the missing methods.  For example, here's the test from the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DiffMeansPermute(HypothesisTest):\n",
    "    \"\"\"Tests a difference in means by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        group1, group2 = data\n",
    "        test_stat = abs(group1.mean() - group2.mean())\n",
    "        return test_stat\n",
    "\n",
    "    def MakeModel(self):\n",
    "        \"\"\"Build a model of the null hypothesis.\n",
    "        \"\"\"\n",
    "        group1, group2 = self.data\n",
    "        self.n, self.m = len(group1), len(group2)\n",
    "        self.pool = numpy.hstack((group1, group2))\n",
    "\n",
    "    def RunModel(self):\n",
    "        \"\"\"Run the model of the null hypothesis.\n",
    "\n",
    "        returns: simulated data\n",
    "        \"\"\"\n",
    "        numpy.random.shuffle(self.pool)\n",
    "        data = self.pool[:self.n], self.pool[self.n:]\n",
    "        return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the test by instantiating a DiffMeansPermute object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = (firsts.prglngth, others.prglngth)\n",
    "ht = DiffMeansPermute(data)\n",
    "p_value = ht.PValue(iters=1000)\n",
    "print('\\nmeans permute pregnancy length')\n",
    "print('p-value =', p_value)\n",
    "print('actual =', ht.actual)\n",
    "print('ts max =', ht.MaxTestStat())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot the sampling distribution of the test statistic under the null hypothesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ht.PlotHist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Difference in standard deviation\n",
    "\n",
    "**Exercize 1**: Write a class named `DiffStdPermute` that extends `DiffMeansPermute` and overrides `TestStatistic` to compute the difference in standard deviations.  Is the difference in standard deviations statistically significant?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here\n",
    "\n",
    "class DiffStdPermute(DiffMeansPermute):\n",
    "    \"\"\"Tests a difference in means by permutation.\"\"\"\n",
    "\n",
    "    def TestStatistic(self, data):\n",
    "        \"\"\"Computes the test statistic.\n",
    "\n",
    "        data: data in whatever form is relevant        \n",
    "        \"\"\"\n",
    "        group1, group2 = data\n",
    "        test_stat = abs(group1.std() - group2.std())\n",
    "        return test_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the code to test your solution to the previous exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = (firsts.prglngth, others.prglngth)\n",
    "ht = DiffStdPermute(data)\n",
    "p_value = ht.PValue(iters=1000)\n",
    "print('\\nstd permute pregnancy length')\n",
    "print('p-value =', p_value)\n",
    "print('actual =', ht.actual)\n",
    "print('ts max =', ht.MaxTestStat())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Difference in birth weights\n",
    "\n",
    "Now let's run DiffMeansPermute again to see if there is a difference in birth weight between first babies and others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = (firsts.totalwgt_lb.dropna(), others.totalwgt_lb.dropna())\n",
    "ht = DiffMeansPermute(data)\n",
    "p_value = ht.PValue(iters=1000)\n",
    "print('\\nmeans permute birthweight')\n",
    "print('p-value =', p_value)\n",
    "print('actual =', ht.actual)\n",
    "print('ts max =', ht.MaxTestStat())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, after 1000 attempts, we never see a sample difference as big as the observed difference, so we conclude that the apparent effect is unlikely under the null hypothesis.  Under normal circumstances, we can also make the inference that the apparent effect is unlikely to be caused by random sampling.\n",
    "\n",
    "One final note: in this case I would report that the p-value is less than 1/1000 or less than 0.001.  I would not report p=0, because  the apparent effect is not impossible under the null hypothesis; just unlikely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part Three\n",
    "\n",
    "In this section, we'll explore the dangers of p-hacking by running multiple tests until we find one that's statistically significant.\n",
    "\n",
    "Suppose we want to compare IQs for two groups of people.  And suppose that, in fact, the two groups are statistically identical; that is, their IQs are drawn from a normal distribution with mean 100 and standard deviation 15.\n",
    "\n",
    "I'll use `numpy.random.normal` to generate fake data I might get from running such an experiment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group1 = numpy.random.normal(100, 15, size=100)\n",
    "group2 = numpy.random.normal(100, 15, size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We expect the mean in both groups to be near 100, but just by random chance, it might be higher or lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group1.mean(), group2.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use DiffMeansPermute to compute the p-value for this fake data, which is the probability that we would see a difference between the groups as big as what we saw, just by chance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = (group1, group2)\n",
    "ht = DiffMeansPermute(data)\n",
    "p_value = ht.PValue(iters=1000)\n",
    "p_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's check the p-value.  If it's less than 0.05, the result is statistically significant, and we can publish it.  Otherwise, we can try again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if p_value < 0.05:\n",
    "    print('Congratulations!  Publish it!')\n",
    "else:\n",
    "    print('Too bad!  Please try again.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can probably see where this is going.  If we play this game over and over (or if many researchers play it in parallel), the false positive rate can be as high as 100%.\n",
    "\n",
    "To see this more clearly, let's simulate 100 researchers playing this game.  I'll take the code we have so far and wrap it in a function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_a_test(sample_size=100):\n",
    "    \"\"\"Generate random data and run a hypothesis test on it.\n",
    "\n",
    "    sample_size: integer\n",
    "\n",
    "    returns: p-value\n",
    "    \"\"\"\n",
    "    group1 = numpy.random.normal(100, 15, size=sample_size)\n",
    "    group2 = numpy.random.normal(100, 15, size=sample_size)\n",
    "    data = (group1, group2)\n",
    "    ht = DiffMeansPermute(data)\n",
    "    p_value = ht.PValue(iters=200)\n",
    "    return p_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run that function 100 times and save the p-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_experiments = 100\n",
    "p_values = numpy.array([run_a_test() for i in range(num_experiments)])\n",
    "sum(p_values < 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On average, we expect to get a false positive about 5 times out of 100.  To see why, let's plot the histogram of the p-values we got."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = numpy.linspace(0, 1, 21)\n",
    "bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyplot.hist(p_values, bins, color=COLOR5)\n",
    "pyplot.axvline(0.05, linewidth=3, color='0.8')\n",
    "pyplot.xlabel('p-value')\n",
    "pyplot.ylabel('count')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution of p-values is uniform from 0 to 1.  So it falls below 5% about 5% of the time.\n",
    "\n",
    "**Exercise:** If the threshold for statistical signficance is 5%, the probability of a false positive is 5%.  You might hope that things would get better with larger sample sizes, but they don't.  Run this experiment again with a larger sample size, and see for yourself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part four\n",
    "\n",
    "In the previous section, we computed the false positive rate, which is the probability of seeing a \"statistically significant\" result, even if there is no statistical difference between groups.\n",
    "\n",
    "Now let's ask the complementary question: if there really is a difference between groups, what is the chance of seeing a \"statistically significant\" result?\n",
    "\n",
    "The answer to this question is called the \"power\" of the test.  It depends on the sample size (unlike the false positive rate), and it also depends on how big the actual difference is.\n",
    "\n",
    "We can estimate the power of a test by running simulations similar to the ones in the previous section.  Here's a version of `run_a_test` that takes the actual difference between groups as a parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_a_test2(actual_diff, sample_size=100):\n",
    "    \"\"\"Generate random data and run a hypothesis test on it.\n",
    "\n",
    "    actual_diff: The actual difference between groups.\n",
    "    sample_size: integer\n",
    "\n",
    "    returns: p-value\n",
    "    \"\"\"\n",
    "    group1 = numpy.random.normal(100, 15, \n",
    "                                 size=sample_size)\n",
    "    group2 = numpy.random.normal(100 + actual_diff, 15, \n",
    "                                 size=sample_size)\n",
    "    data = (group1, group2)\n",
    "    ht = DiffMeansPermute(data)\n",
    "    p_value = ht.PValue(iters=200)\n",
    "    return p_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run it 100 times with an actual difference of 5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_values = numpy.array([run_a_test2(5) for i in range(100)])\n",
    "sum(p_values < 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With sample size 100 and an actual difference of 5, the power of the test is approximately 65%.  That means if we ran this hypothetical experiment 100 times, we'd expect a statistically significant result about 65 times.\n",
    "\n",
    "That's pretty good, but it also means we would NOT get a statistically significant result about 35 times, which is a lot.\n",
    "\n",
    "Again, let's look at the distribution of p-values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyplot.hist(p_values, bins, color=COLOR5)\n",
    "pyplot.axvline(0.05, linewidth=3, color='0.8')\n",
    "pyplot.xlabel('p-value')\n",
    "pyplot.ylabel('count')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the point of this example: if you get a negative result (no statistical significance), that is not always strong evidence that there is no difference between the groups.  It is also possible that the power of the test was too low; that is, that it was unlikely to produce a positive result, even if there is a difference between the groups.\n",
    "\n",
    "**Exercise:** Assuming that the actual difference between the groups is 5, what sample size is needed to get the power of the test up to 80%?  What if the actual difference is 2, what sample size do we need to get to 80%?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}