{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75b2d75c",
   "metadata": {},
   "source": [
    "\n",
    "# Chapter 13\n",
    "\n",
    "# Lab: Multiple Testing\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34e410a6",
   "metadata": {},
   "source": [
    "We include our usual imports seen in earlier labs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1f928b2d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:21.791802Z",
     "iopub.status.busy": "2023-07-26T10:52:21.791660Z",
     "iopub.status.idle": "2023-07-26T10:52:23.360588Z",
     "shell.execute_reply": "2023-07-26T10:52:23.359991Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from ISLP import load_data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12319e0a",
   "metadata": {},
   "source": [
    "We also collect the new imports\n",
    "needed for this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eb4b32aa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.363356Z",
     "iopub.status.busy": "2023-07-26T10:52:23.362843Z",
     "iopub.status.idle": "2023-07-26T10:52:23.365980Z",
     "shell.execute_reply": "2023-07-26T10:52:23.365438Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "from scipy.stats import \\\n",
    "    (ttest_1samp,\n",
    "     ttest_rel,\n",
    "     ttest_ind,\n",
    "     t as t_dbn)\n",
    "from statsmodels.stats.multicomp import \\\n",
    "     pairwise_tukeyhsd\n",
    "from statsmodels.stats.multitest import \\\n",
    "     multipletests as mult_test\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2747e58",
   "metadata": {},
   "source": [
    "## Review of Hypothesis Tests\n",
    "We begin by performing some one-sample $t$-tests.\n",
    "\n",
    "First we create 100 variables, each consisting of 10 observations. The\n",
    "first 50 variables have mean $0.5$ and variance $1$, while the others\n",
    "have mean $0$ and variance $1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e12ac0cd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.368340Z",
     "iopub.status.busy": "2023-07-26T10:52:23.368164Z",
     "iopub.status.idle": "2023-07-26T10:52:23.371202Z",
     "shell.execute_reply": "2023-07-26T10:52:23.370779Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(12)\n",
    "X = rng.standard_normal((10, 100))\n",
    "true_mean = np.array([0.5]*50 + [0]*50)\n",
    "X += true_mean[None,:]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70d37233",
   "metadata": {},
   "source": [
    "To begin, we use `ttest_1samp()`  from the\n",
    "`scipy.stats` module to test $H_{0}: \\mu_1=0$, the null\n",
    "hypothesis that the first variable has mean zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "04d0f49e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.373512Z",
     "iopub.status.busy": "2023-07-26T10:52:23.373319Z",
     "iopub.status.idle": "2023-07-26T10:52:23.378383Z",
     "shell.execute_reply": "2023-07-26T10:52:23.377960Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9307442156164141"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = ttest_1samp(X[:,0], 0)\n",
    "result.pvalue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf83426f",
   "metadata": {},
   "source": [
    "The $p$-value comes out to 0.931, which is not low enough to\n",
    "reject the null hypothesis at level $\\alpha=0.05$.  In this case,\n",
    "$\\mu_1=0.5$, so the null hypothesis is false. Therefore, we have made\n",
    "a Type II error by failing to reject the null hypothesis when the null\n",
    "hypothesis is false. \n",
    "\n",
    "We now test $H_{0,j}: \\mu_j=0$ for $j=1,\\ldots,100$. We compute the\n",
    "100 $p$-values, and then construct a vector recording whether the\n",
    "$j$th $p$-value is less than or equal to 0.05, in which case we reject\n",
    "$H_{0j}$, or greater than 0.05, in which case we do not reject\n",
    "$H_{0j}$, for $j=1,\\ldots,100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d1f0c695",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.380819Z",
     "iopub.status.busy": "2023-07-26T10:52:23.380641Z",
     "iopub.status.idle": "2023-07-26T10:52:23.408797Z",
     "shell.execute_reply": "2023-07-26T10:52:23.408278Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "p_values = np.empty(100)\n",
    "for i in range(100):\n",
    "   p_values[i] = ttest_1samp(X[:,i], 0).pvalue\n",
    "decision = pd.cut(p_values,\n",
    "                  [0, 0.05, 1],\n",
    "                  labels=['Reject H0',\n",
    "                          'Do not reject H0'])\n",
    "truth = pd.Categorical(true_mean == 0,\n",
    "                       categories=[True, False],\n",
    "                       ordered=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d8e0d96",
   "metadata": {},
   "source": [
    "Since this is a simulated data set, we can create a $2 \\times 2$ table\n",
    "similar to  Table 13.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7a9594a0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.411256Z",
     "iopub.status.busy": "2023-07-26T10:52:23.411053Z",
     "iopub.status.idle": "2023-07-26T10:52:23.422652Z",
     "shell.execute_reply": "2023-07-26T10:52:23.422286Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>H0</th>\n",
       "      <th>True</th>\n",
       "      <th>False</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decision</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Reject H0</th>\n",
       "      <td>5</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Do not reject H0</th>\n",
       "      <td>45</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "H0                True  False\n",
       "Decision                     \n",
       "Reject H0            5     15\n",
       "Do not reject H0    45     35"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(decision,\n",
    "            truth,\n",
    "     rownames=['Decision'],\n",
    "     colnames=['H0'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9610c817",
   "metadata": {},
   "source": [
    "Therefore, at level $\\alpha=0.05$, we reject 15 of the 50 false\n",
    "null hypotheses, and we incorrectly reject 5 of the true null\n",
    "hypotheses. Using the notation from Section 13.3, we have\n",
    "$V=5$, $S=15$, $U=45$ and $W=35$.\n",
    "We have set $\\alpha=0.05$, which means that we expect to reject around\n",
    "5% of the true null hypotheses. This is in line with the $2 \\times 2$\n",
    "table above, which indicates that we rejected $V=5$ of the $50$ true\n",
    "null hypotheses.\n",
    "\n",
    "In the simulation above, for the false null hypotheses, the ratio of\n",
    "the mean to the standard deviation was only $0.5/1 = 0.5$. This\n",
    "amounts to quite a weak signal, and it resulted in a high number of\n",
    "Type II errors.  Let’s instead simulate data with a stronger signal,\n",
    "so that the ratio of the mean to the standard deviation for the false\n",
    "null hypotheses equals $1$. We make only 10 Type II errors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "25f7fc5d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.425157Z",
     "iopub.status.busy": "2023-07-26T10:52:23.424798Z",
     "iopub.status.idle": "2023-07-26T10:52:23.455752Z",
     "shell.execute_reply": "2023-07-26T10:52:23.455052Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>H0</th>\n",
       "      <th>True</th>\n",
       "      <th>False</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decision</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Reject H0</th>\n",
       "      <td>2</td>\n",
       "      <td>40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Do not reject H0</th>\n",
       "      <td>48</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "H0                True  False\n",
       "Decision                     \n",
       "Reject H0            2     40\n",
       "Do not reject H0    48     10"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_mean = np.array([1]*50 + [0]*50)\n",
    "X = rng.standard_normal((10, 100))\n",
    "X += true_mean[None,:]\n",
    "for i in range(100):\n",
    "   p_values[i] = ttest_1samp(X[:,i], 0).pvalue\n",
    "decision = pd.cut(p_values,\n",
    "                  [0, 0.05, 1],\n",
    "                  labels=['Reject H0',\n",
    "                          'Do not reject H0'])\n",
    "truth = pd.Categorical(true_mean == 0,\n",
    "                       categories=[True, False],\n",
    "                       ordered=True)\n",
    "pd.crosstab(decision,\n",
    "            truth,\n",
    "            rownames=['Decision'],\n",
    "            colnames=['H0'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb70c597",
   "metadata": {},
   "source": [
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6953d33",
   "metadata": {},
   "source": [
    "## Family-Wise Error Rate\n",
    "Recall from  (13.5)  that if the null hypothesis is true\n",
    "for each of $m$ independent hypothesis tests, then the FWER is equal\n",
    "to $1-(1-\\alpha)^m$. We can use this expression to compute the FWER\n",
    "for $m=1,\\ldots, 500$ and $\\alpha=0.05$, $0.01$, and $0.001$.\n",
    "We plot the FWER for these values of $\\alpha$ in order to\n",
    "reproduce  Figure 13.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "369b5bd3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.458631Z",
     "iopub.status.busy": "2023-07-26T10:52:23.458434Z",
     "iopub.status.idle": "2023-07-26T10:52:23.710141Z",
     "shell.execute_reply": "2023-07-26T10:52:23.709635Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43736 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = np.linspace(1, 501)\n",
    "fig, ax = plt.subplots()\n",
    "[ax.plot(m,\n",
    "         1 - (1 - alpha)**m,\n",
    "         label=r'$\\alpha=%s$' % str(alpha))\n",
    "         for alpha in [0.05, 0.01, 0.001]]\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('Number of Hypotheses')\n",
    "ax.set_ylabel('Family-Wise Error Rate')\n",
    "ax.legend()\n",
    "ax.axhline(0.05, c='k', ls='--');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a81479e",
   "metadata": {},
   "source": [
    "As discussed previously, even for moderate values of $m$ such as $50$,\n",
    "the FWER exceeds $0.05$ unless $\\alpha$ is set to a very low value,\n",
    "such as $0.001$.  Of course, the problem with setting $\\alpha$ to such\n",
    "a low value is that we are likely to make a number of Type II errors:\n",
    "in other words, our power is very low.\n",
    "\n",
    "We now conduct a one-sample $t$-test for each of the first five\n",
    "managers in the  \n",
    "`Fund`   dataset, in order to test the null\n",
    "hypothesis that the $j$th fund manager’s mean return equals zero,\n",
    "$H_{0,j}: \\mu_j=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9ce7a19f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.713438Z",
     "iopub.status.busy": "2023-07-26T10:52:23.713214Z",
     "iopub.status.idle": "2023-07-26T10:52:23.754484Z",
     "shell.execute_reply": "2023-07-26T10:52:23.754038Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00620236, 0.91827115, 0.01160098, 0.6005396 , 0.75578151])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fund = load_data('Fund')\n",
    "fund_mini = Fund.iloc[:,:5]\n",
    "fund_mini_pvals = np.empty(5)\n",
    "for i in range(5):\n",
    "    fund_mini_pvals[i] = ttest_1samp(fund_mini.iloc[:,i], 0).pvalue\n",
    "fund_mini_pvals\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7561e3a3",
   "metadata": {},
   "source": [
    "The $p$-values are low for Managers One and Three, and high for the\n",
    "other three managers.  However, we cannot simply reject $H_{0,1}$ and\n",
    "$H_{0,3}$, since this would fail to account for the multiple testing\n",
    "that we have performed. Instead, we will conduct Bonferroni’s method\n",
    "and Holm’s method to control the FWER.\n",
    "\n",
    "To do this, we use the `multipletests()`  function from the\n",
    "`statsmodels` module (abbreviated to `mult_test()`). Given the $p$-values,\n",
    "for methods like Holm and Bonferroni the function outputs\n",
    "adjusted $p$-values,  which\n",
    "can be thought of as a new set of $p$-values that have been corrected\n",
    "for multiple testing. If the adjusted $p$-value for a given hypothesis\n",
    "is less than or equal to $\\alpha$, then that hypothesis can be\n",
    "rejected while maintaining a FWER of no more than $\\alpha$. In other\n",
    "words, for such methods, the adjusted $p$-values resulting from the `multipletests()`\n",
    "function can simply be compared to the desired FWER in order to\n",
    "determine whether or not to reject each hypothesis. We will later\n",
    "see that we can use the same function to control FDR as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b608e46",
   "metadata": {},
   "source": [
    "The `mult_test()` function takes $p$-values and a `method` argument, as well as an optional\n",
    "`alpha` argument. It returns the  decisions (`reject` below)\n",
    "as well as the adjusted $p$-values (`bonf`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "de6cffed",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.756665Z",
     "iopub.status.busy": "2023-07-26T10:52:23.756493Z",
     "iopub.status.idle": "2023-07-26T10:52:23.759438Z",
     "shell.execute_reply": "2023-07-26T10:52:23.759028Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True, False, False, False, False])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reject, bonf = mult_test(fund_mini_pvals, method = \"bonferroni\")[:2]\n",
    "reject\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5135c6b9",
   "metadata": {},
   "source": [
    "The $p$-values `bonf` are simply the `fund_mini_pvalues` multiplied by 5 and truncated to be less than\n",
    "or equal to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0de71500",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.761459Z",
     "iopub.status.busy": "2023-07-26T10:52:23.761287Z",
     "iopub.status.idle": "2023-07-26T10:52:23.764317Z",
     "shell.execute_reply": "2023-07-26T10:52:23.763952Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.03101178, 1.        , 0.05800491, 1.        , 1.        ]),\n",
       " array([0.03101178, 1.        , 0.05800491, 1.        , 1.        ]))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bonf, np.minimum(fund_mini_pvals * 5, 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f0bc112",
   "metadata": {},
   "source": [
    "Therefore, using Bonferroni’s method, we are able to reject the null hypothesis only for Manager\n",
    "One while controlling FWER at $0.05$.\n",
    "\n",
    "By contrast, using Holm’s method, the adjusted $p$-values indicate\n",
    "that we can  reject the null\n",
    "hypotheses for Managers One and Three at a FWER of $0.05$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f7e87bdb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.766793Z",
     "iopub.status.busy": "2023-07-26T10:52:23.766606Z",
     "iopub.status.idle": "2023-07-26T10:52:23.810152Z",
     "shell.execute_reply": "2023-07-26T10:52:23.809736Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ True, False,  True, False, False]),\n",
       " array([0.03101178, 1.        , 0.04640393, 1.        , 1.        ]))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mult_test(fund_mini_pvals, method = \"holm\", alpha=0.05)[:2]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f762fecd",
   "metadata": {},
   "source": [
    "As discussed previously, Manager One seems to perform particularly\n",
    "well, whereas Manager Two has poor performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e88be376",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.812384Z",
     "iopub.status.busy": "2023-07-26T10:52:23.812232Z",
     "iopub.status.idle": "2023-07-26T10:52:23.815971Z",
     "shell.execute_reply": "2023-07-26T10:52:23.815576Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Manager1    3.0\n",
       "Manager2   -0.1\n",
       "Manager3    2.8\n",
       "Manager4    0.5\n",
       "Manager5    0.3\n",
       "dtype: float64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fund_mini.mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88dbf0a6",
   "metadata": {},
   "source": [
    "Is there evidence of a meaningful difference in performance between\n",
    "these two managers?  We can check this by performing a  paired $t$-test  using the `ttest_rel()` function\n",
    "from `scipy.stats`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "41149af6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.817997Z",
     "iopub.status.busy": "2023-07-26T10:52:23.817852Z",
     "iopub.status.idle": "2023-07-26T10:52:23.821301Z",
     "shell.execute_reply": "2023-07-26T10:52:23.820976Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.038391072368079586"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ttest_rel(fund_mini['Manager1'],\n",
    "          fund_mini['Manager2']).pvalue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aca6122",
   "metadata": {},
   "source": [
    "The test results in a $p$-value of 0.038,\n",
    "suggesting a statistically significant difference.\n",
    "\n",
    "However, we decided to perform this test only after examining the data\n",
    "and noting that Managers One and Two had the highest and lowest mean\n",
    "performances.  In a sense, this means that we have implicitly\n",
    "performed ${5 \\choose 2} = 5(5-1)/2=10$ hypothesis tests, rather than\n",
    "just one, as discussed in  Section 13.3.2.  Hence, we use the\n",
    "`pairwise_tukeyhsd()`  function from\n",
    "`statsmodels.stats.multicomp` to apply Tukey’s method\n",
    "  in order to adjust for multiple testing.  This function takes\n",
    "as input a fitted *ANOVA*  regression model, which is\n",
    "essentially just a linear regression in which all of the predictors\n",
    "are qualitative.  In this case, the response consists of the monthly\n",
    "excess returns achieved by each manager, and the predictor indicates\n",
    "the manager to which each return corresponds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "61aabda7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:23.823367Z",
     "iopub.status.busy": "2023-07-26T10:52:23.823190Z",
     "iopub.status.idle": "2023-07-26T10:52:24.308715Z",
     "shell.execute_reply": "2023-07-26T10:52:24.308222Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiple Comparison of Means - Tukey HSD, FWER=0.05\n",
      "===================================================\n",
      "group1 group2 meandiff p-adj   lower  upper  reject\n",
      "---------------------------------------------------\n",
      "     1      2     -3.1 0.1862 -6.9865 0.7865  False\n",
      "     1      3     -0.2 0.9999 -4.0865 3.6865  False\n",
      "     1      4     -2.5 0.3948 -6.3865 1.3865  False\n",
      "     1      5     -2.7 0.3152 -6.5865 1.1865  False\n",
      "     2      3      2.9 0.2453 -0.9865 6.7865  False\n",
      "     2      4      0.6 0.9932 -3.2865 4.4865  False\n",
      "     2      5      0.4 0.9986 -3.4865 4.2865  False\n",
      "     3      4     -2.3  0.482 -6.1865 1.5865  False\n",
      "     3      5     -2.5 0.3948 -6.3865 1.3865  False\n",
      "     4      5     -0.2 0.9999 -4.0865 3.6865  False\n",
      "---------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "returns = np.hstack([fund_mini.iloc[:,i] for i in range(5)])\n",
    "managers = np.hstack([[i+1]*50 for i in range(5)])\n",
    "tukey = pairwise_tukeyhsd(returns, managers)\n",
    "print(tukey.summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0084fc5",
   "metadata": {},
   "source": [
    "The `pairwise_tukeyhsd()` function provides confidence intervals\n",
    "for the difference between each pair of managers (`lower` and\n",
    "`upper`), as well as a $p$-value. All of these quantities have\n",
    "been adjusted for multiple testing. Notice that the $p$-value for the\n",
    "difference between Managers One and Two has increased from $0.038$ to\n",
    "$0.186$, so there is no longer clear evidence of a difference between\n",
    "the managers’ performances.  We can plot the confidence intervals for\n",
    "the pairwise comparisons using the `plot_simultaneous()` method\n",
    "of `tukey`. Any pair of intervals that don’t overlap indicates a significant difference at the nominal level of 0.05. In this case,\n",
    "no differences are considered significant as reported in the table above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cbcad4de",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.311170Z",
     "iopub.status.busy": "2023-07-26T10:52:24.310960Z",
     "iopub.status.idle": "2023-07-26T10:52:24.410742Z",
     "shell.execute_reply": "2023-07-26T10:52:24.410212Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21100 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "tukey.plot_simultaneous(ax=ax);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6278d13c",
   "metadata": {},
   "source": [
    "## False Discovery Rate\n",
    "Now we perform hypothesis tests for all 2,000 fund managers in the\n",
    "`Fund`  dataset. We perform a one-sample $t$-test\n",
    "of   $H_{0,j}: \\mu_j=0$, which states that the\n",
    "$j$th fund manager’s mean return is zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b5842190",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.413638Z",
     "iopub.status.busy": "2023-07-26T10:52:24.413455Z",
     "iopub.status.idle": "2023-07-26T10:52:24.827343Z",
     "shell.execute_reply": "2023-07-26T10:52:24.826609Z"
    }
   },
   "outputs": [],
   "source": [
    "fund_pvalues = np.empty(2000)\n",
    "for i, manager in enumerate(Fund.columns):\n",
    "    fund_pvalues[i] = ttest_1samp(Fund[manager], 0).pvalue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80fc2fcc",
   "metadata": {},
   "source": [
    "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 hypotheses that are actually false positives.\n",
    "The `multipletests()` function (abbreviated `mult_test()`) can be used to carry out the Benjamini--Hochberg procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7c9d8bed",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.830214Z",
     "iopub.status.busy": "2023-07-26T10:52:24.830015Z",
     "iopub.status.idle": "2023-07-26T10:52:24.834953Z",
     "shell.execute_reply": "2023-07-26T10:52:24.834464Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.08988921, 0.991491  , 0.12211561, 0.92342997, 0.95603587,\n",
       "       0.07513802, 0.0767015 , 0.07513802, 0.07513802, 0.07513802])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fund_qvalues = mult_test(fund_pvalues, method = \"fdr_bh\")[1]\n",
    "fund_qvalues[:10]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f73096d",
   "metadata": {},
   "source": [
    "The  *q-values* output by the\n",
    "Benjamini--Hochberg procedure can be interpreted as the smallest FDR\n",
    "threshold at which we would reject a particular null hypothesis. For\n",
    "instance, a $q$-value of $0.1$ indicates that we can reject the\n",
    "corresponding null hypothesis at an FDR of 10% or greater, but that\n",
    "we cannot reject the null hypothesis at an FDR below 10%.\n",
    "\n",
    "If we control the FDR at 10%, then for how many of the fund managers can we reject $H_{0,j}: \\mu_j=0$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bfa39f7c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.837627Z",
     "iopub.status.busy": "2023-07-26T10:52:24.837419Z",
     "iopub.status.idle": "2023-07-26T10:52:24.840834Z",
     "shell.execute_reply": "2023-07-26T10:52:24.840316Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "146"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(fund_qvalues <= 0.1).sum()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb44c8d",
   "metadata": {},
   "source": [
    "We find that 146 of the 2,000 fund managers have a $q$-value below\n",
    "0.1; therefore, we are able to conclude that 146 of the fund managers\n",
    "beat the market at an FDR of 10%.  Only about 15 (10% of 146) of\n",
    "these fund managers are likely to be false discoveries.\n",
    "\n",
    "By contrast, if we had instead used Bonferroni’s method to control the\n",
    "FWER at level $\\alpha=0.1$, then we would have failed to reject any\n",
    "null hypotheses!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "70b69b47",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.843470Z",
     "iopub.status.busy": "2023-07-26T10:52:24.843237Z",
     "iopub.status.idle": "2023-07-26T10:52:24.846889Z",
     "shell.execute_reply": "2023-07-26T10:52:24.846276Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(fund_pvalues <= 0.1 / 2000).sum()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a969f4",
   "metadata": {},
   "source": [
    "Figure 13.6 displays the ordered\n",
    "$p$-values, $p_{(1)} \\leq p_{(2)} \\leq \\cdots \\leq p_{(2000)}$, for\n",
    "the  `Fund`  dataset, as well as the threshold for rejection by the\n",
    "Benjamini--Hochberg procedure.  Recall that the Benjamini--Hochberg\n",
    "procedure identifies the largest $p$-value such that $p_{(j)}<qj/m$,\n",
    "and rejects all hypotheses for which the $p$-value is less than or\n",
    "equal to $p_{(j)}$. In the code below, we implement the\n",
    "Benjamini--Hochberg procedure ourselves, in order to illustrate how it\n",
    "works. We first order the $p$-values. We then identify all $p$-values\n",
    "that satisfy $p_{(j)}<qj/m$ (`sorted_set_`).  Finally, `selected_`\n",
    "is a boolean array indicating which $p$-values\n",
    " are less than or equal to the largest\n",
    "$p$-value in `sorted_[sorted_set_]`. Therefore, `selected_` indexes the\n",
    "$p$-values rejected by the Benjamini--Hochberg procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4c0ddea1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.849388Z",
     "iopub.status.busy": "2023-07-26T10:52:24.849238Z",
     "iopub.status.idle": "2023-07-26T10:52:24.852847Z",
     "shell.execute_reply": "2023-07-26T10:52:24.852354Z"
    }
   },
   "outputs": [],
   "source": [
    "sorted_ = np.sort(fund_pvalues)\n",
    "m = fund_pvalues.shape[0]\n",
    "q = 0.1\n",
    "sorted_set_ = np.where(sorted_ < q * np.linspace(1, m, m) / m)[0]\n",
    "if sorted_set_.shape[0] > 0:\n",
    "    selected_ = fund_pvalues < sorted_[sorted_set_].max()\n",
    "    sorted_set_ = np.arange(sorted_set_.max())\n",
    "else:\n",
    "    selected_ = []\n",
    "    sorted_set_ = []\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddeb3900",
   "metadata": {},
   "source": [
    "We now reproduce  the middle panel of Figure 13.6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0314eac9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:24.855721Z",
     "iopub.status.busy": "2023-07-26T10:52:24.855334Z",
     "iopub.status.idle": "2023-07-26T10:52:25.156072Z",
     "shell.execute_reply": "2023-07-26T10:52:25.155438Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28100 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(np.arange(0, sorted_.shape[0]) + 1,\n",
    "           sorted_, s=10)\n",
    "ax.set_yscale('log')\n",
    "ax.set_xscale('log')\n",
    "ax.set_ylabel('P-Value')\n",
    "ax.set_xlabel('Index')\n",
    "ax.scatter(sorted_set_+1, sorted_[sorted_set_], c='r', s=20)\n",
    "ax.axline((0, 0), (1,q/m), c='k', ls='--', linewidth=3);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83416f4a",
   "metadata": {},
   "source": [
    "## A Re-Sampling Approach\n",
    "Here, we implement the re-sampling approach to hypothesis testing\n",
    "using the  `Khan`  dataset, which we investigated in\n",
    "Section 13.5.  First, we merge the training and\n",
    "testing data, which results in observations on 83 patients for\n",
    "2,308 genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b59b8137",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:25.158998Z",
     "iopub.status.busy": "2023-07-26T10:52:25.158769Z",
     "iopub.status.idle": "2023-07-26T10:52:25.243020Z",
     "shell.execute_reply": "2023-07-26T10:52:25.242512Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2    29\n",
       "4    25\n",
       "3    18\n",
       "1    11\n",
       "Name: Y, dtype: int64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Khan = load_data('Khan')      \n",
    "D = pd.concat([Khan['xtrain'], Khan['xtest']])\n",
    "D['Y'] = pd.concat([Khan['ytrain'], Khan['ytest']])\n",
    "D['Y'].value_counts()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5534c8d4",
   "metadata": {},
   "source": [
    "There are four classes of cancer. For each gene, we compare the mean\n",
    "expression in the second class (rhabdomyosarcoma) to the mean\n",
    "expression in the fourth class (Burkitt’s lymphoma).  Performing a\n",
    "standard two-sample $t$-test  \n",
    "using `ttest_ind()`  from `scipy.stats` on the $11$th\n",
    "gene produces a test-statistic of -2.09 and an associated $p$-value\n",
    "of 0.0412, suggesting modest evidence of a difference in mean\n",
    "expression levels between the two cancer types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "96fb2f61",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:25.245359Z",
     "iopub.status.busy": "2023-07-26T10:52:25.245172Z",
     "iopub.status.idle": "2023-07-26T10:52:25.250139Z",
     "shell.execute_reply": "2023-07-26T10:52:25.249644Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.0936330736768185, 0.04118643782678394)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D2 = D[lambda df:df['Y'] == 2]\n",
    "D4 = D[lambda df:df['Y'] == 4]\n",
    "gene_11 = 'G0011'\n",
    "observedT, pvalue = ttest_ind(D2[gene_11],\n",
    "                              D4[gene_11],\n",
    "                              equal_var=True)\n",
    "observedT, pvalue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3131124e",
   "metadata": {},
   "source": [
    "However, this $p$-value relies on the assumption that under the null\n",
    "hypothesis of no difference between the two groups, the test statistic\n",
    "follows a $t$-distribution with $29+25-2=52$ degrees of freedom.\n",
    "Instead of using this theoretical null distribution, we can randomly\n",
    "split the 54 patients into two groups of 29 and 25, and compute a new\n",
    "test statistic.  Under the null hypothesis of no difference between\n",
    "the groups, this new test statistic should have the same distribution\n",
    "as our original one.  Repeating this process 10,000 times allows us to\n",
    "approximate the null distribution of the test statistic.  We compute\n",
    "the fraction of the time that our observed test statistic exceeds the\n",
    "test statistics obtained via re-sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "fdc229fa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:25.252510Z",
     "iopub.status.busy": "2023-07-26T10:52:25.252319Z",
     "iopub.status.idle": "2023-07-26T10:52:26.208474Z",
     "shell.execute_reply": "2023-07-26T10:52:26.207894Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0398"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = 10000\n",
    "Tnull = np.empty(B)\n",
    "D_ = np.hstack([D2[gene_11], D4[gene_11]])\n",
    "n_ = D2[gene_11].shape[0]\n",
    "D_null = D_.copy()\n",
    "for b in range(B):\n",
    "    rng.shuffle(D_null)\n",
    "    ttest_ = ttest_ind(D_null[:n_],\n",
    "                       D_null[n_:],\n",
    "                       equal_var=True)\n",
    "    Tnull[b] = ttest_.statistic\n",
    "(np.abs(Tnull) > np.abs(observedT)).mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7fc4557",
   "metadata": {},
   "source": [
    "This fraction, 0.0398,\n",
    "is our re-sampling-based $p$-value.\n",
    "It is almost identical to the $p$-value of  0.0412 obtained using the theoretical null distribution.\n",
    "We can plot  a histogram of the re-sampling-based test statistics in order to reproduce  Figure 13.7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e3894695",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:26.210995Z",
     "iopub.status.busy": "2023-07-26T10:52:26.210826Z",
     "iopub.status.idle": "2023-07-26T10:52:26.458716Z",
     "shell.execute_reply": "2023-07-26T10:52:26.458144Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 42000 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "ax.hist(Tnull,\n",
    "        bins=100,\n",
    "        density=True,\n",
    "        facecolor='y',\n",
    "        label='Null')\n",
    "xval = np.linspace(-4.2, 4.2, 1001)\n",
    "ax.plot(xval,\n",
    "        t_dbn.pdf(xval, D_.shape[0]-2),\n",
    "        c='r')\n",
    "ax.axvline(observedT,\n",
    "           c='b',\n",
    "           label='Observed')\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"Null Distribution of Test Statistic\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bd21158",
   "metadata": {},
   "source": [
    "The re-sampling-based null distribution is almost identical to the theoretical null distribution, which is displayed in red.\n",
    "\n",
    "Finally, we implement the plug-in re-sampling FDR approach outlined in\n",
    "Algorithm 13.4. Depending on the speed of your\n",
    "computer, calculating the FDR for all 2,308 genes in the `Khan`\n",
    "dataset may take a while.  Hence, we will illustrate the approach on a\n",
    "random subset of 100 genes.  For each gene, we first compute the\n",
    "observed test statistic, and then produce 10,000 re-sampled test\n",
    "statistics. This may take a few minutes to run.  If you are in a rush,\n",
    "then you could set `B` equal to a smaller value (e.g. `B=500`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "3b7392cb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T10:52:26.461754Z",
     "iopub.status.busy": "2023-07-26T10:52:26.461349Z",
     "iopub.status.idle": "2023-07-26T12:21:43.650925Z",
     "shell.execute_reply": "2023-07-26T12:21:43.650584Z"
    }
   },
   "outputs": [],
   "source": [
    "m, B = 100, 10000\n",
    "idx = rng.choice(Khan['xtest'].columns, m, replace=False)\n",
    "T_vals = np.empty(m)\n",
    "Tnull_vals = np.empty((m, B))\n",
    "\n",
    "for j in range(m):\n",
    "    col = idx[j]\n",
    "    T_vals[j] = ttest_ind(D2[col],\n",
    "                          D4[col],\n",
    "                          equal_var=True).statistic\n",
    "    D_ = np.hstack([D2[col], D4[col]])\n",
    "    D_null = D_.copy()\n",
    "    for b in range(B):\n",
    "        rng.shuffle(D_null)\n",
    "        ttest_ = ttest_ind(D_null[:n_],\n",
    "                           D_null[n_:],\n",
    "                           equal_var=True)\n",
    "        Tnull_vals[j,b] = ttest_.statistic\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b92df1b",
   "metadata": {},
   "source": [
    "Next, we compute the number of rejected null hypotheses $R$, the\n",
    "estimated number of false positives $\\widehat{V}$, and the estimated\n",
    "FDR, for a range of threshold values $c$ in\n",
    "Algorithm 13.4. The threshold values are chosen\n",
    "using the absolute values of the test statistics from the 100 genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "cac15616",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T12:21:43.653105Z",
     "iopub.status.busy": "2023-07-26T12:21:43.652856Z",
     "iopub.status.idle": "2023-07-26T12:21:43.895651Z",
     "shell.execute_reply": "2023-07-26T12:21:43.895292Z"
    }
   },
   "outputs": [],
   "source": [
    "cutoffs = np.sort(np.abs(T_vals))\n",
    "FDRs, Rs, Vs = np.empty((3, m))\n",
    "for j in range(m):\n",
    "   R = np.sum(np.abs(T_vals) >= cutoffs[j])\n",
    "   V = np.sum(np.abs(Tnull_vals) >= cutoffs[j]) / B\n",
    "   Rs[j] = R\n",
    "   Vs[j] = V\n",
    "   FDRs[j] = V / R\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6779ea0",
   "metadata": {},
   "source": [
    "Now, for any given FDR, we can find the genes that will be\n",
    "rejected. For example, with FDR controlled at 0.1, we reject 15 of the\n",
    "100 null hypotheses. On average, we would expect about one or two of\n",
    "these genes (i.e. 10% of 15) to be false discoveries.  At an FDR of\n",
    "0.2, we can reject the null hypothesis for 28 genes, of which we\n",
    "expect around six to be false discoveries.\n",
    "\n",
    "The variable `idx` stores which\n",
    "genes were included in our 100 randomly-selected genes. Let’s look at\n",
    "the genes whose estimated FDR is less than 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9661eb10",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T12:21:43.897943Z",
     "iopub.status.busy": "2023-07-26T12:21:43.897797Z",
     "iopub.status.idle": "2023-07-26T12:21:43.901002Z",
     "shell.execute_reply": "2023-07-26T12:21:43.900637Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['G0097',\n",
       " 'G0129',\n",
       " 'G0182',\n",
       " 'G0714',\n",
       " 'G0812',\n",
       " 'G0941',\n",
       " 'G0982',\n",
       " 'G1020',\n",
       " 'G1022',\n",
       " 'G1090',\n",
       " 'G1320',\n",
       " 'G1634',\n",
       " 'G1697',\n",
       " 'G1853',\n",
       " 'G1854',\n",
       " 'G1994',\n",
       " 'G2017',\n",
       " 'G2115',\n",
       " 'G2193']"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(idx[np.abs(T_vals) >= cutoffs[FDRs < 0.1].min()])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "001e3fc1",
   "metadata": {},
   "source": [
    "At an FDR threshold of 0.2, more genes are selected, at the cost of having a higher expected\n",
    "proportion of false discoveries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "18ad4900",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T12:21:43.902778Z",
     "iopub.status.busy": "2023-07-26T12:21:43.902684Z",
     "iopub.status.idle": "2023-07-26T12:21:43.905535Z",
     "shell.execute_reply": "2023-07-26T12:21:43.905228Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['G0097',\n",
       " 'G0129',\n",
       " 'G0158',\n",
       " 'G0182',\n",
       " 'G0242',\n",
       " 'G0552',\n",
       " 'G0679',\n",
       " 'G0714',\n",
       " 'G0751',\n",
       " 'G0812',\n",
       " 'G0908',\n",
       " 'G0941',\n",
       " 'G0982',\n",
       " 'G1020',\n",
       " 'G1022',\n",
       " 'G1090',\n",
       " 'G1240',\n",
       " 'G1244',\n",
       " 'G1320',\n",
       " 'G1381',\n",
       " 'G1514',\n",
       " 'G1634',\n",
       " 'G1697',\n",
       " 'G1768',\n",
       " 'G1853',\n",
       " 'G1854',\n",
       " 'G1907',\n",
       " 'G1994',\n",
       " 'G2017',\n",
       " 'G2115',\n",
       " 'G2193']"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(idx[np.abs(T_vals) >= cutoffs[FDRs < 0.2].min()])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8767f70c",
   "metadata": {},
   "source": [
    "The next line  generates  Figure 13.11, which is similar\n",
    "to  Figure 13.9,\n",
    "except that it is based on only  a subset of the genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "28c276b6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T12:21:43.907278Z",
     "iopub.status.busy": "2023-07-26T12:21:43.907150Z",
     "iopub.status.idle": "2023-07-26T12:21:43.988689Z",
     "shell.execute_reply": "2023-07-26T12:21:43.988352Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 28364 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(Rs, FDRs, 'b', linewidth=3)\n",
    "ax.set_xlabel(\"Number of Rejections\")\n",
    "ax.set_ylabel(\"False Discovery Rate\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4b5d621",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
