notebooks/Wilcoxon rank-sum test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wilcoxon rank-sum test\n",
    "```bash\n",
    "by Jeffrey Chijioke-Uche, MSIS, MSIT, CPDS (IBM Sr. Solution Architect, Hybrid Cloud & Multicloud)\n",
    "IBM Information Technology PhD Scholar at Harvard University & Walden University \n",
    "``` \n",
    "\n",
    "In this notebook, I show a numerical example and give a brief explanation of the Wilcoxon rank-sum test.   \n",
    "\n",
    "## Contents\n",
    "1. [Wilcoxon rank-sum test](#WRS)    \n",
    "    1.1 [Ranking the data](#ranking)    \n",
    "    1.2 [Rank sums](#rank-sums)   \n",
    "    1.3 [Test statistic](#test-statistic)     \n",
    "    1.4 [z-score](#z-score)\n",
    "    1.5 [Scipy example](#scipy)   \n",
    "2. [Mann-Whitney test](#MW)   \n",
    "    2.1 [Scipy example](#MW_scipy)   \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='WRS'></a>\n",
    "## Wilcoxon rank-sum test\n",
    "As for other non-parametric tests, the Wilcoxon rank-sum test analyses the ranks of data instead of the bare data.    \n",
    "This is equivalent to the [Mann-Whitney test](#MW).  \n",
    "\n",
    "\n",
    "### Assumptions    \n",
    "This is a non-parametric test, i.e., an assumption-free test.   \n",
    "\n",
    "\n",
    "### Hypothesis\n",
    "|Set|\tNull hypothesis\t| Alternative hypothesis |No. of tails |\n",
    "|---|:-----------------:|:----------------------:|:-----------:|\n",
    "|1\t| $\\mu_1 - \\mu_2 = 0$ | $\\mu_1 - \\mu_2 \\neq 0$ |\t2|\n",
    "\n",
    "This is a two-tailed test.    \n",
    "\n",
    "### Test method\n",
    "\n",
    "* Mean of rank-sum   \n",
    "$\\overline{W}_s = n_1(n_1+n_2+1)/2$\n",
    "\n",
    "\n",
    "* Standard error    \n",
    "SE$_\\overline{W_s}= \\sqrt{n_1 n_2 (n_1 + n_2 + 1)/12}$\n",
    "\n",
    "\n",
    "* $z$-score    \n",
    "$z = \\frac{\\displaystyle X - \\overline{X}}{\\displaystyle s} = \n",
    "\\frac{\\displaystyle W_s - \\overline{W_s}}\n",
    "{\\displaystyle \\textrm{SE}_{\\overline{W_s}}}$\n",
    "\n",
    "\n",
    "\n",
    "### Usage   \n",
    "* This is used to compare the means of two distributions.    \n",
    "* The non-parametric equivalent of the independent t-test, please see my previous notebook [Difference between means](https://github.com/trangel/stats-with-python/blob/master//Difference%20between%20means.ipynb).   \n",
    "\n",
    "### Further reading    \n",
    "* Chapter 15 of Discovering Statistics Using SPSS by Andy Field.   \n",
    "\n",
    "### Numerical example    \n",
    "To start with this notebook, let's use data from table 15.1 of Discovering Statistics Using SPSS by Andy Field.    \n",
    "This data is a drug experiment with 2 groups of 10 participants.\n",
    "A group is giving a drug (Ecstasy) for depresion and the other are only allowed to drink alcohol, the levels of depression are then meassured using the Beck Depression Inventory the day after and midweek. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Drug</th>\n",
       "      <th>BDI Sunday</th>\n",
       "      <th>BDI Wednesday</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Participant</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>15</td>\n",
       "      <td>28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>35</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>16</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>18</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>19</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>17</td>\n",
       "      <td>32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>27</td>\n",
       "      <td>27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>16</td>\n",
       "      <td>29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>13</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>20</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>16</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>15</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>20</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>15</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>16</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>13</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>14</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>19</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>18</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>18</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Drug  BDI Sunday  BDI Wednesday\n",
       "Participant                                    \n",
       "1            Ecstasy          15             28\n",
       "2            Ecstasy          35             35\n",
       "3            Ecstasy          16             35\n",
       "4            Ecstasy          18             24\n",
       "5            Ecstasy          19             39\n",
       "6            Ecstasy          17             32\n",
       "7            Ecstasy          27             27\n",
       "8            Ecstasy          16             29\n",
       "9            Ecstasy          13             36\n",
       "10           Ecstasy          20             35\n",
       "11           Alcohol          16              5\n",
       "12           Alcohol          15              6\n",
       "13           Alcohol          20             30\n",
       "14           Alcohol          15              8\n",
       "15           Alcohol          16              9\n",
       "16           Alcohol          13              7\n",
       "17           Alcohol          14              6\n",
       "18           Alcohol          19             17\n",
       "19           Alcohol          18              3\n",
       "20           Alcohol          18             10"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from IPython.display import Markdown, display\n",
    "def printmd(string):\n",
    "    display(Markdown(string))\n",
    "\n",
    "# Read dataset locally\n",
    "#df = pd.read_csv(\"../data/data-for-drug-experiment.csv\",index_col=0)\n",
    "    \n",
    "# Read dataset from url:\n",
    "import io\n",
    "import requests\n",
    "url=\"https://github.com/schijioke-uche/data-analysis-with-python-an-spss-model/blob/3b9249b72e078eca78a0c38c5aa215a67f21656f/data/data-for-drug-experiment.csv\"\n",
    "s=requests.get(url).content\n",
    "df=pd.read_csv(io.StringIO(s.decode('utf-8')),index_col=0)\n",
    "    \n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ranking'></a>\n",
    "## Ranking the data\n",
    "\n",
    "First we need to compute the ranks for each of the distributions.   \n",
    "Ranking the data is done as follows: finding the lowest score and giving it a rank of 1, the next lowest score is given a rank of 2, etc...   \n",
    "In case of tied-ranks (repeated scores), the rank is averaged.    \n",
    "\n",
    "To rank the data we start by sorting the data.    \n",
    "Then we treat tied-ranks    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rank_distribution(dist,column):\n",
    "    \"\"\"Get distribution ranks for non-parametric test\n",
    "    Column: str, column name to get ranking\n",
    "    Input: dataframe with distribution values.\n",
    "    Output: dataframe including ranks\n",
    "    \"\"\"\n",
    "    # Sort values\n",
    "    dist.sort_values(by=column,inplace=True)\n",
    "    ## Convert series to dataframe, to add ranks\n",
    "    #dist = dist.to_frame()\n",
    "    dist['Potential rank']=0\n",
    "    jj=0\n",
    "    for ii in dist.index:\n",
    "        jj = jj + 1\n",
    "        dist.loc[ii,'Potential rank']=jj\n",
    "        \n",
    "    # Get actual ranks by averaging ranks of repeated scores  \n",
    "    # rank table with average (mean) values of ranks group by BDI scores:   \n",
    "    rank = dist.groupby(column).agg({'Potential rank':'mean'})\n",
    "    rank.columns=['Rank']\n",
    "    #\n",
    "    dist = dist.join(rank,on=column)\n",
    "\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Drug</th>\n",
       "      <th>BDI Sunday</th>\n",
       "      <th>Potential rank</th>\n",
       "      <th>Rank</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Participant</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>13</td>\n",
       "      <td>1</td>\n",
       "      <td>1.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>13</td>\n",
       "      <td>2</td>\n",
       "      <td>1.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>14</td>\n",
       "      <td>3</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Ecstasy</td>\n",
       "      <td>15</td>\n",
       "      <td>4</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>15</td>\n",
       "      <td>6</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Drug  BDI Sunday  Potential rank  Rank\n",
       "Participant                                           \n",
       "16           Alcohol          13               1   1.5\n",
       "9            Ecstasy          13               2   1.5\n",
       "17           Alcohol          14               3   3.0\n",
       "1            Ecstasy          15               4   5.0\n",
       "14           Alcohol          15               5   5.0\n",
       "12           Alcohol          15               6   5.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BDI_Sunday = df[['Drug','BDI Sunday']].copy()\n",
    "BDI_Sunday = rank_distribution(BDI_Sunday,'BDI Sunday')\n",
    "\n",
    "BDI_Sunday.head(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the first two rows have the same BDI and have a rank of 1.5 (the average of 1 and 2) and similar for rows 4 to 6 the rank is 5.0, the average of 4, 5 and 6.\n",
    "\n",
    "We get ranks for the Wednesday's samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Drug</th>\n",
       "      <th>BDI Wednesday</th>\n",
       "      <th>Potential rank</th>\n",
       "      <th>Rank</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Participant</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Alcohol</td>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                Drug  BDI Wednesday  Potential rank  Rank\n",
       "Participant                                              \n",
       "19           Alcohol              3               1   1.0\n",
       "11           Alcohol              5               2   2.0\n",
       "17           Alcohol              6               3   3.5\n",
       "12           Alcohol              6               4   3.5\n",
       "16           Alcohol              7               5   5.0\n",
       "14           Alcohol              8               6   6.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BDI_Wednesday = df[['Drug','BDI Wednesday']].copy()\n",
    "\n",
    "BDI_Wednesday = rank_distribution(BDI_Wednesday,'BDI Wednesday')\n",
    "\n",
    "BDI_Wednesday.head(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rank-sums'></a>\n",
    "## Rank - sums  \n",
    "Let's compute the rank sums $W_s$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Rank sum</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Drug</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alcohol</th>\n",
       "      <td>90.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ecstasy</th>\n",
       "      <td>119.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Rank sum\n",
       "Drug             \n",
       "Alcohol      90.5\n",
       "Ecstasy     119.5"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For Sunday's BDI:\n",
    "\n",
    "rank_sum = BDI_Sunday[['Drug','Rank']].groupby('Drug').agg('sum')\n",
    "rank_sum.columns=['Rank sum']\n",
    "rank_sum_sunday = rank_sum\n",
    "rank_sum\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id ='test-statistic'></a>\n",
    "### Test-statistic\n",
    "When the groups have unequal numbers of participants in them then $W_s$ is simply the sum of ranks in the group that contains the fewer people. When the group sizes are equal it's the value of the smaller summed rank.    \n",
    "In this case, the group sizes are equal, and hence the lowest of the rank sum (90.5) is our test statistic $W_s$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$W_s$ 90.5"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "W = rank_sum.values.min()\n",
    "printmd('$W_s$ {}'.format(W))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='z-score'></a>\n",
    "### $z$-score    \n",
    "Remember    \n",
    "$z = \\frac{\\displaystyle X - \\overline{X}}{\\displaystyle s} = \n",
    "\\frac{\\displaystyle W_s - \\overline{W_s}}\n",
    "{\\displaystyle \\textrm{SE}_{\\overline{W_s}}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "** Sunday **"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$W_s$ 90.5"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$W_s$ mean 105.0"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "SE 13.23"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$z$-score -1.1"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "n1 = 10\n",
    "n2 = 10\n",
    "W_mean = n1 *(n1 + n2 + 1.0)/2.0\n",
    "SE  = (n1*n2 * (n1+n2+1.0)/12.0)**0.5\n",
    "\n",
    "# Using the z-score equation the Sunday score is:   \n",
    "z = (W - W_mean)/SE\n",
    "\n",
    "printmd('** Sunday **')\n",
    "printmd('$W_s$ {}'.format(W))\n",
    "printmd('$W_s$ mean {}'.format(W_mean))\n",
    "printmd('SE {}'.format(round(SE,2)))\n",
    "printmd('$z$-score {}'.format(round(z,2)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to look for the probability for this $z$-score on a table, or using the CDF for a t-distribution as below.   \n",
    "Remember that we get the z-distribution using the t-distribution with large degrees of freedom.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "p-value 0.2731"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "DF=10000  # set a large degrees of freedom.\n",
    "pvalue = 2* stats.t.cdf(z,DF) # get the P value from the CDF of a t-distribution\n",
    "# Multiply by two, since this is two-tailed test\n",
    "\n",
    "printmd('p-value {}'.format(round(pvalue,4)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "** Wednesday **"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$W_s$ 59.0"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$W_s$ mean 105.0"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "SE 13.23"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$z$-score -3.48"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "p-value 0.0005"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute now the Wednesday z-score\n",
    "# Using the z-score equation the Sunday score is:   \n",
    "\n",
    "rank_sum = BDI_Wednesday[['Drug','Rank']].groupby('Drug').agg('sum')\n",
    "rank_sum.columns=['Rank sum']\n",
    "\n",
    "W = rank_sum.values.min()\n",
    "\n",
    "z = (W - W_mean)/SE\n",
    "\n",
    "# Get the P-value as above\n",
    "DF=10000  # set a large degrees of freedom.\n",
    "pvalue = 2* stats.t.cdf(z,DF) # get the P value from the CDF of a t-distribution\n",
    "# Multiply by two, since this is two-tailed test\n",
    "\n",
    "\n",
    "printmd('** Wednesday **')\n",
    "printmd('$W_s$ {}'.format(W))\n",
    "printmd('$W_s$ mean {}'.format(W_mean))\n",
    "printmd('SE {}'.format(round(SE,2)))\n",
    "printmd('$z$-score {}'.format(round(z,2)))\n",
    "printmd('p-value {}'.format(round(pvalue,4)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='scipy'></a>\n",
    "## Using scipy   \n",
    "\n",
    "Let's run the same test using python's scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "** Wednesday **"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$z$-score -3.48"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "p-value 0.0005"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "** Sunday **"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$z$-score -1.1"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "p-value 0.273"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "# Let's do the Wednesday sample as an example:   \n",
    "# SELECT 'BDI Wednesday' FROM df WHERE 'Drug'='Alcohol'\n",
    "x = df.loc[df['Drug']=='Alcohol']['BDI Wednesday'].values\n",
    "# SELECT 'BDI Wednesday' FROM df WHERE 'Drug'='Ecstasy'\n",
    "y = df.loc[df['Drug']=='Ecstasy']['BDI Wednesday'].values\n",
    "\n",
    "results = stats.ranksums(x, y)\n",
    "z = results[0]\n",
    "pvalue = results[1]\n",
    "printmd('** Wednesday **')\n",
    "printmd('$z$-score {}'.format(round(z,2)))\n",
    "printmd('p-value {}'.format(round(pvalue,4)))\n",
    "\n",
    "# Let's do the Sunday sample :   \n",
    "# SELECT 'BDI Sunday' FROM df WHERE 'Drug'='Alcohol'\n",
    "x = df.loc[df['Drug']=='Alcohol']['BDI Sunday'].values\n",
    "# SELECT 'BDI Sunday' FROM df WHERE 'Drug'='Ecstasy'\n",
    "y = df.loc[df['Drug']=='Ecstasy']['BDI Sunday'].values\n",
    "\n",
    "results = stats.ranksums(x, y)\n",
    "z = results[0]\n",
    "pvalue = results[1]\n",
    "printmd('** Sunday **')\n",
    "printmd('$z$-score {}'.format(round(z,2)))\n",
    "printmd('p-value {}'.format(round(pvalue,4)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='MW'></a>\n",
    "## Mann-Whitney test     \n",
    "\n",
    "As mentioned before, the Mann-Whitney test is equivalent to the Wilcoxon rank-sum test (above). This is based on a test statistic $U$,\n",
    "\n",
    "$ U = n_1 n_2 + \\frac{\\displaystyle n_1(n_1+1)}{\\displaystyle 2} - R_1$,   \n",
    "\n",
    "where $n_i$ is the sample size for the $i^\\textrm{th}$ group and $R_1$ is the sum of ranks for group 1.   \n",
    "\n",
    "In this case, the rank sums for Sunday are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Rank sum</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Drug</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Alcohol</th>\n",
       "      <td>90.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ecstasy</th>\n",
       "      <td>119.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Rank sum\n",
       "Drug             \n",
       "Alcohol      90.5\n",
       "Ecstasy     119.5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rank_sum_sunday"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "**Sunday**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "Rank sum for group 1 (Ecstasy data) is 119.5"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "Test statistic, U = 35.5"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R1=rank_sum_sunday['Rank sum'][1]\n",
    "printmd('**Sunday**')\n",
    "printmd('Rank sum for group 1 (Ecstasy data) is {}'.format(R1))\n",
    "\n",
    "n1=10\n",
    "n2=10\n",
    "U = n1*n2 + n1*(n1+1)/2 - R1\n",
    "printmd('Test statistic, U = {}'.format(U))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a id='MW_scipy'></a>\n",
    "## Mann-Whitney test with scipy   \n",
    "Let's perform the test with Python's scipy   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MannwhitneyuResult(statistic=35.5, pvalue=0.14304106361646457)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's do the Sunday sample :   \n",
    "# SELECT 'BDI Sunday' FROM df WHERE 'Drug'='Alcohol'\n",
    "x = df.loc[df['Drug']=='Alcohol']['BDI Sunday'].values\n",
    "# SELECT 'BDI Sunday' FROM df WHERE 'Drug'='Ecstasy'\n",
    "y = df.loc[df['Drug']=='Ecstasy']['BDI Sunday'].values\n",
    "\n",
    "stats.mannwhitneyu(x, y, use_continuity=True, alternative=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pvalue for Sunday is 0.14 (higher than a significance of 0.05) in agreement with the Wilcoxon rank-sum test, above."
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}