Chapter16_Confidence_intervals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting up a fancy stylesheet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Arvo' sans-serif;\n",
       "    line-height: 130%;\n",
       "    font-size: 115%;\n",
       "    width:1000px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "css_file = 'style.css'\n",
    "HTML(open(css_file, 'r').read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting up the required python &#8482; environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import bayes_mvs\n",
    "#from math import factorial\n",
    "import scikits.bootstrap as bs\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from warnings import filterwarnings\n",
    "sns.set_style('white')\n",
    "sns.set_context('paper', font_scale = 2.0, rc = {'lines.linewidth': 1.5, 'figure.figsize' : (10, 8)})\n",
    "%matplotlib inline\n",
    "filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importing the mock dataset and checking if it imported correctly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>File</th>\n",
       "      <th>Age</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Delay</th>\n",
       "      <th>Stay</th>\n",
       "      <th>ICU</th>\n",
       "      <th>RVD</th>\n",
       "      <th>CD4</th>\n",
       "      <th>HR</th>\n",
       "      <th>Temp</th>\n",
       "      <th>CRP</th>\n",
       "      <th>WCC</th>\n",
       "      <th>HB</th>\n",
       "      <th>Rupture</th>\n",
       "      <th>Histo</th>\n",
       "      <th>Comp</th>\n",
       "      <th>MASS</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>38</td>\n",
       "      <td>Female</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>97</td>\n",
       "      <td>35.2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10.49</td>\n",
       "      <td>10.4</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>32</td>\n",
       "      <td>Male</td>\n",
       "      <td>6</td>\n",
       "      <td>10</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>57</td>\n",
       "      <td>109</td>\n",
       "      <td>38.8</td>\n",
       "      <td>45.3</td>\n",
       "      <td>7.08</td>\n",
       "      <td>19.8</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>19</td>\n",
       "      <td>Female</td>\n",
       "      <td>1</td>\n",
       "      <td>16</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>120</td>\n",
       "      <td>36.3</td>\n",
       "      <td>10.7</td>\n",
       "      <td>13.00</td>\n",
       "      <td>8.7</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   File  Age  Gender  Delay  Stay ICU  RVD  CD4   HR  Temp   CRP    WCC    HB  \\\n",
       "0     1   38  Female      3     6  No   No  NaN   97  35.2   NaN  10.49  10.4   \n",
       "1     2   32    Male      6    10  No  Yes   57  109  38.8  45.3   7.08  19.8   \n",
       "2     3   19  Female      1    16  No   No  NaN  120  36.3  10.7  13.00   8.7   \n",
       "\n",
       "  Rupture Histo Comp  MASS  \n",
       "0      No   Yes  Yes     5  \n",
       "1      No    No  Yes     8  \n",
       "2      No    No   No     3  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('MOOC_Mock.csv')\n",
    "data.head(3) #Looking at the first 3 rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Confidence intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>When conducting research by scrutinizing a sample set of individuals taken from a larger population, we create statistics.  We then use these statistics to infer parameter values for the greater population.  It can be used to guide patient treatment and management, by inferring that our results are applicable to that larger (patient) population.</p>\n",
    "<p>However, there is no way that we can say with absolute confidence that our statistic is exactly what would to be found in the general population.  Instead we can place an arbitrary confidence around our test statistics, and infer that we are sure that the real (larger population) parameters are within that arbitrary confidence that we placed in our statistic(s).</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence interval around mean of a sample variable set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Most textbooks begin with the chapter on confidence intervals (CI's) by using an example where the population mean and standard deviation is known.  This is almost never the case in healthcare research, so let's have a look at what is commonly observed.</p>\n",
    "<p>Most of the time we are stuck with a set of variable values from our sample set and we can calculate the mean, standard deviation and standard error of that set.</p>\n",
    "<p>In our example dataset we might look at the white cell count (WCC) on admission.  We could divide our sample patients into two groups, say those with and without retroviral disease.</p>\n",
    "<p>In the code below we make two new DataFrames, with each DataFrame object attached to a new computer variable.<p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RVD_neg = data[data['RVD'] == 'No'] # Only including the rows in the original DataFrame that have a Yes in the RVD column\n",
    "RVD_neg_appx = RVD_neg[RVD_neg['Histo'] == 'Yes'] # Extracting only those who actually had appendicitis\n",
    "RVD_pos = data[data['RVD'] == 'Yes'] # Doing the same but with No in the RVD column\n",
    "RVD_pos_appx = RVD_pos[RVD_pos['Histo'] == 'Yes'] # Extracting only those who actually had appendicitis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Let's look at the first rows in each new DataFrame using the *.head()* function.  Leaving the argument blank will result in a default of *5* rows.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>File</th>\n",
       "      <th>Age</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Delay</th>\n",
       "      <th>Stay</th>\n",
       "      <th>ICU</th>\n",
       "      <th>RVD</th>\n",
       "      <th>CD4</th>\n",
       "      <th>HR</th>\n",
       "      <th>Temp</th>\n",
       "      <th>CRP</th>\n",
       "      <th>WCC</th>\n",
       "      <th>HB</th>\n",
       "      <th>Rupture</th>\n",
       "      <th>Histo</th>\n",
       "      <th>Comp</th>\n",
       "      <th>MASS</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>38</td>\n",
       "      <td>Female</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>97</td>\n",
       "      <td>35.2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>10.49</td>\n",
       "      <td>10.4</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>8</td>\n",
       "      <td>19</td>\n",
       "      <td>Male</td>\n",
       "      <td>2</td>\n",
       "      <td>13</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>102</td>\n",
       "      <td>38.5</td>\n",
       "      <td>224.6</td>\n",
       "      <td>8.99</td>\n",
       "      <td>15.5</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>9</td>\n",
       "      <td>46</td>\n",
       "      <td>Male</td>\n",
       "      <td>2</td>\n",
       "      <td>9</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>107</td>\n",
       "      <td>36.2</td>\n",
       "      <td>385.1</td>\n",
       "      <td>10.26</td>\n",
       "      <td>15.5</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10</td>\n",
       "      <td>19</td>\n",
       "      <td>Male</td>\n",
       "      <td>7</td>\n",
       "      <td>12</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>105</td>\n",
       "      <td>36.6</td>\n",
       "      <td>123</td>\n",
       "      <td>22.19</td>\n",
       "      <td>13.7</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>11</td>\n",
       "      <td>33</td>\n",
       "      <td>Male</td>\n",
       "      <td>2</td>\n",
       "      <td>14</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>NaN</td>\n",
       "      <td>104</td>\n",
       "      <td>36.8</td>\n",
       "      <td>NaN</td>\n",
       "      <td>16.20</td>\n",
       "      <td>12.7</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    File  Age  Gender  Delay  Stay  ICU RVD  CD4   HR  Temp    CRP    WCC  \\\n",
       "0      1   38  Female      3     6   No  No  NaN   97  35.2    NaN  10.49   \n",
       "7      8   19    Male      2    13   No  No  NaN  102  38.5  224.6   8.99   \n",
       "8      9   46    Male      2     9   No  No  NaN  107  36.2  385.1  10.26   \n",
       "9     10   19    Male      7    12  Yes  No  NaN  105  36.6    123  22.19   \n",
       "10    11   33    Male      2    14   No  No  NaN  104  36.8    NaN  16.20   \n",
       "\n",
       "      HB Rupture Histo Comp  MASS  \n",
       "0   10.4      No   Yes  Yes     5  \n",
       "7   15.5      No   Yes  Yes     5  \n",
       "8   15.5      No   Yes   No     1  \n",
       "9   13.7      No   Yes  Yes     5  \n",
       "10  12.7      No   Yes   No     9  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVD_neg_appx.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>File</th>\n",
       "      <th>Age</th>\n",
       "      <th>Gender</th>\n",
       "      <th>Delay</th>\n",
       "      <th>Stay</th>\n",
       "      <th>ICU</th>\n",
       "      <th>RVD</th>\n",
       "      <th>CD4</th>\n",
       "      <th>HR</th>\n",
       "      <th>Temp</th>\n",
       "      <th>CRP</th>\n",
       "      <th>WCC</th>\n",
       "      <th>HB</th>\n",
       "      <th>Rupture</th>\n",
       "      <th>Histo</th>\n",
       "      <th>Comp</th>\n",
       "      <th>MASS</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>28</td>\n",
       "      <td>Female</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>491</td>\n",
       "      <td>115</td>\n",
       "      <td>37.1</td>\n",
       "      <td>51.6</td>\n",
       "      <td>21.98</td>\n",
       "      <td>13.4</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7</td>\n",
       "      <td>36</td>\n",
       "      <td>Male</td>\n",
       "      <td>6</td>\n",
       "      <td>14</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>203</td>\n",
       "      <td>95</td>\n",
       "      <td>37.5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>16.21</td>\n",
       "      <td>17.7</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>20</td>\n",
       "      <td>33</td>\n",
       "      <td>Female</td>\n",
       "      <td>1</td>\n",
       "      <td>19</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>203</td>\n",
       "      <td>123</td>\n",
       "      <td>37.5</td>\n",
       "      <td>46.5</td>\n",
       "      <td>8.00</td>\n",
       "      <td>10.4</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>32</td>\n",
       "      <td>39</td>\n",
       "      <td>Male</td>\n",
       "      <td>2</td>\n",
       "      <td>10</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>102</td>\n",
       "      <td>100</td>\n",
       "      <td>36.8</td>\n",
       "      <td>208</td>\n",
       "      <td>20.43</td>\n",
       "      <td>12.9</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>33</td>\n",
       "      <td>33</td>\n",
       "      <td>Male</td>\n",
       "      <td>0</td>\n",
       "      <td>21</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>266</td>\n",
       "      <td>100</td>\n",
       "      <td>39.7</td>\n",
       "      <td>194.1</td>\n",
       "      <td>7.31</td>\n",
       "      <td>10.7</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    File  Age  Gender  Delay  Stay ICU  RVD  CD4   HR  Temp    CRP    WCC  \\\n",
       "4      5   28  Female      3     3  No  Yes  491  115  37.1   51.6  21.98   \n",
       "6      7   36    Male      6    14  No  Yes  203   95  37.5    NaN  16.21   \n",
       "19    20   33  Female      1    19  No  Yes  203  123  37.5   46.5   8.00   \n",
       "31    32   39    Male      2    10  No  Yes  102  100  36.8    208  20.43   \n",
       "32    33   33    Male      0    21  No  Yes  266  100  39.7  194.1   7.31   \n",
       "\n",
       "      HB Rupture Histo Comp  MASS  \n",
       "4   13.4      No   Yes   No     7  \n",
       "6   17.7     Yes   Yes  Yes     7  \n",
       "19  10.4      No   Yes  Yes     8  \n",
       "31  12.9     Yes   Yes   No     1  \n",
       "32  10.7     Yes   Yes  Yes     3  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVD_pos_appx.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the two **RVD** columns contain only *No* and *Yes* respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>The easiest way to calculate the confidence intervals is using the *scipy.stats* function *bayes_mvs*.  It actually gives nine results, three rows of three values each.  The first row gives the mean and then the lower as well as the upper bound of the stated confidence interval.</p>\n",
    "<p>The second row does the same but for the variance. The last row calculates the values for the standard deviation.</p>\n",
    "<p>Let's first describe our two DataFrames then calculate the CI's.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    78.000000\n",
       "mean     14.044103\n",
       "std       4.347939\n",
       "min       6.610000\n",
       "25%      11.270000\n",
       "50%      13.655000\n",
       "75%      16.305000\n",
       "max      26.400000\n",
       "Name: WCC, dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVD_neg_appx['WCC'].describe() # Describing the WCC column in the RVD_neg DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    40.000000\n",
       "mean     15.773250\n",
       "std       4.934148\n",
       "min       2.950000\n",
       "25%      12.475000\n",
       "50%      15.575000\n",
       "75%      19.920000\n",
       "max      24.890000\n",
       "Name: WCC, dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RVD_pos_appx['WCC'].describe() # Describing the WCC column in the RVD_pos DataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Now for the CI's.  The *bayes_mvs()* function takes two arguments, the first is the values themselves and the second is the required confidence interval. In this example I want the 95% CI's.  I also make sure to discard all non-numerical values.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Mean(statistic=14.044102564102564, minmax=(13.063793818374444, 15.024411309830684)),\n",
       " Variance(statistic=19.408694495726504, minmax=(14.110883384251707, 26.648893927036038)),\n",
       " Std_dev(statistic=4.3908697719337093, minmax=(3.7564455785025959, 5.1622566700074142)))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bayes_mvs(RVD_neg_appx['WCC'].dropna(), 0.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Mean(statistic=15.773250000000001, minmax=(14.195232891628422, 17.35126710837158)),\n",
       " Variance(statistic=25.661807500000002, minmax=(16.336646621396056, 40.140096800828552)),\n",
       " Std_dev(statistic=5.0316399443801929, minmax=(4.041861776631662, 6.33562126399839)))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bayes_mvs(RVD_pos_appx['WCC'].dropna(), 0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Let's put down in words what we have found from the calculations above.  We could state in a paper that there were 78 patients in the RVD negative group with acute appendicitis.  They had a mean WCC of 14.04 on admission, with a 95% confidence interval of 13.06 to 15.02.  There were 40 RVD positive patients with an average admission WCC of 15.77 (95% CI 14.20 to 17.35). (Note the mathematical rounding).</p>\n",
    "<p>However, what do we mean by this 95% CI? Well, by our statistical anlaysis, we found a mean WCC of 15.77 in RVD positive patients.  If want want to infer this on the larger population we will have to say that according to our analysis we think that the true population admission WCC value out there for all RVD positive patients is between two values.  We can state those two values (lower and upper) with a set level of confidence.  Here, we are 95% confident that in the larger population out there, we think the true value lies between 14.20 and 17.35.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Let's drop our CI to 80%.  What do you expect to happen?  Well, the lower and upper values will be closer to the mean, making us less confident that our *prediction* (inference) is correct.  Let's see:</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Mean(statistic=15.773250000000001, minmax=(14.756206821106424, 16.790293178893577)),\n",
       " Variance(statistic=25.661807500000002, minmax=(18.742423588894685, 33.674780516226335)),\n",
       " Std_dev(statistic=5.0316399443801929, minmax=(4.329252081929936, 5.8029975457711798)))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bayes_mvs(RVD_pos_appx['WCC'], 0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, the values are now in a much narrower range of 14.76 to 16.79.  We are much less confident about this, though!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What lies behind these calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Let's start with the equation:</p>\n",
    "$$ \\overline { x } \\pm t{ s }_{ SE } $$\n",
    "<p>Now *t* is the t-statistic.  The t-statistic is calculated (or read from a table) and depends on the value of confidence desired (actually alpha value) and the degrees of freedom.</p>\n",
    "<p>The *s<sub>SE</sub>* is the standard error of the mean or...</p>\n",
    "$$ \\frac { s }{ \\sqrt { n }  }  $$\n",
    "<p>...the standard deviation over the square root of the sample size.</p>\n",
    "<p>Given the degrees of freedom and alpha value, a t-distribution is calculated.  It looks like a normal distribution and is different for any given parameters.  Therefore, when we state a confidence interval, we calculate what middle area will be covered by that confidence, i.e. we'll cover 0.95 of the total area under the curve and calculate where on the x-axis these cut-offs fall (the t-statistic).  This we multiply by the standard error of the mean, and subtract as well as add to the mean to get the lower and upper bounds of the required CI.  Neat!</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrapping a confidence interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>There is another method to calculate the CI.  It involves bootstrapping or resampling.  We are interested in infereing population parameters from statistics.  We would like to say that the real population parameter is between such and such a value and we are confident about this *guess* at a given level.</p>\n",
    "<p>The only knowledge we have about the population, though, is our sample set.  Actually we have nothing else.  What we good do is take random samples from our sample set.  We need to take exactely the same number of random samples from our sample set.  In order not to have identical resamples, we allow a given sample value taken to be *thrown back* into the pool and have a chance of being resampled again, and again.  We can do thos thousanfs of times.  With the dawn of the age of fast computers, this is no problem.  Now we have as big a sample set as we want.  The hope is that this is a good refelction of the population.  From it we can form a distribution curve and calculate the CI<sup>s</sup>.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 14.23725,  17.23825])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bs.ci(RVD_pos_appx['WCC'].dropna(), alpha = 0.05, n_samples = 10000)\n",
    "# Note that we state the alpha value directly and here I use 10000 resamples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>I can now given the CI<sup>s</sup> as the values in the array above.  Note, it will be different everytime you do it!</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}