permutation_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Permutation test session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib as plt\n",
    "import itertools\n",
    "from scipy import stats\n",
    "from scipy.special import comb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_t_statistic(data, ind, var):\n",
    "    d = data.loc[data.index.isin(ind), var].values\n",
    "    p = data.loc[~data.index.isin(ind), var].values\n",
    "    return stats.ttest_ind(d,p)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean_difference(data, ind, var):\n",
    "    d = data.loc[data.index.isin(ind), var].values\n",
    "    p = data.loc[~data.index.isin(ind), var].values\n",
    "    return np.mean(d) - np.mean(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Outline\n",
    "\n",
    "- [Experimental situation](#1)\n",
    "- [Intuitive idea](#2)\n",
    "- [Permutation test in action: a t-test example](#3)\n",
    "- [Permutation test in action: another statistic](#4)\n",
    "- [Montecarlo or approximated permutation test](#5)\n",
    "- [Permutation test and AB testing](#6)\n",
    "- [(Brief) Introduction to Bootstrap](#7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='1'></a>\n",
    "## Experimental situation\n",
    "\n",
    "We are part of a team in charge of analysing a randomised study which aim is to check if there exists a relationship between a drug (*D*) and the level of ethinyl estradiol (*EE*), which is a estrogen medication.\n",
    "\n",
    "To carry out the study, a (random) sample of women is obtained. All of them are taken contraceptive with ethinyl estradiol (*EE*). We want to study the influence of the drug *D* on the levels of *EE*. Such a level is measured by the variable *area under the curve* (*AUC*). That is: *EEAUC*.\n",
    "\n",
    "In total, there are 16 women. 8 receive the drug *D* and 8 receives a placebus *P*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('./data/women.csv')\n",
    "print(data.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data.treat.value_counts())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.groupby('treat')['EEAUC'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.boxplot(by='treat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could proceed with *t-test* to check if there exists a difference between both populations (*p-value*, bla bla bla). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = data.loc[data.treat == 'D']['EEAUC'].values\n",
    "p = data.loc[data.treat == 'P']['EEAUC'].values\n",
    "t_estimator, p_value = stats.ttest_ind(d,p)\n",
    "print(f'The value of the estimator is {t_estimator}. The p-value is {p_value}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative way of analysing this data is by means of **permutation test**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='2'></a>\n",
    "# Intuitive idea\n",
    "\n",
    "If the two samples that we are comparing (*D* and *P*) have the variable *AUC* equal, belonging to *D* or *P* is simply a label and it has no relation with the fact that AUC is large or small systematically. With this statement we are **establishing $H_0$**.\n",
    "\n",
    "As a consequence, it is equally probable to observe the previous sample and any permutation of the sample in which the values of *AUC* belong to *D* or *P* in another way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data = data.copy()\n",
    "new_data['treat'] = np.random.permutation(new_data.treat)\n",
    "print('Original data:\\n', data.head(), '\\n\\n Permuation:\\n', new_data.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A method of pure hard strength consists in:\n",
    "\n",
    "- Calculating the statistic on the original data\n",
    "- Numbering all possible permutatuins of the data\n",
    "- Computing a proper statistic for each permutation \n",
    "- Counting how many times this value is more extrem than the original one\n",
    "\n",
    "In our case, there are *16!* permutations. This is not feaseable at all!!!!. However, since some of the *16!* perms are repeated, we could avoid some of them (Why??). In the end, we need:\n",
    "\n",
    "$$\\frac{16!}{8!\\cdot8!} = 12870$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='3'></a>\n",
    "# Permutation test in action: a t-test example\n",
    "\n",
    "The way to proceed with the permutation test is the following one:\n",
    "\n",
    "1. Calculation of statistic for the original data\n",
    "2. List of all possible combinations\n",
    "3. Calculation of statistic for each one of these indexes resortings\n",
    "4. Finally, *p-value*: proportion of times that the statistic computed from permutations is equal or more extrem than the one compuded from the original data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. t-statistic for the original data\n",
    "d = data.loc[data.treat == 'D']['EEAUC'].values\n",
    "p = data.loc[data.treat == 'P']['EEAUC'].values\n",
    "t_orig, _ = stats.ttest_ind(d,p)\n",
    "print(f'The t-statistic for the original data is {t_orig}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. List of all possible combinations\n",
    "all_combinations_d = itertools.combinations(range(len(data.index)), sum(data.treat == 'D'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does itertools.combinations do?\n",
    "all_comb = itertools.combinations(range(5), 3)\n",
    "for comb in all_comb:\n",
    "    print(comb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. Calculation of statistic for each one of these indexes resortings\n",
    "result = [get_t_statistic(data=data, ind=sample, var='EEAUC') for sample in all_combinations_d]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Total number of permutations: {len(result)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pyplot.hist(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. Finally, p-value\n",
    "sum(result>=t_orig)/len(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is this *p-value* different from the *p-value* of the *t-test*???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details><summary>CLICK ME</summary>\n",
    "<p>\n",
    "The p-value we have just obtained refers to unilateral test in which we consider the alternative hypothesis:\n",
    "    \n",
    "$$H_1 : D > P$$\n",
    "\n",
    "In order to perform the bilateral test with $H_1 : D \\ne P$, the most appropriate statistic is the absolute value of t-Student because if both negative and positive differences are far away from zero, there are evidences against the null hypothesis.\n",
    "    \n",
    "</p>\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Place the code here\n",
    "# # 4.1 Finally, p-value\n",
    "sum(np.abs(result)>=t_orig)/len(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='4'></a>\n",
    "# Permutation test in action: another statistic\n",
    "\n",
    "we have just seen a use case of permutation test. However, we could avoid it because we could use *t-test* to analyse the hypothesis test. \n",
    "\n",
    "**What are the benefits of permutation tests?**\n",
    "\n",
    "One of the advantages of permutation tests is that we can use the statistic that best expresses the type of differences we want to show, and we are not forced to use a statistic which distribution under null hypothesis is mathematically easy to determine and with good statistical prorperties. \n",
    "\n",
    "The permutation mechanism allows to ingnore the necessity of good mathematical behavior.\n",
    "Computation of p-value is easy\n",
    "\n",
    "In this part we are going to analyse the same data but, instead of using *t-statistic*, we are going to use the following one:\n",
    "\n",
    "$$diff = mean(d) - mean(p)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_combinations_d = itertools.combinations(range(len(data.index)), sum(data.treat == 'D'))\n",
    "result = [get_mean_difference(data=data, ind=sample, var='EEAUC') for sample in all_combinations_d]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pyplot.hist(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index_d = data.loc[data.treat == 'D'].index\n",
    "mean_original = get_mean_difference(data=data, ind=index_d, var='EEAUC')\n",
    "print(f'The original statistic is: {mean_original}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(np.abs(result)>=mean_original)/len(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can get a CI for the t-statistic:\n",
    "level = 0.05\n",
    "ci = np.quantile(result, [level/2, 1-level/2])\n",
    "print(f'{100*(1-level)}% CI: ({ci[0]}, {ci[1]})')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='4'></a>\n",
    "# Montecarlo or approximated permutation test\n",
    "\n",
    "Exact permutation test as the previous one is impossible to be computed for medium or large sample sizes. The most usual solution is doing an aproximated permutation test or Montecarlo test. \n",
    "\n",
    "In this approach, a random sample of permutations is simulated, usually very large but smaller than the whole 'population' of possible permutations, which can be huge. Proceeding in this way, a *p-value* estimation is obtained which can be supposed very precise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='5'></a>\n",
    "# Permutation test and AB testing\n",
    "We decide to perform an *AB test* in order to check the performance of our new algorithm (**A group**). For this experiment, the following decisions are taken:\n",
    "- Split the population randomly with 50% of probability\n",
    "- Same sample size in both groups (100 visits)\n",
    "\n",
    "Once the experiment is finished, we take the data and analyse it. We have to decide if the new algorithm has a better performace, which means **having a higher click rate**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab_result = pd.read_csv('./data/ab.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab_result.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab_result.groupby('group')['click'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We are going to perform a permutation test to analyse if the \"labels\" are exchangeable. The equivalent parametric hypothesis test would be following one:\n",
    "\n",
    "$$\\left\\{\n",
    "\\begin{array}{ll}\n",
    "      H_0 & \\mu_A=\\mu_B \\\\\n",
    "      H_1 & \\mu_A \\ne \\mu_B \\\\\n",
    "\\end{array} \n",
    "\\right. $$\n",
    "\n",
    "Instef uf using traditional hypothesis test (*binomial test* is this case),  we are going to use permutation test. So, we are testing whether the labels are exchangeable or not. As the statistic we will use the mean difference. Notice that if we wanted to perform the exact permutation test we would need a lot of permutatitions:\n",
    "\n",
    "$$\\frac{200!}{100!\\cdot100!}.$$\n",
    "\n",
    "We are going to use *Montecarlo permutation test* with 9999 permutations. When using this approximation, an unbiased estimator of *p-value* is the following one:\n",
    "\n",
    "$$\\frac{\\text{simulations} >= \\text{original statistic} \\bf{+1}}{\\text{total_perms} \\bf{+1}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we calculate the statistic for the original sample\n",
    "ind = ab_result.loc[ab_result.group == 'A'].index\n",
    "orig_stat = get_mean_difference(data=ab_result, ind=ind, var='click')\n",
    "print(f'The original statistic is: {orig_stat}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate permutations\n",
    "total_perms = 9999\n",
    "sample_size = len(ab_result)\n",
    "len_a = sum(ab_result.group == 'A')\n",
    "ind_list = [np.random.choice(sample_size, len_a) for i in range(total_perms)]\n",
    "result = [get_mean_difference(data=ab_result, ind=ind, var='click') for ind in ind_list]\n",
    "print(f'Total permutations: {len(result)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value = (sum(np.abs(result) >=orig_stat) +1)/(len(result)+1)\n",
    "print(f'p-value: {p_value}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given such a *p-value*. Do you think you could have used result to estimate a CI for the the t-test?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pyplot.hist(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can get a CI for the t-statistic:\n",
    "level = 0.05\n",
    "ci = np.quantile(result, [level/2, 1-level/2])\n",
    "print(f'{100*(1-level)}% CI: ({ci[0]}, {ci[1]})')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name='7'></a>\n",
    "## (Brief) Introduction to Bootstrap\n",
    "\n",
    "Bootstrap method is a technique that provides a general methodology useful to determine the sample distribution of a statistic. Like permutation test, is a resampling method where **replacement** is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_boots = 1000\n",
    "a_values = ab_result.loc[ab_result.group == 'A', 'click'].values\n",
    "b_values = ab_result.loc[ab_result.group == 'B', 'click'].values\n",
    "diff_mean = []\n",
    "for i in range(n_boots):\n",
    "    ind_a = np.random.choice(len(a_values), len(a_values), replace=True)\n",
    "    ind_b = np.random.choice(len(b_values), len(b_values), replace=True)\n",
    "    a_mean = np.mean(a_values[ind_a])\n",
    "    b_mean = np.mean(b_values[ind_b])\n",
    "    diff_mean.append(a_mean-b_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pyplot.hist(diff_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can get a CI for the means difference:\n",
    "level = 0.05\n",
    "ci = np.quantile(diff_mean, [level/2, 1-level/2])\n",
    "print(f'{100*(1-level)}% CI: ({ci[0]}, {ci[1]})')"
   ]
  }
 ],
 "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": 4
}