Chapter_13_sec_13.6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 13.6 Lab: Multiple Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats as st\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "from statsmodels.sandbox.stats.multicomp import TukeyHSDResults\n",
    "\n",
    "import json\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.6.1 Review of Hypothesis Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(21)\n",
    "X = np.random.normal(loc=0.0, scale=1.0, size=(10, 100))\n",
    "offset = 0.5\n",
    "X[:,:50] = X[:,:50] + offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here I used scipy. During google search, I came across bioinfokit module, could explore more. \n",
    "result=st.ttest_1samp(a = X[:, 0], popmean = 0)\n",
    "print(result.pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us run the same t-test for all 100 columns\n",
    "p_values = []\n",
    "decision = []\n",
    "for i in range(100):\n",
    "    result=st.ttest_1samp(a = X[:, i], popmean = 0)\n",
    "    p_values.append(result.pvalue)\n",
    "    if result.pvalue < 0.05:\n",
    "        decision.append('Reject H0')\n",
    "    else:\n",
    "        decision.append('Do not reject H0')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# after computing the p-values, we can use the ground truth to evaluate the performance\n",
    "ground_truth = np.repeat(['Reject H0', 'Do not reject H0'], [50, 50], axis=0)\n",
    "labels = ['Reject H0', 'Do not reject H0']\n",
    "cm = confusion_matrix (ground_truth, decision, labels=labels)\n",
    "print(cm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "cax = ax.matshow(cm)\n",
    "fig.colorbar(cax)\n",
    "ax.set_xticklabels([''] + labels)\n",
    "ax.set_yticklabels([''] + labels)\n",
    "plt.xlabel('One sample t-test')\n",
    "plt.ylabel('Ground truth')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could make the offset larger (from 0.5 to 1) and see the change to the confusion matrix\n",
    "offset = 1\n",
    "X[:,:50] = X[:,:50] + offset\n",
    "\n",
    "p_values = []\n",
    "decision = []\n",
    "for i in range(100):\n",
    "    result=st.ttest_1samp(a = X[:, i], popmean = 0)\n",
    "    p_values.append(result.pvalue)\n",
    "    if result.pvalue < 0.05:\n",
    "        decision.append('Reject H0')\n",
    "    else:\n",
    "        decision.append('Do not reject H0')\n",
    "\n",
    "\n",
    "ground_truth = np.repeat(['Reject H0', 'Do not reject H0'], [50, 50], axis=0)\n",
    "labels = ['Reject H0', 'Do not reject H0']\n",
    "cm = confusion_matrix (ground_truth, decision, labels=labels)\n",
    "print(cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.6.2 The Family-Wise Error Rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = range(500)\n",
    "fwe1 = list(map(lambda x:1 - pow(1 - 0.05,x),m))\n",
    "fwe2 = list(map(lambda x:1 - pow(1 - 0.01,x),m))\n",
    "fwe3 = list(map(lambda x:1 - pow(1 - 0.001,x),m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(m, fwe1, label = \"0.05\")\n",
    "plt.plot(m, fwe2, label = \"0.01\")\n",
    "plt.plot(m, fwe3, label = \"0.001\")\n",
    "plt.xlabel('Number of tests in log scale')\n",
    "plt.ylabel('FWE')\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\"\"\" \n",
    "We see that setting α = 0.05 results in a high FWER even for moderate m. \n",
    "With α = 0.01, we can test no more than five null hypotheses before the FWER exceeds 0.05. \n",
    "Only for very small values, such as α = 0.001, do we manage to ensure a small FWER, \n",
    "at least for moderately-sized m.\n",
    "\n",
    "Of course, the problem with setting α to such a low value is that we are likely to \n",
    "make a number of Type II errors: in other words, our power is very low.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fund = pd.read_csv('data/Fund.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fund.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will do the one sample t test for the first manager\n",
    "result=st.ttest_1samp(a = Fund['Manager1'], popmean = 0)\n",
    "print(result.pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_values = []\n",
    "manager_number = 5 \n",
    "\n",
    "for i in range(manager_number):\n",
    "    result=st.ttest_1samp(a = Fund.iloc[:,i], popmean = 0)\n",
    "    p_values.append(result.pvalue)\n",
    "\n",
    "print(p_values)\n",
    "\n",
    "\"\"\" \n",
    "The p-values are low for Managers One and Three, and high for the other three managers. \n",
    "However, we cannot simply reject H01 and H03, since this would fail to account for \n",
    "the multiple testing that we have performed. \n",
    "Instead, we will conduct Bonferroni’s method and Holm’s method to control the FWER.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could bonferroni to adjust the raw p-values and take care of family wise error rate\n",
    "reject, p_values_corrected, alphacSidak, alphacBonf = multipletests(p_values, method = 'bonferroni')\n",
    "print(p_values_corrected)\n",
    "\"\"\" \n",
    "Therefore, using Bonferroni’s method, \n",
    "we are able to reject the null hypothesis only for Manager One while controlling the FWER at 0.05.\n",
    "This information is also available in the variable reject.\n",
    "\"\"\"\n",
    "print(reject)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonferroni's method is more conservative. We could apply holm's method to control the FWER\n",
    "reject, p_values_corrected, alphacSidak, alphacBonf = multipletests(p_values, method = 'holm')\n",
    "print(p_values_corrected)\n",
    "print(reject)\n",
    "\"\"\" \n",
    "By contrast, using Holm’s method, the adjusted p-values indicate that we can reject the null hypotheses \n",
    "for Both Managers One and Three at a FWER of 0.05.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can see the average for each manager \n",
    "Fund.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# next, we could test whether 2 managers are significantly different. For example Manager 1 and Manager 2\n",
    "result=st.ttest_rel(a = Fund['Manager1'], b = Fund['Manager2'])\n",
    "print(result.pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "However, we decided to perform this test only after examining the data and \n",
    "noting that Managers One and Two had the highest and lowest mean performances. \n",
    "In a sense, this means that we have implicitly performed a manual selection \n",
    "from the 5(5 − 1)/2 = 10 hypothesis tests, rather than just one. \n",
    "Hence, we use Tukey’s method in order to adjust for multiple testing. \n",
    "\"\"\"\n",
    "returns = Fund.iloc[:, :5].to_numpy().flatten(order='F') # we flatten by col (i.e. order='F')\n",
    "manager = np.repeat(['1', '2', '3', '4', '5'], repeats=Fund.shape[0])\n",
    "\n",
    "# perform Tukey's test\n",
    "tukey = pairwise_tukeyhsd(endog=returns, groups=manager, alpha=0.05)\n",
    "\n",
    "print(tukey)\n",
    "\n",
    "\"\"\" \n",
    "Notice that the p-value for the difference between Managers One and Two has increased from 0.038 to 0.186, \n",
    "so there is no longer clear evidence of a difference between the managers’ performances.\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.6.3 The False Discovery Rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_values = []\n",
    "manager_number = Fund.shape[1]\n",
    "\n",
    "for i in range(manager_number):\n",
    "    result=st.ttest_1samp(a = Fund.iloc[:,i], popmean = 0)\n",
    "    p_values.append(result.pvalue)\n",
    "\n",
    "print(p_values[0:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "There are far too many managers to consider trying to control the FWER. \n",
    "Instead, we focus on controlling the FDR: that is, the expected fraction of rejected null \n",
    "hypotheses that are actually false positives. \n",
    "\"\"\"\n",
    "\n",
    "reject, p_values_corrected, alphacSidak, alphacBonf = multipletests(p_values, method = 'fdr_bh')\n",
    "print(p_values_corrected[0:10])\n",
    "\n",
    "\"\"\" \n",
    "The q-values output by the Benjamini-Hochberg procedure can be interpreted as the smallest \n",
    "FDR threshold at which we would reject a particular null hypothesis.\n",
    "\n",
    "For instance, a q-value of 0.1 indicates that we can reject the corresponding null hypothesis\n",
    "at an FDR of 10% or greater, but that we cannot reject the null hypothesis at an FDR below 10%.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we would find that 146 of the 2,000 fund managers have a p_values_corrected below 0.1\n",
    "sum(p_values_corrected <= .1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if we use bonferroni method, we will find None\n",
    "sum(np.array(p_values) <= .1/Fund.shape[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.6.4 A Re-Sampling Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I saved the gene expression data as a json file, in python we could load the json file using the json library\n",
    "# after reading in the data, we can use the data is same as a dictionary, we can use the keys to access the data\n",
    "\n",
    "f = open('./data/Khan.json')\n",
    "Khan = json.load(f)\n",
    "\n",
    "X_train = np.array(Khan['xtrain'])\n",
    "y_train = np.array(Khan['ytrain'])\n",
    "X_test = np.array(Khan['xtest'])\n",
    "y_test = np.array(Khan['ytest'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.concatenate((X_train, X_test), axis=0)\n",
    "y = np.concatenate((y_train, y_test), axis=0)\n",
    "unique, counts = np.unique(y, return_counts=True)\n",
    "print(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x1: take the x for cancer type == 2\n",
    "# x2: take the x for cancer type == 4\n",
    "x1 = x[y==2, :]\n",
    "x2 = x[y==4, :]\n",
    "n1 = x1.shape[0]\n",
    "n2 = x2.shape[0]\n",
    "print(n1)\n",
    "print(n2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# performing a standard two-sample t-test on the 11th (gene_index = 10 in python) gene produces a test-statistic \n",
    "gene_index = 10\n",
    "original_result=st.ttest_ind(a=x1[:,gene_index], b=x2[:,gene_index], equal_var=True)\n",
    "print(original_result.statistic)\n",
    "print(original_result.pvalue)\n",
    "\n",
    "\"\"\" \n",
    "The 2 sample t-test produces a test-statistic of −2.09 and an associated p-value of 0.0412, \n",
    "suggesting modest evidence of a difference in mean expression levels between the two cancer types.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Instead of doing a parameterized 2 sample t-test, we could do a non-parameterized test(i.e. permutation test).\n",
    "we can randomly split the 54 patients (in cancer group 2 and 4) into two groups of 29 and 25 \n",
    "(same as the original split),and compute a new test statistic. \n",
    "Under the null hypothesis of no difference between the groups, this new test statistic should have \n",
    "the same distribution as our original one. \n",
    "Repeating this process many (i.e.10,000) times allows us to approximate the null distribution of the test statistic. \n",
    "We compute the fraction of the time that our observed test statistic exceeds the test statistics obtained \n",
    "via re-sampling.\n",
    "\"\"\"\n",
    "\n",
    "np.random.seed(21)\n",
    "iteration = 10000\n",
    "test_stats = []\n",
    "x_temp = np.concatenate((x1[:,gene_index], x2[:,gene_index]), axis=0)\n",
    "\n",
    "for i in range(iteration):\n",
    "    np.random.shuffle(x_temp)\n",
    "    result_temp = st.ttest_ind(a=x_temp[:n1], b=x_temp[-n2:], equal_var=True)\n",
    "    test_stats.append(result_temp.statistic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean((np.abs(test_stats) >= np.abs(original_result.statistic))))\n",
    "\n",
    "\"\"\" \n",
    "This fraction is our re-sampling-based p-value. It is almost identical to the p-value of 0.0412 \n",
    "obtained using the theoretical null distribution.\n",
    "\n",
    "The reason for this is that the parametrized distribution is a pretty good assumption in this case\n",
    "To see this, we can plot the histogram of the re-sampled statistics vs. parametrized distribution. \n",
    "\n",
    "We could try other genes (i.e. gene_index = 876) to see its theoretical and re-sampling null distributions are \n",
    "quite different\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct the t distribution \n",
    "df = n1 + n2 - 2\n",
    "rv = st.t(df)\n",
    "x = np.linspace(-4.2, 4.2, 1000)\n",
    "\n",
    "\n",
    "plt.hist(test_stats, 100, density=True, facecolor='g', alpha=0.75)\n",
    "plt.plot(x, rv.pdf(x), 'k-', lw=2)\n",
    "plt.xlabel('Null Distribution of Test Statistic')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Histogram of re-sample stats')\n",
    "plt.xlim(-4.2, 4.2)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could do this for 100 and see how FDR works under re-sample \n",
    "# it would be good to use small iterations to make sure the code runs okay \n",
    "num_gene = 100\n",
    "iteration = 500\n",
    "test_stats_matrix = []\n",
    "test_stats_origin = []\n",
    "\n",
    "for j in range(num_gene):\n",
    "    gene_index = j \n",
    "    x_temp = np.concatenate((x1[:,gene_index], x2[:,gene_index]), axis=0)\n",
    "    result_origin = st.ttest_ind(a=x1[:,gene_index], b=x2[:,gene_index], equal_var=True)\n",
    "    test_stats_origin.append(result_origin.statistic)\n",
    "    test_stats = []\n",
    "    for i in range(iteration):\n",
    "        np.random.shuffle(x_temp)\n",
    "        result_temp = st.ttest_ind(a=x_temp[:n1], b=x_temp[-n2:], equal_var=True)\n",
    "        test_stats.append(result_temp.statistic)\n",
    "        \n",
    "    test_stats_matrix.append(test_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_stats_origin_sorted =  np.sort(np.abs(test_stats_origin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs = []\n",
    "Vs = []\n",
    "FDRs = []\n",
    "for j in range(num_gene):\n",
    "    R = np.sum(np.abs(test_stats_origin) >= test_stats_origin_sorted[j])\n",
    "    V = np.sum(np.abs(test_stats_matrix) >= test_stats_origin_sorted[j]) / iteration\n",
    "    Rs.append(R)\n",
    "    Vs.append(V)\n",
    "    FDRs.append(V*1.0/R)\n",
    "\n",
    "Rs = np.array(Rs)\n",
    "Vs = np.array(Vs)\n",
    "FDRs = np.array(FDRs) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(Rs[FDRs <= .1]))\n",
    "print(np.max(Rs[FDRs <= .2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Rs, FDRs, 'k-', lw=2)\n",
    "plt.xlabel('Number of Rejections')\n",
    "plt.ylabel('False Discovery Rate')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 13"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}