tumbleweed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## There is only one test\n",
    "\n",
    "Copyright 2018 Allen Downey\n",
    "\n",
    "License: http://creativecommons.org/licenses/by/4.0/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The tumbleweed problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example that uses the \"[There is only one test](http://allendowney.blogspot.com/2016/06/there-is-still-only-one-test.html)\" framework to solve a problem from the [2019 AP® Statistics Exam](https://secure-media.collegeboard.org/apc/ap19-frq-statistics.pdf)\n",
    "\n",
    "> Tumbleweed, commonly found in the western United States, is the dried structure of certain plants that are\n",
    "blown by the wind. Kochia, a type of plant that turns into tumbleweed at the end of the summer, is a problem for\n",
    "farmers because it takes nutrients away from soil that would otherwise go to more beneficial plants. Scientists\n",
    "are concerned that kochia plants are becoming resistant to the most commonly used herbicide, glyphosate. \n",
    "\n",
    "> In 2014, 19.7 percent of 61 randomly selected kochia plants were resistant to glyphosate. In 2017, 38.5 percent of 52 randomly selected kochia plants were resistant to glyphosate. Do the data provide convincing statistical\n",
    "evidence, at the level of $\\alpha = 0.05$, that there has been an increase in the proportion of all kochia plants that are resistant to glyphosate? \n",
    "\n",
    "Because the problem reports the data in terms of percentages, we have to work backwards to get the actual number of plants in each group that were resistant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.017000000000001"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.197 * 61"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.02"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.385 * 52"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like there were 12 resistant plants in the first group and 20 in the second group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "32"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k1 = 12\n",
    "k2 = 20\n",
    "k = k1 + k2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here are the sample sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "113"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n1 = 61\n",
    "n2 = 52\n",
    "n = n1 + n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step toward hypothesis testing is to define a test statistic that measures the size of the effect.  In this case, the hypothetical effect is an increase in the fraction of resistant plants, that is, a difference in proportions.\n",
    "\n",
    "I'll write a function that takes the data as a parameter and computes the difference in proportions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_stat(data):\n",
    "    \"\"\"Compute the test statistic.\"\"\"\n",
    "    k1, n1, k2, n2 = data\n",
    "    p1 = k1/n1\n",
    "    p2 = k2/n2\n",
    "    return p2 - p1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I pass the data to the function as a tuple of 4 numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(12, 61, 20, 52)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = k1, n1, k2, n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the observed difference in proportions (and check it against the numbers in the problem)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1878940731399748"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "actual = test_stat(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.188"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.385 - 0.197"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need a null hypothesis.  The scientists think there might be an increase in the fraction of resistant plants, so the null hypothesis is that there is no increase.\n",
    "\n",
    "Under that assumption, the two groups have the same proportion of resistant plants, which we can estimate by pooling the data and computing the proportion among all of the observed plants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2831858407079646"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = k / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the fundamental question in hypothesis testing is this: under the null hypothesis, what would be the chance of seeing an apparent effect as big as `actual`?\n",
    "\n",
    "To answer that question, we simulate the null hypothesis and generate \"fake data\", then compute the test statistic of the fake data.\n",
    "\n",
    "I'll use the following function, which simulates flipping `n` coins, each with probability `p` of landing heads, and returns the total number of heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flip(n, p):\n",
    "    \"\"\"Flip n coins with probability p.\n",
    "    \n",
    "    n: number of coins\n",
    "    p: probability of heads\n",
    "    \n",
    "    return: total number of heads\n",
    "    \"\"\"\n",
    "    return np.sum(np.random.random(n) < p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, I'll simulate the two groups and compute the number of resistant plants in each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip(n1, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flip(n2, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can wrap that in a function that returns the results in the same format as the original data, a tuple of 4 numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_model(data):\n",
    "    \"\"\"Simulate the null hypothesis.\"\"\"\n",
    "    k1, n1, k2, n2 = data\n",
    "    p = (k1 + k2) / (n1 + n2)\n",
    "    k1 = flip(n1, p)\n",
    "    k2 = flip(n2, p)\n",
    "    return k1, n1, k2, n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Returning the results in that format makes it possible to pass the \"fake data\" directly to `test_stat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0012610340479193016"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_stat(run_model(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's the result of one simulated experiment.\n",
    "\n",
    "Now we run the experiment 1000 times and collect the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.1056116 ,  0.00693569, -0.00945776, -0.11349306, -0.01229508,\n",
       "        0.10025221,  0.07534678, -0.02301387,  0.00693569, -0.15983607,\n",
       "        0.00126103, -0.16551072,  0.02332913, -0.07786885, -0.15479193,\n",
       "        0.00189155, -0.04224464, -0.07786885,  0.10813367, -0.02868852,\n",
       "        0.20428752,  0.10025221,  0.0759773 ,  0.03972257, -0.04791929,\n",
       "       -0.0037831 , -0.08858764, -0.00662043, -0.08921816,  0.01544767,\n",
       "       -0.06147541,  0.07534678, -0.09993695,  0.11380832, -0.04224464,\n",
       "       -0.12988651, -0.1056116 , -0.00094578,  0.1223203 , -0.18757881,\n",
       "        0.03184111, -0.02301387, -0.04508197,  0.02332913, -0.1056116 ,\n",
       "        0.05390921,  0.02900378, -0.08354351,  0.00409836, -0.06147541,\n",
       "        0.11948298, -0.0501261 ,  0.09741488, -0.03940731,  0.08102144,\n",
       "       -0.00662043,  0.11380832, -0.08575032, -0.08070618, -0.11916772,\n",
       "       -0.01450189, -0.0258512 , -0.0721942 ,  0.110971  ,  0.10025221,\n",
       "       -0.04791929,  0.02332913, -0.04224464,  0.05958386,  0.00977301,\n",
       "       -0.08858764,  0.03972257,  0.02332913, -0.03089533, -0.11633039,\n",
       "       -0.03940731,  0.00409836, -0.07786885, -0.02017654,  0.00693569,\n",
       "        0.02332913,  0.08669609, -0.12988651, -0.1406053 , -0.12137453,\n",
       "        0.04823455,  0.06179067,  0.11664565,  0.02332913,  0.06746532,\n",
       "        0.07030265,  0.07818411,  0.1223203 ,  0.15226986, -0.05296343,\n",
       "        0.16078184, -0.0372005 ,  0.04823455,  0.1223203 , -0.01796974,\n",
       "        0.02616646, -0.09142497, -0.1056116 , -0.05580076,  0.01544767,\n",
       "        0.0425599 , -0.09993695,  0.01765448, -0.04508197,  0.03184111,\n",
       "        0.1223203 , -0.03373266,  0.02332913, -0.11916772,  0.02112232,\n",
       "       -0.11349306,  0.07030265,  0.08102144,  0.06746532,  0.07250946,\n",
       "       -0.15983607, -0.04791929,  0.00977301,  0.06462799, -0.12704918,\n",
       "       -0.00945776,  0.01544767, -0.06431274,  0.14659521,  0.07250946,\n",
       "       -0.04224464, -0.08070618,  0.06462799,  0.04539723, -0.04508197,\n",
       "        0.08606557,  0.00977301,  0.05674653, -0.03940731, -0.05296343,\n",
       "       -0.0721942 ,  0.08669609, -0.05580076,  0.20996217,  0.00189155,\n",
       "       -0.03436318, -0.08070618,  0.07534678, -0.29728878,  0.05107188,\n",
       "        0.02332913, -0.01513241, -0.15479193,  0.02112232,  0.07250946,\n",
       "       -0.11349306,  0.02900378,  0.18505675,  0.03688525,  0.02900378,\n",
       "        0.08385876, -0.0942623 ,  0.01261034,  0.05390921, -0.07503153,\n",
       "        0.09741488, -0.11349306,  0.0425599 , -0.08070618,  0.03751576,\n",
       "        0.02616646, -0.01229508,  0.19861286, -0.02805801,  0.05107188,\n",
       "        0.12452711,  0.10592686, -0.10277427, -0.07503153, -0.03656999,\n",
       "       -0.06431274, -0.02805801,  0.08102144, -0.11916772, -0.0258512 ,\n",
       "       -0.03940731,  0.25409836, -0.01450189,  0.20208071,  0.04823455,\n",
       "        0.09457755,  0.12515763,  0.06462799,  0.07250946, -0.03656999,\n",
       "       -0.01733922, -0.00094578, -0.08070618,  0.04539723,  0.07313997,\n",
       "       -0.03089533,  0.05107188, -0.03152585,  0.02900378,  0.00693569,\n",
       "        0.09741488,  0.02112232,  0.00977301,  0.01261034, -0.11412358,\n",
       "        0.22068096, -0.01229508, -0.11349306, -0.00094578,  0.03467844,\n",
       "        0.23139975, -0.05075662, -0.0037831 , -0.12988651, -0.0942623 ,\n",
       "        0.10025221,  0.09741488,  0.10876419, -0.07503153,  0.08953342,\n",
       "       -0.02868852,  0.05327869,  0.07534678,  0.05107188,  0.0425599 ,\n",
       "        0.02900378, -0.10277427, -0.03152585,  0.00977301,  0.07250946,\n",
       "       -0.04224464,  0.06179067,  0.03467844, -0.10844893,  0.03467844,\n",
       "        0.17150063,  0.00693569,  0.03972257,  0.13871375,  0.00977301,\n",
       "       -0.13776797,  0.00977301, -0.12421185, -0.00094578, -0.06147541,\n",
       "       -0.09142497, -0.01229508, -0.11633039, -0.14407314,  0.03184111,\n",
       "       -0.03152585,  0.12168979,  0.04823455,  0.09457755,  0.03972257,\n",
       "       -0.0721942 ,  0.05044136,  0.09174023,  0.02616646,  0.09457755,\n",
       "        0.10245902, -0.10844893,  0.00126103,  0.14218159,  0.00409836,\n",
       "       -0.13556116, -0.03656999, -0.04003783,  0.10876419,  0.04035309,\n",
       "        0.110971  , -0.16834805, -0.05580076,  0.0889029 ,  0.13587642,\n",
       "        0.04539723, -0.00945776, -0.0258512 ,  0.10592686, -0.05580076,\n",
       "       -0.1406053 ,  0.07818411, -0.0942623 , -0.08070618, -0.0721942 ,\n",
       "       -0.08070618,  0.06462799,  0.06179067, -0.03089533, -0.01513241,\n",
       "        0.15226986, -0.00945776,  0.00977301,  0.10876419,  0.10592686,\n",
       "        0.15510719, -0.02080706, -0.00157629, -0.08354351, -0.09142497,\n",
       "       -0.08070618, -0.02017654, -0.06651955,  0.04539723,  0.03467844,\n",
       "       -0.03940731, -0.11633039,  0.10025221,  0.1223203 , -0.10277427,\n",
       "        0.01481715,  0.00977301,  0.05895334, -0.11412358,  0.04035309,\n",
       "       -0.04508197, -0.19546028,  0.00126103, -0.08070618,  0.01828499,\n",
       "        0.28121059,  0.01828499,  0.02616646, -0.08070618,  0.07534678,\n",
       "        0.02900378,  0.05107188,  0.01544767,  0.06746532,  0.06746532,\n",
       "       -0.04508197,  0.06179067,  0.06462799,  0.02616646, -0.01513241,\n",
       "       -0.01733922,  0.03184111,  0.00693569, -0.02868852,  0.06179067,\n",
       "        0.02616646,  0.06179067, -0.06431274,  0.07534678, -0.14911728,\n",
       "       -0.00094578, -0.15762926, -0.05580076, -0.11349306,  0.03530895,\n",
       "       -0.05359395, -0.15479193, -0.04224464,  0.00977301, -0.11916772,\n",
       "        0.0889029 ,  0.01544767, -0.09142497, -0.00441362, -0.12704918,\n",
       "        0.05895334,  0.00693569,  0.1223203 , -0.00945776, -0.09993695,\n",
       "       -0.04224464,  0.01828499, -0.11349306,  0.11664565, -0.08070618,\n",
       "        0.1223203 , -0.00945776, -0.02017654, -0.16551072, -0.11633039,\n",
       "       -0.10498108,  0.1443884 ,  0.00977301, -0.14691047, -0.03152585,\n",
       "        0.13240858, -0.14911728, -0.07503153, -0.07503153,  0.16078184,\n",
       "        0.02616646,  0.03184111,  0.08669609, -0.00662043,  0.09174023,\n",
       "       -0.00945776, -0.09709962, -0.07503153, -0.11916772, -0.00945776,\n",
       "       -0.1056116 ,  0.04539723,  0.06746532,  0.08669609, -0.12704918,\n",
       "        0.13303909, -0.06998739,  0.10025221, -0.08070618, -0.00945776,\n",
       "        0.00126103, -0.03940731,  0.04823455, -0.00945776,  0.07030265,\n",
       "        0.02616646, -0.00725095, -0.12988651,  0.03467844, -0.08575032,\n",
       "       -0.11633039,  0.01828499,  0.09741488, -0.03940731, -0.04508197,\n",
       "        0.03404792,  0.00126103,  0.00409836,  0.07818411,  0.04539723,\n",
       "       -0.11349306, -0.13556116,  0.14722573, -0.03940731, -0.02868852,\n",
       "       -0.05359395, -0.03940731,  0.10592686, -0.03940731,  0.12168979,\n",
       "        0.07534678,  0.01261034,  0.04539723, -0.01513241, -0.07786885,\n",
       "       -0.06651955, -0.20901639, -0.1056116 ,  0.06179067, -0.09205549,\n",
       "       -0.11349306,  0.13303909, -0.08070618,  0.10308953, -0.00094578,\n",
       "       -0.15479193, -0.12988651, -0.09993695, -0.01513241, -0.05296343,\n",
       "       -0.03656999, -0.04224464, -0.05863808, -0.02017654,  0.0425599 ,\n",
       "        0.04035309, -0.12704918, -0.00662043,  0.03688525, -0.06431274,\n",
       "        0.03751576, -0.0258512 , -0.00441362,  0.00693569,  0.19356873,\n",
       "       -0.05863808, -0.11065574, -0.06147541,  0.10025221, -0.01796974,\n",
       "        0.08102144,  0.02900378, -0.02017654, -0.01513241,  0.07030265,\n",
       "        0.00126103, -0.02301387,  0.0425599 ,  0.14155107,  0.01828499,\n",
       "        0.09457755, -0.04224464,  0.00693569, -0.17622951,  0.11380832,\n",
       "        0.08102144,  0.10308953,  0.04539723,  0.02332913,  0.02900378,\n",
       "       -0.03940731,  0.06746532, -0.08921816,  0.06462799, -0.0721942 ,\n",
       "        0.08385876,  0.17213115, -0.0721942 ,  0.01197982, -0.0258512 ,\n",
       "       -0.0964691 ,  0.11664565,  0.0425599 , -0.04224464, -0.12421185,\n",
       "        0.23203026, -0.04508197,  0.12799496,  0.01261034,  0.03467844,\n",
       "       -0.12421185,  0.06179067, -0.04508197,  0.1686633 , -0.01229508,\n",
       "       -0.07786885,  0.11380832, -0.0721942 ,  0.07818411,  0.02616646,\n",
       "       -0.08354351, -0.08354351, -0.03656999,  0.03688525, -0.07282472,\n",
       "       -0.01513241, -0.11349306, -0.00662043, -0.05296343,  0.07818411,\n",
       "        0.05107188, -0.03940731, -0.01229508,  0.02332913, -0.00945776,\n",
       "        0.14155107, -0.00945776, -0.0037831 ,  0.00409836,  0.22351828,\n",
       "       -0.01166456, -0.12704918,  0.02900378, -0.0258512 ,  0.07030265,\n",
       "        0.02616646, -0.19262295,  0.23486759, -0.00945776, -0.0721942 ,\n",
       "       -0.04508197,  0.00472888,  0.03972257, -0.01796974, -0.13272383,\n",
       "        0.06746532, -0.00094578, -0.18978562,  0.04823455, -0.03436318,\n",
       "       -0.00945776,  0.06179067, -0.16551072, -0.00441362,  0.04476671,\n",
       "       -0.06368222,  0.16298865,  0.17938209,  0.00977301,  0.02332913,\n",
       "        0.1443884 , -0.04508197,  0.01261034,  0.01261034,  0.06179067,\n",
       "       -0.05359395,  0.04823455, -0.11633039, -0.0258512 , -0.21469105,\n",
       "       -0.08858764,  0.02332913,  0.00693569, -0.15983607, -0.05359395,\n",
       "       -0.01229508,  0.0075662 ,  0.03972257,  0.09741488, -0.05296343,\n",
       "       -0.01229508, -0.02301387, -0.01229508, -0.07786885,  0.0889029 ,\n",
       "       -0.13556116, -0.03089533, -0.12200504, -0.1982976 , -0.04508197,\n",
       "       -0.07786885, -0.06368222,  0.05895334, -0.01796974, -0.06147541,\n",
       "        0.19577554,  0.1223203 ,  0.07818411, -0.00441362, -0.04508197,\n",
       "       -0.05580076, -0.0258512 , -0.00945776,  0.08102144, -0.06147541,\n",
       "       -0.04791929, -0.11916772, -0.02301387, -0.04791929,  0.00693569,\n",
       "        0.01261034,  0.04539723,  0.04539723,  0.11948298,  0.08669609,\n",
       "        0.10876419, -0.0258512 , -0.0037831 , -0.00157629, -0.0501261 ,\n",
       "        0.11664565,  0.16078184,  0.08102144,  0.00472888, -0.07503153,\n",
       "        0.02900378,  0.09237074,  0.11664565,  0.03751576,  0.0425599 ,\n",
       "       -0.16551072,  0.10813367,  0.10308953, -0.07786885, -0.16834805,\n",
       "       -0.0258512 ,  0.00693569,  0.01765448, -0.12704918, -0.02301387,\n",
       "        0.08102144,  0.12168979, -0.04791929, -0.09709962,  0.04319042,\n",
       "        0.10025221, -0.20680958,  0.09174023, -0.07786885, -0.02301387,\n",
       "        0.03688525,  0.03751576, -0.16834805, -0.02868852,  0.11664565,\n",
       "       -0.0372005 ,  0.09741488, -0.05296343, -0.15479193, -0.08921816,\n",
       "       -0.12421185,  0.16078184, -0.1056116 , -0.02868852,  0.11664565,\n",
       "       -0.07503153,  0.02395965, -0.00662043,  0.08102144, -0.19325347,\n",
       "        0.06179067, -0.04224464,  0.00693569, -0.09142497, -0.08291299,\n",
       "        0.03184111, -0.06715006, -0.04791929,  0.15510719,  0.09237074,\n",
       "       -0.0037831 , -0.07503153, -0.03436318,  0.01544767,  0.04823455,\n",
       "        0.15510719, -0.00662043, -0.05075662, -0.00662043,  0.04539723,\n",
       "        0.0425599 , -0.05296343,  0.09741488,  0.14155107,  0.01261034,\n",
       "        0.00126103, -0.07503153,  0.05611602,  0.06462799, -0.06431274,\n",
       "       -0.08354351, -0.03940731, -0.05863808, -0.04508197,  0.01765448,\n",
       "        0.08448928, -0.03152585,  0.02616646, -0.10781841, -0.02868852,\n",
       "        0.04319042,  0.01481715,  0.0075662 , -0.0942623 , -0.14911728,\n",
       "        0.02616646,  0.04035309, -0.04224464, -0.05359395, -0.06651955,\n",
       "       -0.0721942 , -0.09142497, -0.04728878,  0.13303909, -0.01796974,\n",
       "        0.06462799, -0.02080706, -0.1519546 , -0.04508197,  0.01828499,\n",
       "       -0.04508197,  0.08953342, -0.04224464,  0.10025221,  0.18505675,\n",
       "        0.09741488,  0.05107188,  0.05044136,  0.00977301, -0.11349306,\n",
       "        0.06179067,  0.07534678, -0.02017654, -0.04224464,  0.09741488,\n",
       "       -0.04224464, -0.07786885, -0.11916772, -0.03152585,  0.1686633 ,\n",
       "       -0.1740227 , -0.09709962, -0.08354351,  0.07313997,  0.11664565,\n",
       "       -0.07786885, -0.13556116, -0.01513241,  0.02900378,  0.00977301,\n",
       "       -0.04791929,  0.04823455, -0.0721942 ,  0.01765448,  0.09237074,\n",
       "       -0.00094578,  0.08669609,  0.04823455, -0.07786885, -0.00662043,\n",
       "       -0.12704918,  0.05327869,  0.16361917, -0.03940731,  0.11948298,\n",
       "        0.05327869,  0.04539723,  0.02900378, -0.2203657 ,  0.05895334,\n",
       "       -0.02017654, -0.01229508,  0.13303909, -0.02080706,  0.14092055,\n",
       "       -0.12704918,  0.07534678,  0.05390921,  0.02900378, -0.21752837,\n",
       "        0.18221942, -0.05296343, -0.04791929,  0.08385876,  0.04539723,\n",
       "       -0.01513241, -0.14911728,  0.11948298, -0.08638083, -0.01796974,\n",
       "       -0.02301387,  0.05611602,  0.19640605,  0.01261034, -0.13556116,\n",
       "        0.18505675,  0.09741488,  0.02332913,  0.00977301, -0.0721942 ,\n",
       "        0.02616646,  0.09457755,  0.03972257,  0.07030265,  0.0204918 ,\n",
       "        0.09741488,  0.09174023, -0.03436318, -0.06715006, -0.09142497,\n",
       "        0.13871375,  0.08953342, -0.00882724, -0.00157629, -0.01796974,\n",
       "       -0.05075662,  0.08165195, -0.09993695,  0.22351828,  0.11664565,\n",
       "       -0.00945776, -0.06431274, -0.06147541, -0.0372005 ,  0.07313997,\n",
       "       -0.06431274,  0.15794451,  0.00977301, -0.13839849, -0.11065574,\n",
       "       -0.0721942 ,  0.08385876, -0.05580076, -0.22320303,  0.08385876,\n",
       "       -0.18757881,  0.03688525,  0.12452711, -0.11065574,  0.07030265,\n",
       "        0.10308953, -0.15762926,  0.01261034, -0.01229508,  0.00126103,\n",
       "       -0.00945776,  0.01544767,  0.05044136, -0.05296343,  0.06967213,\n",
       "        0.10813367, -0.00094578, -0.02017654,  0.09457755,  0.02332913,\n",
       "       -0.00094578,  0.00977301, -0.13839849, -0.10781841,  0.10025221,\n",
       "        0.10592686,  0.03184111,  0.05107188, -0.0258512 ,  0.01765448,\n",
       "        0.13020177, -0.05359395, -0.03940731, -0.17906683, -0.12421185,\n",
       "        0.05390921,  0.0425599 , -0.08070618, -0.01796974, -0.05863808,\n",
       "        0.02332913,  0.07313997, -0.11065574, -0.11349306, -0.07786885,\n",
       "        0.10025221, -0.01229508, -0.2203657 ,  0.0204918 ,  0.14659521,\n",
       "       -0.08858764, -0.04508197, -0.15762926,  0.06746532, -0.08575032,\n",
       "        0.02616646,  0.02900378,  0.08669609, -0.21469105, -0.0501261 ,\n",
       "       -0.02301387,  0.07250946, -0.04003783,  0.03467844,  0.06179067,\n",
       "       -0.06147541, -0.00315259, -0.1056116 ,  0.00472888, -0.03436318,\n",
       "       -0.02364439,  0.08102144,  0.03467844, -0.06715006,  0.00189155,\n",
       "        0.13587642, -0.07503153,  0.08385876,  0.00126103,  0.08102144,\n",
       "       -0.12421185, -0.01513241, -0.10498108, -0.07503153, -0.18190416,\n",
       "        0.05390921, -0.03940731, -0.06431274, -0.07503153, -0.0037831 ,\n",
       "        0.03972257,  0.01261034,  0.07250946, -0.0037831 ,  0.07818411,\n",
       "        0.12452711,  0.17150063, -0.02301387, -0.02868852, -0.05296343,\n",
       "        0.03184111, -0.02868852, -0.03940731,  0.00126103, -0.09993695,\n",
       "       -0.01229508, -0.00094578, -0.06431274, -0.18474149,  0.06179067,\n",
       "       -0.14911728,  0.05390921,  0.13587642, -0.00945776,  0.10025221,\n",
       "        0.01261034, -0.08858764, -0.08575032,  0.09174023, -0.11633039,\n",
       "        0.08385876,  0.03184111, -0.06715006, -0.01229508, -0.0258512 ,\n",
       "        0.02332913, -0.00945776, -0.03656999,  0.07818411, -0.01229508,\n",
       "       -0.0258512 ,  0.01765448,  0.02900378,  0.00977301,  0.00126103])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_stat_dist = np.array([test_stat(run_model(data)) \n",
    "                           for i in range(1000)])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is the \"sampling distribution of the test statistic under the null hypothesis\".\n",
    "\n",
    "The mean of this distribution is close to zero, which is not suprising because it is based on the assumption that there is actually no difference between the groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0010378310214375767"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(test_stat_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the distribution of test statistics looks like.\n",
    "\n",
    "The vertical gray line shows the observed effect size, `actual`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XncW1Wdx/HPlx1ZZGlBCoUCgoqClanoDCooyLCIMI6IuAGiRQcGVFyqgwIqA8wIMzI6KIpSkEVkWCqgglWqoAVaZCsIFKhQW9uytWURafnNH+eEXtKb5D5LnuTJ832/XnklubnL7yQ3+eWcc++5igjMzMzqrdLpAMzMrDs5QZiZWSknCDMzK+UEYWZmpZwgzMyslBOEmZmVcoLoMEknSvpRfrylpKckrToE2z1M0g2F509J2maQ1v0lSd/Pj8dJCkmrDdK6h+w9qtvuppJ+I2mppNMrLjNH0p7tjm04Ke4b/Vh20PbR/pJ0vaSPdTKGoTRiE4Skt0j6naTFkh6XdKOkN3Yypoh4OCLWjYjlHdj2uhHxYLN5JO0uaW6Fdf17RAzKl6j+R7aD79FE4FFg/Yg4rv5FSedK+nq7Np6T7CsHYT0v/iFp97Jl+0vVfaPsh7jKPjqcNPrzVL8vSRovaaakZ/L9+KGKcUQmCEnrA1cB/wNsBGwOnAQ818m4esFg1RS60FbA3eEzS60f+vu9kLQGcCXwI2BDYDJwZZ7efhEx4m7ABODJJq9vC/wKeIz0r/ECYIPC63OAzwF3AE8D5wCbAj8DlgK/BDbM844DgvQPdB4wHziusK4TgR/Vzbtafn498DXgxrzea4FRhWU/Avwpx/nlHNeeDcq0MTAFWALcnNd7Q+H1AF6ZH+8L3J23+Wfgs8A6wLPAC8BT+TYmx38paQdeAnysQZkalf9c4OuF57sDc/Pj8/P2ns3b+3zJezQml+txYDbw8br39hLgvFyWWcCEJp/7PwC3AIvz/T8UYnwe+FuOY8+65SbWvf7Twn7yWdJ+shj4MbBWYbl3AbcBTwK/A3ZqENdvcpmfzus/uNXywBfyZ7cUuBfYA9g7x/h8Xs/tDbZXeVngcOCePO+DwJF5erP9pbZvrEXabx7LZbiF9D06GVgO/DUv962SfXRt4HTS/r8YuAFYu6Qsh1HYz0vWcy7wbeDqXIabgG0L874T+GPexreAacDHCq9/NJf/CeAXwFZ12zkKuB94qCS2cRT25bLvBLBX/ixUeP1hYO8h+a0cio102w1YP++Uk4F9yD/mhddfmXeMNYHRpC/ofxdenwNMzzvz5sBC4FbgDXmZXwEn1O0EF+UvzY7AIvKPDK0TxAPA9vkLcT1wan5th/zleQuwBvAN0pe3UYK4mPRjuQ7wurzTNUoQ84G35scbAjvnx7uTf7wLy52Yt3sgqUa6doMyNSr/i1+Gsm1Ql/RK3qNpwP+SfmzG53XvUYjtr6SEtypwCjC9wfuzEelL/mFgNeCQ/HzjsjhLll/p9Rz7zaQfxo1IPySfyK/tTNpv3pRjOzTPv2aD9b/4+bRaHngV8AgwpvCebVu/vzXYTp+WBfYj/aESsBvwDK33l9q+cSTwU+BluQx/R2rCg7Svf6zRe0D6Ub+e9P1blZTcV3rvqJYgHgd2yZ/7BcDF+bVRpD897wVWBz4NLKvFRdrnZwOvycseD/yubjvX5c++LHmNo3WC+DTws7rXr6LwJ6udtxHZxBQRS0g/rAF8D1gkaYqkTfPrsyPiuoh4LiIWAWeQdv6i/4mIBRHxZ+C3wE0R8YeIeA64nJQsik6KiKcj4k7gh6QfoCp+GBH3RcSzpB/4Wvvje0n/VG+IiL8BX8nlWUnu0P1n4Cs5hrtIybGR54EdJK0fEU9ExK0tYvx9RFwRES/kOMv0t/wNSRpL+hy/EBF/jYjbgO+TfuRrboiIayL1WZwPvL7B6vYD7o+I8yNiWURcRPrnuP8AwzwzIuZFxOOkH8Pa5/dx4LsRcVNELI+IyaQmzjdXXG+z5ZeTEsUOklaPiDkR8UDF9fZp2Yi4OiIeiGQaqZb71orbep5Us31lLsPM/N1sStIqpH/ux0bEn/Oyv8vfvf64LCJujohlpARR+4z2JTUrXhoRzwP/DfylsNyRwCkRcU9e9t+B8ZK2KsxzSkQ83uR7AfCopCdrN+ADhdfWJdVeihYD6/W5lP0wIhMEQP5QD4uILUj/qMeQdgAkbSLpYkl/lrSEVA0eVbeKBYXHz5Y8X7du/kcKj/+Ut1dFcYd8prDeMcV1RsQzpFpRmdGkfzj1MTTyz6Qvx58kTZP09y1ifKTF6/Xz9KX8zYwBHo+IpXXr3rzwvP79W6tBe/AYVn5P6tfVH40+v62A4+p+GMZS/X1puHxEzAY+Rfq3vjDvy5XW29dlJe0jaXo+0ONJ0n5T/11p5HxSs8zFkuZJ+g9Jq1dYbhSpxlg16bVS9TsWvHQ/3gr4ZuH9f5xUkyruM1W+G6MiYoPaDbiw8NpTpBaPovVJzWFtN2ITRFFE/JFUrXtdnnQK6d/4ThGxPvAh0gc/EGMLj7cktccPxHxgi9oTSWuT/o2VWUSqGtfHUCoibomIA4BNgCtINRdoUENpMr2oUfmfJjUx1LyiD+ueB2wkqfhvaktS81lfzSN94Yv6sq4q70HRI8DJxR+GiHhZrrkMePmIuDAi3kIqUwCnVY2z6rKS1gT+j9S8uWn+cbuGFd+VptuKiOcj4qSI2IHURPQuUr9aq2UfJTUdbtuqLNTtX5Lq969m5lPYbyWJl+7Hj5D6XIqfwdoR8bvCPH3dL+rNAnbK267ZKU9vuxGZICS9WtJxkrbIz8eSmjym51nWI2XuJyVtTuqQHqgvS3qZpNeSOvZ+PMD1XQrsL+kf8hENJ9EgieXmlcuAE3MMO5DarFciaQ1JH5T08lytXkJqdoBUS9pY0sv7EW+j8t8G7Ctpo/zl/VTdcguA0mPfI+IRUufsKZLWkrQTcASpmaCvrgG2l/QBSatJOpjUz3NVxeUbxtnA94BPSHqTknUk7VeX7Jqtv+Hykl4l6R35B/yvpBpt8TMcl5tpVtLHZdcgNUctApZJ2ofUqVqMueH+IuntknbMTaBLSE1OxW01+txfAH4AnCFpjKRVJf19jrne7cBr86Gia5FqRlVdnZd9T651HsNL/8B8B/hi3qeR9HJJB/Vh/VVcT3pPjpG0pqSj8/RfDfJ2So3IBEGqnr0JuEnS06TEcBdQO779JFIn4GLSTnLZIGxzGqlDayrwjYi4diAri4hZwL+SOp/nk8q0kMaH6h5Nqjr/hVRb+mGT1X8YmJOb1z5BqkHValoXAQ/manVfmokalf980pd4Dqn9uj5xngIcn7f32ZL1HkLq7JtH6vs5ISKu60NcAETEY6R/sMeRmuo+D7wrIh6tuIpzSO32T0q6osL2ZpD6Eb5F6gyfTepQbeREYHJe//taLL8mcCrpn/ZfSDXBL+XXfpLvH5NU1rdUednctHcMqYb5BKntfEqhjK32l1eQ/ugsIXXgTyM15wJ8E3ivpCcknVkS52eBO0lHPj1OquWs9HsWEfcBXyUdWXg/6WinSvJnfxDp/XgM2I50RGHt9cvzdi/O35W7SAe9DJrcv3ggqWb1JKnv5cA8ve2UmtWsXSSNAx4CVs8dWe3azrqkHWi7iHioXdsxs5FjpNYgeoKk/XOzzTqkduA7Sf/EzcwGzAlieDuA1LQyj1T9fX+4Smhmg8RNTGZmVso1CDMzKzWsB1YbNWpUjBs3rtNhmNkw8thj6XzSjTdudNpQ75s5c+ajETG61XzDOkGMGzeOGTNmdDoMMxtGzj33XAAOO+ywjsbRSZKajaTwIjcxmZlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JtSxCSxkr6taR7JM2SdGyefqLShXhuy7d9C8t8UdJsSfdK+sd2xWZmZq218zyIZaTrpt6ax7ifKak2DPN/RcQ3ijPnaxS8H3gt6UpOv5S0fb6WgZmZDbG21SAiYn7tWsZ53Ph7aH75xgNIFwt/Lg9XPZt0IXEzM+uAITmTOl8T4Q3ATcCuwNGSPgLMINUyniAlj+mFxeYy8OsBm3XMuElXd2S7c07dryPbtd7T9k7qfCGb/wM+FRFLgLNI15IdT7oS2um1WUsWX2moWUkTJc2QNGPRokVtitrMzNqaICStTkoOF0TEZQARsSAilufryn6PFc1Ic3npBcG3YMWF7V8UEWdHxISImDB6dMuxpszMrJ/aeRSTSNfpvScizihM36ww2z+RruMK6Vq2788X5t6adAGcm9sVn5mZNdfOPohdgQ8Dd0q6LU/7EnCIpPGk5qM5wJEAETFL0iXA3aQjoI7yEUxmZp3TtgQRETdQ3q9wTZNlTgZObldMZmZWnc+kNjOzUk4QZmZWalhfUc6slU6di2DWC1yDMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqXaliAkjZX0a0n3SJol6dg8fSNJ10m6P99vmKdL0pmSZku6Q9LO7YrNzMxaa2cNYhlwXES8BngzcJSkHYBJwNSI2A6Ymp8D7ANsl28TgbPaGJuZmbXQtgQREfMj4tb8eClwD7A5cAAwOc82GTgwPz4AOC+S6cAGkjZrV3xmZtbckPRBSBoHvAG4Cdg0IuZDSiLAJnm2zYFHCovNzdPq1zVR0gxJMxYtWtTOsM3MRrS2JwhJ6wL/B3wqIpY0m7VkWqw0IeLsiJgQERNGjx49WGGamVmdlglC0rGS1s+dyOdIulXSXlVWLml1UnK4ICIuy5MX1JqO8v3CPH0uMLaw+BbAvKoFMTOzwVWlBvHR/M9/L2A0cDhwaquFJAk4B7gnIs4ovDQFODQ/PhS4sjD9IzkRvRlYXGuKMjOzobdahXlqTT/7Aj+MiNvzj38ruwIfBu6UdFue9iVScrlE0hHAw8BB+bVr8jZmA8+QEpGZmXVIlQQxU9K1wNbAFyWtB7zQaqGIuIHyfgWAPUrmD+CoCvGYmdkQqJIgjgDGAw9GxDOSNsb/7s3Mel7LBBERL0haAOwgqUpCMTOzHtDyB1/SacDBwN3A8jw5gN+0MS4zM+uwKjWCA4FXRcRz7Q7GzMy6R5XDXB8EVm93IGZm1l2q1CCeAW6TNBV4sRYREce0LSozM+u4KgliSr6ZmdkIUuUopsmS1gC2z5PujYjn2xuWmZl1WpWjmHYnDcs9h3Ti21hJh0aEj2IyM+thVZqYTgf2ioh7ASRtD1wE/F07AzMzs86qchTT6rXkABAR9+GjmszMel6VGsQMSecA5+fnHwRmti8kMzPrBlUSxCdJg+gdQ+qD+A3wv+0MyszMOq/KUUzPAWfkm5l1uXGTru7Iduecul9Htmvt0zBBSLokIt4n6U7KL/25U1sjMzOzjmpWgzg2379rKAIxM7Pu0vAopsLlPv8lIv5UvAH/MjThmZlZp1Q5zPWdJdP2GexAzMysuzTrg/gkqaawjaQ7Ci+tB9zY7sDMzKyzmvVBXAj8DDgFmFSYvjQiHm9rVGZm1nENE0RELAYWA4cASNoEWAtYV9K6EfHw0IRoZmad0LIPQtL+ku4HHgKmkQbt+1mb4zIzsw6r0kn9deDNwH0RsTWwB+6DMDPreVUSxPMR8RiwiqRVIuLXwPg2x2VmZh1WZSymJyWtSxqD6QJJC4Fl7Q3LzMw6rUoN4gDgWeDTwM+BB4D92xmUmZl1XpXB+p4uPJ3cxljMzKyLNDtRbiklg/SRhvyOiFi/bVGZmVnHNTsPYr2hDMTMzLpLyyYmSVuWTfeJcmZmva3KUUzFq4+sBWwN3Au8ti0RmZlZV6jSSb1j8bmknYEj2xaRmZl1hSo1iJeIiFslvbHVfJJ+QLrY0MKIeF2ediLwcWBRnu1LEXFNfu2LwBHAcuCYiPhFX2Oz7tWpy2CaWf9V6YP4TOHpKsDOrPiBb+Zc4FvAeXXT/ysivlG3jR2A95OarcYAv5S0fUQsr7AdMzNrgyonyq1XuK1J6pM4oNVCEfEboOqw4AcAF0fEcxHxEDAb2KXismZm1gZV+iBOGuRtHi3pI8AM4LiIeALYHJhemGdunrYSSROBiQBbbll6gJWZmQ2CKsN9T5B0uaRbJd1Ru/Vze2cB25IG+5sPnF7bTMm8ZSfpERFnR8SEiJgwevTofoZhZmatVOmkvgD4HHAn8MJANhYRC2qPJX0PuCo/nQuMLcy6BTBvINsyM7OBqZIgFkXElMHYmKTNImJ+fvpPwF358RTgQklnkDqptwNuHoxtmplZ/1RJECdI+j4wFXiuNjEiLmu2kKSLgN2BUZLmAicAu0saT2o+mkM+nyIiZkm6BLibNJT4UT6Cycyss6okiMOBVwOrs6KJKYCmCSIiDimZfE6T+U8GTq4Qj5mZDYEqCeL19WdTm5lZ76tyHsT0fCKbmZmNIFVqEG8BDpX0EKkPonY9iJ3aGpmZmXVUlQSxd9ujMDOzrtPsinLrR8QSYOkQxmNmZl2iWQ3iQtJorDNJRy0Vz3YOYJs2xmVmZh3W7JKj78r3Ww9dOGZm1i2qjMW0q6R18uMPSTqj0WVIzcysd1Q5zPUs4BlJrwc+D/wJOL+tUZmZWcdVSRDLIiJI12z4ZkR8k3RtCDMz62FVDnNdmi8H+iHgbZJWJQ27YWZmPaxKDeJg0glyR0TEX0gX8vnPtkZlZmYdV+WKcn8Bzig8f5iVrzNtZmY9pkoNwszMRiAnCDMzK9UwQUiamu9PG7pwzMysWzTrg9hM0m7AuyVdzEuH2iAibm1rZGZm1lHNEsRXgEnAFhQ6qbMA3tGuoMzMrPOajcV0KXCppC9HxNeGMCYzM+sCVQ5z/ZqkdwNvy5Ouj4ir2huWmZl1WpXB+k4BjgXuzrdj8zQzM+thVYba2A8YHxEvAEiaDPwB+GI7AzMzs86qeh7EBoXHL29HIGZm1l2q1CBOAf4g6dekQ13fhmsPZmY9r0on9UWSrgfeSEoQX8jjM5mZWQ+rUoMgIuYDU9oci5mZdRGPxWRmZqWcIMzMrFTTBCFpFUl3DVUwZmbWPZomiHzuw+2SthyieMzMrEtU6aTeDJgl6Wbg6drEiHh326IyM7OOq5IgTmp7FGZm1nVadlJHxDRgDrB6fnwL0PJaEJJ+IGlhsQ9D0kaSrpN0f77fME+XpDMlzZZ0h6Sd+10iMzMbFFUG6/s4cCnw3Txpc+CKCus+F9i7btokYGpEbAdMzc8B9gG2y7eJwFkV1m9mZm1U5TDXo4BdgSUAEXE/sEmrhSLiN8DjdZMPACbnx5OBAwvTz4tkOrCBpM0qxGZmZm1SJUE8FxF/qz2RtBrpinL9sWk+K7t2dnYt0WwOPFKYb26ethJJEyXNkDRj0aJF/QzDzMxaqZIgpkn6ErC2pHcCPwF+OshxqGRaaRKKiLMjYkJETBg9evQgh2FmZjVVEsQkYBFwJ3AkcA1wfD+3t6DWdJTvF+bpc4Gxhfm2AOb1cxtmZjYIqozm+kK+SNBNpH/190ZEf5uYpgCHAqfm+ysL04+WdDHwJmBxrSnKzMw6o2WCkLQf8B3gAVJT0NaSjoyIn7VY7iJgd2CUpLnACaTEcImkI4CHgYPy7NcA+wKzgWeAw/tVGjMzGzRVTpQ7HXh7RMwGkLQtcDXQNEFExCENXtqjZN4gHS1lZmZdokofxMJacsgeZEXfgZmZ9aiGNQhJ78kPZ0m6BriE1AdxEOlsajMz62HNmpj2LzxeAOyWHy8CNmxbRGZm1hUaJoiIcEexmdkIVuUopq2BfwXGFef3cN9mZr2tylFMVwDnkM6efqG94ZiZWbeokiD+GhFntj0SMzPrKlUSxDclnQBcCzxXmxgRLa8JYWZmw1eVBLEj8GHgHaxoYor83MzMelSVBPFPwDbFIb/NzKz3VTmT+nZgg3YHYmZm3aVKDWJT4I+SbuGlfRA+zNXMrIdVSRAntD0KMzPrOlWuBzFtKAIxM7PuUuVM6qWsuPznGsDqwNMRsX47AzMzs86qUoNYr/hc0oHALm2LyMyGpXGTru7Ytuecul/Htt3LqhzF9BIRcQU+B8LMrOdVaWJ6T+HpKsAEVjQ5mZlZj6pyFFPxuhDLgDnAAW2JxszMukaVPghfF8LMbARqdsnRrzRZLiLia22Ix8zMukSzGsTTJdPWAY4ANgacIMzMelizS46eXnssaT3gWOBw4GLg9EbLmZlZb2jaByFpI+AzwAeBycDOEfHEUARmZmad1awP4j+B9wBnAztGxFNDFpWZmXVcsxPljgPGAMcD8yQtybelkpYMTXhmZtYpzfog+nyWtZmZ9Q4nATMzK1XlTGrrEZ0cTM3Mhh/XIMzMrJQThJmZlXKCMDOzUh3pg5A0B1gKLAeWRcSEfFLej4FxpBFj3+eT8szMOqeTNYi3R8T4iJiQn08CpkbEdsDU/NzMzDqkm5qYDiAN50G+P7CDsZiZjXidShABXCtppqSJedqmETEfIN9v0qHYzMyMzp0HsWtEzJO0CXCdpD9WXTAnlIkAW265ZbviMzMb8TpSg4iIefl+IXA5sAuwQNJmAPl+YYNlz46ICRExYfTo0UMVspnZiDPkCULSOvn6EkhaB9gLuAuYAhyaZzsUuHKoYzMzsxU60cS0KXC5pNr2L4yIn0u6BbhE0hHAw8BBHYjNzMyyIU8QEfEg8PqS6Y8Bewx1PGZmVq6bDnM1M7Mu4gRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1KrdToAM7OBGjfp6srz7r3GYwCc2IdlGplz6n4DXkc3c4LogL7szGZmneImJjMzK+UEYWZmpZwgzMyslBOEmZmVcoIwM7NSThBmZlbKCcLMzEr5PAgzs37q5DlNQ3GSnmsQZmZWygnCzMxKdV2CkLS3pHslzZY0qdPxmJmNVF2VICStCnwb2AfYAThE0g6djcrMbGTqtk7qXYDZEfEggKSLgQOAuwd7Qx4wz8ysuW5LEJsDjxSezwXeVJxB0kRgYn76lKR7+7GdUcCj/Yqwe/VamXqtPNB7ZRqW5flu85eHTZl0WuVZy8q0VZUFuy1BqGRavORJxNnA2QPaiDQjIiYMZB3dptfK1Gvlgd4rU6+VB1ymel3VB0GqMYwtPN8CmNehWMzMRrRuSxC3ANtJ2lrSGsD7gSkdjsnMbETqqiamiFgm6WjgF8CqwA8iYlYbNjWgJqou1Wtl6rXyQO+VqdfKAy7TSygiWs9lZmYjTrc1MZmZWZdwgjAzs1IjIkFI2kjSdZLuz/cblsyzlaSZkm6TNEvSJzoRa1UVyzRe0u9zee6QdHAnYq2iSnnyfD+X9KSkq4Y6xqpaDRcjaU1JP86v3yRp3NBHWV2F8rxN0q2Slkl6bydi7KsKZfqMpLvz92aqpErnDXRKhfJ8QtKd+ffthsojVEREz9+A/wAm5ceTgNNK5lkDWDM/XheYA4zpdOwDLNP2wHb58RhgPrBBp2Pvb3nya3sA+wNXdTrmBvGtCjwAbJP3qduBHerm+RfgO/nx+4EfdzruAZZnHLATcB7w3k7HPEhlejvwsvz4kz3wGa1fePxu4OdV1j0iahCk4Tom58eTgQPrZ4iIv0XEc/npmnR/7apKme6LiPvz43nAQmD0kEXYNy3LAxARU4GlQxVUP7w4XExE/A2oDRdTVCzrpcAekspOEu0GLcsTEXMi4g7ghU4E2A9VyvTriHgmP51OOierW1Upz5LC03WoOwG5kW7/ERwsm0bEfIB8v0nZTJLGSrqDNNzHaflHtVtVKlONpF1I/y4eGILY+qNP5eliZcPFbN5onohYBiwGNh6S6PquSnmGm76W6QjgZ22NaGAqlUfSUZIeINXWj6my4q46D2IgJP0SeEXJS/9WdR0R8Qiwk6QxwBWSLo2IBYMVY18NRpnyejYDzgcOjYiO/csbrPJ0uZbDxVScp1sMp1irqlwmSR8CJgC7tTWigalUnoj4NvBtSR8AjgcObbXinkkQEbFno9ckLZC0WUTMzz+WC1usa56kWcBbSU0AHTEYZZK0PnA1cHxETG9TqJUM5mfUxaoMF1ObZ66k1YCXA48PTXh91ovD31Qqk6Q9SX9edis0P3ejvn5GFwNnVVnxSGlimsKKbHkocGX9DJK2kLR2frwhsCvQn5Fih0qVMq0BXA6cFxE/GcLY+qNleYaJKsPFFMv6XuBXkXsPu1AvDn/TskyS3kAa+PXdEdHtf1aqlGe7wtP9gPsrrbnTPfBD1Mu/MTA1vylTgY3y9AnA9/PjdwJ3kI4AuAOY2Om4B6FMHwKeB24r3MZ3Ovb+lic//y2wCHiW9M/pHzsde0lZ9gXuI/X3/Fue9lXSjw3AWsBPgNnAzcA2nY55gOV5Y/4sngYeA2Z1OuZBKNMvgQWF782UTsc8wPJ8E5iVy/Jr4LVV1uuhNszMrNRIaWIyM7M+coIwM7NSThBmZlbKCcLMzEo5QZiZWSknCGs7ScsLo+TenkfKXCW/NkHSmfnxmpJ+mec9WNJb8zK31c5R6TaSvppPqBo2JB1YHM1zOJbBhoYPc7W2k/RURKybH28CXAjcGBEn1M33ZtIYWLvl598BboqIH1bcjkj79HAZNK5UO8uRz9z+Pmk03I6NEmDDg2sQNqQinZU6EThaye6SrsqJ40fA+FxjOBJ4H/AVSRcASPqcpFvyGP0n5WnjJN0j6X+BW4GxkvZSug7GrZJ+IqmWnOZIOilPv1PSq/P0dSX9ME+7Q9I/5+ml6ymSdK7yNRAarb9u/sMkXal0XYt7JZ3QpByH5PXcJem0wjqeknR63s5USaPz9PGSpucyXJ5HBEDS9ZL+XdI04AuuwNo4AAADKUlEQVSk4Z7/M7/P29aVYQ9Jf8jb/YGkNVu8d7vl9dyWl1tvALuHdRknCBtyEfEgad/bpDBtIfAx4LcRMT4ivksaLuBzEfFBSXsB25GGNh4P/J2kt+XFX0UaTuQNpLN5jwf2jIidgRnAZwqbfzRPPwv4bJ72ZWBxROwYETsBv5I0qsV6Gilbf71dgA/mchwkaUJJOZ4HTgPeked7o6TaEOjrALfm7UwDajWx84Av5DLcWZgO6Togu0XEyax4X8dHxIuj+0paCzgXODgidiSN1fbJFmX7LHBURIwnjV32bIX3yIYJJwjrlL5e/2CvfPsD6R/2q0kJA+BPsWIgwjcDOwA3SrqNNOZR8Wpgl+X7maQL3QDsCXy7NkNEPFFhPY2Urb/edRHxWEQ8m+d/S0k53ghcHxGLIg0JfgFQS4gvAD/Oj38EvEXSy0lJYFqePrkwP4X5m3kV8FBE3NdgHWVluxE4Q9IxefvLKmzHhomeGc3Vhg9J2wDLSSO2vqbqYsApuWZRXNc4Uq2hON91EXFIg/XURuVczor9X5QPyd1sPY2Urb9e/bZqz+vLUVWVjsSnW8/ScpsrlS0iTpV0NWksoOmS9oyIP1bYlg0DrkHYkMrt5d8BvhV9O0LiF8BHC/0Jm+d+i3rTgV0lvTLP9zJJ27dY97XA0YUYN+zneqp6p9I1uNcmXTnvxpJ5bgJ2kzRK0qrAIaTmJEjf29q1nz8A3BARi4EnJL01T/9wYf56S4GyvoI/AuNqZW6xDgAkbRsRd0bEaaRmuJX6XWz4cg3ChsLauZlmdWAZ6eJFZ/RlBRFxraTXAL9PB/nwFGm02uV18y2SdBhwUa2DldSXcB+NfZ10IZW78vpOiojL+rGeqm4gvQevBC6MiBm5JlQsx3xJXySNvCngmoioDYH+NPBaSTNJV6M7OE8/FPiOpJcBDwKHN9j+xcD3crNQLdEQEX+VdDjwE6WjnW4hJfNmPiXp7aT37W66+8pr1kc+zNVsCOWkMyEijm41b5N1vHjYsFk7uYnJzMxKuQZhZmalXIMwM7NSThBmZlbKCcLMzEo5QZiZWSknCDMzK/X/xVaAlKTR1lMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(test_stat_dist)\n",
    "plt.axvline(actual, color='gray')\n",
    "    \n",
    "plt.xlabel('Difference in proportions')\n",
    "plt.ylabel('Number of simulations')\n",
    "plt.title('Sampling distribution of the test statistic under H0');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the probability that the test statistic, under the null hypothesis, exceeds the observed difference.\n",
    "\n",
    "This probability is the \"p-value\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.015"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_value = np.mean(test_stat_dist >= actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated p-value is less than the specified threshold, $\\alpha = 0.05$, so the observed difference is considered statistically significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Different test statistics\n",
    "\n",
    "When I present the \"There is only one test\" framework, many of the questions I get are about choosing the test statistic.  In many cases, it seems like an arbitrary choice.\n",
    "\n",
    "And in some ways it is, although rather than \"arbitrary\", I would say it is a modeling decision informed by the context.  In this example, the researchers think there might be an *increase* in the proportion of resistant plants, so I chose a test statistic that measures an increase.\n",
    "\n",
    "If the researchers expected a *change* in resistance without specifying an increase or decrease, I might use the absolute value of the difference as a test statistic, because the absolute value measures magnitude without regard to sign.\n",
    "\n",
    "Using the *difference* in proportions was also a choice; another possibility is to use the *ratio* of proportions.  In the context of epidemiology, it is common to use [risk ratios](https://en.wikipedia.org/wiki/Risk_ratio) to compare proportions.\n",
    "\n",
    "One benefit of the computational framework is that it is easy to try a different test statistic (or a different model of the null hypothesis) and see what effect it has on the results.\n",
    "\n",
    "Here's a function that computes the ratio of proportions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_stat(data):\n",
    "    k1, n1, k2, n2 = data\n",
    "    p1 = k1/n1\n",
    "    p2 = k2/n2\n",
    "    return p2 / p1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the observed risk ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9551282051282053"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "actual = test_stat(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the simulated experiments with this test statistic, and compute the p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.016"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_stat_dist = np.array([test_stat(run_model(data)) \n",
    "                       for i in range(1000)])\n",
    "\n",
    "p_value = np.mean(test_stat_dist >= actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two test statistics yield pretty much the same p-values.  That's good; it means that the interpretation of the results is not sensitive to our modeling choices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}