solutions/ch10q07.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv('../data/snodgrass.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>author</th>\n",
       "      <th>total</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.225</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.229</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>twain</td>\n",
       "      <td>0.217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.196</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.210</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.202</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>snodgrass</td>\n",
       "      <td>0.201</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       author  total\n",
       "0       twain  0.225\n",
       "1       twain  0.262\n",
       "2       twain  0.217\n",
       "3       twain  0.240\n",
       "4       twain  0.230\n",
       "5       twain  0.229\n",
       "6       twain  0.235\n",
       "7       twain  0.217\n",
       "8   snodgrass  0.209\n",
       "9   snodgrass  0.205\n",
       "10  snodgrass  0.196\n",
       "11  snodgrass  0.210\n",
       "12  snodgrass  0.202\n",
       "13  snodgrass  0.207\n",
       "14  snodgrass  0.224\n",
       "15  snodgrass  0.223\n",
       "16  snodgrass  0.220\n",
       "17  snodgrass  0.201"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"2\" halign=\"left\">total</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>author</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>snodgrass</th>\n",
       "      <td>0.209700</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>twain</th>\n",
       "      <td>0.231875</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              total      \n",
       "               mean count\n",
       "author                   \n",
       "snodgrass  0.209700    10\n",
       "twain      0.231875     8"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby('author').aggregate(['mean', 'count'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$t_{obs}$$ is \n",
    "\n",
    "$$|.2097 - .231875| = .022175$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_obs = .022175"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutation test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full set of observations in $$X_1,\\ldots,X_{10}, Y_1,\\ldots,Y_{8}$$\n",
    "\n",
    "For all 18! permutations of the data, compute\n",
    "\n",
    "$$|(\\bar{X}) - (\\bar{Y})|$$\n",
    "\n",
    "where $$(\\bar{X})$$ is the mean of the first 10 observations, and $$(\\bar{Y})$$ is the mean of the last 8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from itertools import permutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = list(data['total'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.5\n",
      "0.0\n",
      "-1.5\n",
      "1.5\n",
      "0.0\n",
      "1.5\n"
     ]
    }
   ],
   "source": [
    "perms = (p for p in permutations([1,2,3]))\n",
    "for p in perms:\n",
    "    print(np.mean(p[:2]) - np.mean(p[2:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# there are too many permutations to count all of them, so we set some limit N\n",
    "# on the total number of observed permutations\n",
    "\n",
    "\n",
    "def permutation_test(obs, m, c, N):\n",
    "    i = 0 # counter for T>c\n",
    "    \n",
    "    # adjust m so there's at least one observation in each sample\n",
    "    m = min(m, len(obs)-1)\n",
    "\n",
    "    # evaluate permutations\n",
    "    for j in range(N):\n",
    "        p = np.random.permutation(obs)\n",
    "        # calculate stat\n",
    "        T = np.mean(p[:m]) - np.mean(p[m:])\n",
    "        # count if > c\n",
    "        if T > c:\n",
    "            i += 1\n",
    "            \n",
    "    # i/N approximates the p-value as N -> N!\n",
    "    return i/N\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00023"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "permutation_test(obs, 10, t_obs, N=100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wald test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Wald test of two means can be constructed as (see example 10.7)\n",
    "\n",
    "$$ W = \\frac{\\bar{X}-\\bar{Y}}{\\sqrt{\\frac{s_x^2}{m}+\\frac{s_y^2}{n}}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = data[data['author']=='snodgrass']['total']\n",
    "Y = data[data['author']=='twain']['total']\n",
    "\n",
    "x_bar = np.mean(X)\n",
    "m = len(X)\n",
    "s_x = np.std(X)\n",
    "\n",
    "y_bar = np.mean(Y)\n",
    "n = len(Y)\n",
    "s_y = np.std(Y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = (x_bar - y_bar)/np.sqrt(s_x**2/m + s_y**2/n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we just need to test that $$W > z_{\\alpha/2}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(W) > 1.96"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the p-value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.9963324780729e-05"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm.cdf(W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the CI we need $$\\hat{\\theta}$$, which here is $$\\bar{X}-\\bar{Y}$$, and se_hat, which is $$\\sqrt{\\frac{s_x^2}{m}+\\frac{s_y^2}{n}}$$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "theta_hat = x_bar-y_bar\n",
    "se_hat = np.sqrt(s_x**2/m + s_y**2/n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.022175, -0.010931838189370394, -0.033418161810629607)"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(theta_hat, theta_hat+2*se_hat, theta_hat-2*se_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "Our p-value from both tests is very small, which is strong evidence for rejecting the null that the means of the two samples are the same.\n",
    "\n",
    "Additionally, the confidence intervals do not overlap, which again is evidence against $$H_0$$"
   ]
  }
 ],
 "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": 2
}