Multiple Comparisons.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiple Comparisons Problem\n",
    "\n",
    "The multiple comparisons problem occurs when... making multiple comparisons. Pretty much any time you're working with multiple p-values.\n",
    "\n",
    "The more comparisons you make, the more likely you are to run into a false positive result where you see a p-value that is \"statistically significant\".\n",
    "\n",
    "Let's say you set a test with $alpha = 0.05$, meaning anytime your p-value is less than alpha, you call it statistically significant. However, instead of running the test on just 1 metric, you're looking at 2 metrics. What's the probability you would see atleast 1 false positive result even if none existed?\n",
    "\n",
    "To answer this question it's useful to understand the definition of a frequentist __p-value__: It's the probability of obtaining results atleast as extreme as the observed results __given that the null hypothesis is true__.\n",
    "\n",
    "Using this definition, its easy to see that if the null hypothesis were true and $alpha = 0.05$, then 5% of the time the data generating process from the null hypothesis would result in a statistically significant result by random chance alone.\n",
    "\n",
    "We can further quantify this false positive rate with a simple formula\n",
    "\n",
    "$$\n",
    "\\text{FPR} = (1-\\alpha)^n\n",
    "$$\n",
    "\n",
    "where $n$ is the number of metrics you are evaluating. The simulations below will confirm this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Does the multiple comparisons problem apply to the p values in a regression model?\n",
    "\n",
    "Of course it does. And its a more elegant way to illustrate the problem than running 20 t-tests. Let's simulate it out below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_metrics = 100\n",
    "N = 10000\n",
    "X = np.array([np.random.normal(0,1,size=N) for _ in range(n_metrics)])\n",
    "y = np.random.normal(0,1, size=N) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'p-value')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5wcZZ3v8c8vYUgGCIyQhEsSSPCQICeYAJMsJl7CTcJFCXJXXOOF6LKIupJjUFQuKtmFFVRYVlZRFj2SXVBOFoIsV0VQySSASCQaEM1MAgksgUgSmYTf+aNqQs9Mz0x3T1fXU1Xf9+s1r0x31XQ/1V15fvX8nkuZuyMiIsU1JO0CiIhIuhQIREQKToFARKTgFAhERApOgUBEpOB2SLsA1Ro5cqSPHz8+7WKIiGTKsmXLXnD3UeW2ZS4QjB8/nra2trSLISKSKWb2p762KTUkIlJwCgQiIgWnQCAiUnAKBCIiBadAICJScAoEIiIFl1ggMLMbzGydmf22j+1mZt80s1Vm9hszOzSpsoiISN+SbBF8H5jdz/bjgAPin3nAdQmWRURE+pDYhDJ3/7mZje9nl5OAf/fohgi/MrMWM9vb3dcmVSYpjtse7eCKu1ayZsNm9mlpZv6xk5hzyJi0iyUSpDRnFo8BVpc8bo+fUyDIiFAr29se7eDCHz/B5s5tAHRs2MyFP34CIIjySWVCPb/yKBOdxWY2z8zazKxt/fr1aRdHeKOy7diwGeeNyva2RzvSLhpX3LVyexDosrlzG1fctTKlEkm1Qj6/8ijNQNABjCt5PDZ+rhd3v97dW929ddSosmsmSYOFXNmu2bC5quclPCGfX3mUZiBYDPxtPHrocOBl9Q9kR8iV7T4tzVU9L+EJ+fzKoySHj/4I+CUwyczazeyjZvYJM/tEvMsS4BlgFfBvwLlJlUXqL+TKdv6xk2huGtrtueamocw/dlJKJZJqhXx+5VGSo4bOGmC7A3+f1PtLsuYfO6lbhyyEU9l2dSiqozG7Qj6/8ihz9yOQMIRe2c45ZEwwZZHqDfb80oij6lh0YZ4dra2trhvTlKeTX5KUlfOr5/BhiFoTl7/v4CDL2yhmtszdW8ttK3SLICsndiU0dl6SlKXzq78RR4Mpa57qi54yMY8gCXkbp6zhdpKkLJ1fSYw4ylt90VNhA0GWTuxK1HLy3/ZoBzMX3seEBXcwc+F9uTmppf6yNJwziRFHeasveipsIMjSiV2Jak/+vF/hSH1laThnEsOH81Zf9FTYQJClE7sS1Z78eb/CkfrK0tyMOYeM4fL3HcyYlmYMGNPSPOiO4rzVFz0VtrM4b+OUqx1ul/crHKmv0IcL91Tv4cN5qy96KmwgyNqJXYlqTv59WprpKFPp5+UKR+qvyHMz8lhflNI8goLSWGuRYtE8Aukl71c4IlI5BYICK3JTX0TeUNhRQyIiElGLQEQaIs9LNGSdAoGIJC5LaxUVkQKBiCQuqYXg8iatVpMCgYgkThMYB5Zmq0mBQKTOlAvvrZ4TGPP6+abZatKoIZE60mJ+5dVrraI8f75ptpoUCETqSIv5lVevheDy/PmmubCdUkMiPQwm9aBceN/qMYExj59v1/nWsWEzBpQu+tOohe0UCDIir3nR0Ay2w06L+SUrb59vz/PNYXswGNPA/+dKDWVAnvOioRls6iFL6/ZnUd4+33LnW1cQeGjBkQ272FOLIAM0BrtxBpt60GJ+yWr055t0SzyUVJcCQQaEcrIUQT1SD1rML1mN+nwbMa4/lFSXUkMZkPfb5IUkb6kHqV0tacLbHu1g5sL7mLDgDmYuvG/A9G0o55sCQQaEcrIUQRL3u5VsqrYlXktfXijnm1JDGaC8c2MptSNQfdqm1r68EM43BYKMCOFkEekSynDmJMtR7Q3rs9yXp0AgIlUJZUnppMtRbUs8lI7fWiR683ozmw18AxgKfMfdF/bYvi9wI9AS77PA3Zf095q6eb1IumYuvK9shdc19j3tcnSVpdGtlJ6BCaIWRCh9TKncvN7MhgLXAscA7cBSM1vs7itKdrsI+A93v87MDgKWAOOTKpOIDGygdEsoKZD+3i+NVkqW+/KSTA1NB1a5+zMAZnYzcBJQGggc2DX+fTdgTYLlEambUHLk9VZJuiWUFEhf5ehSOtSzUd9VVvvykhw+OgZYXfK4PX6u1MXA2WbWTtQa+GS5FzKzeWbWZmZt69evT6KsIhXL85IflYydD2U4c7ly9NT13eTxu6qntOcRnAV8393HAscDN5lZrzK5+/Xu3ururaNGjWp4IUVK5Xkp5ErSPqGMfS8tR1+GmuX2u6qnJFNDHcC4ksdj4+dKfRSYDeDuvzSz4cBIYF2C5RIZlFBy5EmoNO0TSgqkqxx9ddT2DAJd8vBd1VOSLYKlwAFmNsHMdgTOBBb32OfPwFEAZvYWYDig3I8EqWv5gL7G2WVhmOBAQkn7VKuvVkpfrYU8fFf1lFiLwN23mtl5wF1EQ0NvcPcnzexSoM3dFwOfBf7NzD5D1HE815MczypSo3JXnKWyUFlWIssjX/pqpVQzKayoEp1HkATNI8i3UEfjhDZmXSoX6jnVaKnMIxCpVigzVsvpK6ds0NBJVFK9UPozQpb2qCGR7UIejZP3pcCrXT5Z8kWBQIIR8micrHaiViLP8yKkMgoEEoyQr7pDGTufhJBbYtIY6iOQuqhHh1y1y/42Wl5zzSG3xKQxFAgKIOlRE/Xq5M3y0MUsC2XtIEmPAkHONWIkTq13Zionr1fdIQu9JSbJUx9BzjUi/6vUQrbluf9DKqMWQQZVk+ppRCWt1EL2qSVWbGoRZEy1Q/0aMRInz0MrRYpAgSBjqk31NKKSVmpBJNuUGqpSaVpmt+YmzGDDps6GjXCpNtXTqJE4Si2IZJcCQRV6jsDZsLlz+7ZGrYtTSz5elbSkQYu9ZYdSQ1Uol5Yp1YjZmMrHSxZo2YpsUSCoQiUjbZIeMql8vGSBlq3IFqWGqtBXWqbnPklTqqf+QkljhFKOwcra3JK8fO61UougCuXSMqWUosmmUNIYoZSjHkJeQLCnPH3utVIgqELPtExLcxNv2qlJKZqMCyWNEUo56iFLfVl5+txrpdRQlZSWyZ9Q0hiDLUdI6Y0sLSAYyvefJgUCKbxQlsgYTDlCvM1nVi6a+vrcnehe1aEGsHpSaigAuk1gukJJYwymHEpv1K6/vr+i9BcoEKRMHVXpC2VI7mDKofRG7Uo/93KKEFCVGkpZPdfyl9qFksaopBzl+gJCSW9lVdfnPmHBHXiZ7XkPqAoECQhtmegQhdSxmSV99QWcctgYbl3WoZvLDFJRA6pSQ3UW4jLRjTZQn4fSYbXrqwV5/1Prg0hvZV0o/UWNphZBnVWb6snbbQIrGb2idFjt+mtBhpLeyrJQh70m3YJWIKizUJeJbpRKKvmipsPqoaipi0YKLaA2YmiwAkGdFX2Z6EoqeVVmtctbC1IG1ogWtPoI6qyoOcYulfR5FP0zGoxQhrpK4zSiBa0WQZ3lLdVTrUquWIv+GQ1WnlqQMrBGtKDNvdyo2Tq9uNls4BvAUOA77r6wzD6nAxcTzeh+3N3f399rtra2eltbWwKllXrR0FCR+unZRwDRxVW1LUEzW+bureW2JdYiMLOhwLXAMUA7sNTMFrv7ipJ9DgAuBGa6+0tmNjqp8kjj6IpVpH4a0YJOMjU0HVjl7s8AmNnNwEnAipJ9zgGudfeXANx9XYLlERHJpKQvrpLsLB4DrC553B4/V2oiMNHMHjKzX8WppF7MbJ6ZtZlZ2/r16xMqrohIMaU9amgH4ABgFnAW8G9m1tJzJ3e/3t1b3b111KhRDS6iiEi+JRkIOoBxJY/Hxs+VagcWu3unu/8R+D1RYBARkQZJso9gKXCAmU0gCgBnAj1HBN1G1BL4npmNJEoVPZNgmQpPI3pEpKfEAoG7bzWz84C7iIaP3uDuT5rZpUCbuy+Ot73bzFYA24D57v5iUmUquhDvYiUi6atoHoGZ7QR8FtjX3c+Jh31Ocvfbky5gT5pHULuZC+8rOzFlTEszDy04MoUSiUij9DePoNI+gu8BfwXeFj/uAL5Sh7JJA2mxNxEpp9LU0Jvd/QwzOwvA3TeZmSVYLkmAFnsrJvULyUAqbRG8ZmbNRMtAYGZvJmohSIZosbfi0U2ApBKVtgi+DPwUGGdmPwRmAnOTKpQkI++LvenKtzfdBEgqUVEgcPe7zWw5cDhgwKfc/YVESyaJyOs6QBoRVZ76haQSFQUCM3tn/OvG+N+DzAx3/3kyxZJQZOUqW1e+5YXaL9TZ2Ul7eztbtmxJtRx5NHz4cMaOHUtTU1PFf1Npamh+6fsQLSi3DNCYwxzL0lW2rnzLC/WOZu3t7YwYMYLx48ejcSf14+68+OKLtLe3M2HChIr/rqLOYnd/T8nPMcBk4KUayyoZ0d9VdmgquTNaEYV6R7MtW7awxx57KAjUmZmxxx57VN3SqnVmcTvwlhr/VjIiS1fZoV75hiDUfiEFgWTU8rlW2kfwLeKho0StiKnA8qrfTTIl1PxyOXkfESWSpErnEbQR9QksA34JfM7dz06sVBKErM07mHPIGB5acCR/XHgCDy04UkFAqvKxj32MFStWDLxjGWvWrOHUU0/d/viss87irW99K1dddRVf+tKXuOeee+pVTACuvvpqNm3aVLfXS/SexUnQWkONlZVRQ5Itv/vd73jLW/KZXX7uued4+9vfzqpVqxJ7j/Hjx9PW1sbIkSPLbi/3+da81pCZPWFmvynz84SZ/abWg5Ds0FW25NGrr77KCSecwJQpU5g8eTKLFi1i1qxZdF1kfve732XixIlMnz6dc845h/POOw+AuXPncv755zNjxgz2339/brnlFgCeffZZJk+eDMC73/1uOjo6mDp1Kg8++CBz587dvt/SpUuZMWMGU6ZMYfr06WzcuJFnn32Wd7zjHRx66KEceuihPPzwwwA88MADzJo1i1NPPZUDDzyQD3zgA7g73/zmN1mzZg1HHHEERxxxRF0+j4H6CE6sy7uIiPRn1qzez51+Opx7LmzaBMcf33v73LnRzwsvQElaBoAHHuj37X7605+yzz77cMcddwDw8ssvc9111wFRmueyyy5j+fLljBgxgiOPPJIpU6Zs/9u1a9fyi1/8gqeeeor3vve93VJCAIsXL+bEE0/kscceA6KgAvDaa69xxhlnsGjRIqZNm8Yrr7xCc3Mzo0eP5u6772b48OH84Q9/4KyzztoekB599FGefPJJ9tlnH2bOnMlDDz3E+eefz9e//nXuv//+PlsE1eq3ReDuf+rvpy4lEBFpsIMPPpi7776bz33uczz44IPstttu27c98sgjvOtd72L33XenqamJ0047rdvfzpkzhyFDhnDQQQfx/PPPV/yeK1euZO+992batGkA7Lrrruywww50dnZyzjnncPDBB3Paaad166eYPn06Y8eOZciQIUydOpVnn312cAfeh0pHDR0OfItoyOiORDeaedXdd02kVCIyKJnr2+nvCn6nnfrfPnLkgC2AniZOnMjy5ctZsmQJF110EUcddVTFfzts2LDtv9ejj/Wqq65izz335PHHH+f1119n+PDhZd9r6NChbN26ddDvV06lo4auIbql5B+AZuBjwLWJlEikQW57tIOZC+9jwoI7mLnwvtysyKkVRwe2Zs0adtppJ84++2zmz5/P8uVvjIafNm0aP/vZz3jppZfYunUrt956a13ec9KkSaxdu5alS5cCsHHjRrZu3crLL7/M3nvvzZAhQ7jpppvYtm3bAK8EI0aMYOPGjQPuV6mKb17v7quAoe6+zd2/B8yuWylEGizPlWWWZoSn5YknnmD69OlMnTqVSy65hIsuumj7tjFjxvD5z3+e6dOnM3PmTMaPH98tdVSrHXfckUWLFvHJT36SKVOmcMwxx7BlyxbOPfdcbrzxRqZMmcJTTz3FzjvvPOBrzZs3j9mzZ9ets7jSW1X+HDga+A7wHLAWmOvuU/r9wwRo+Gi2hJqiyPNtOycsuINy/6sN+OPCExpdnLJCHz76l7/8hV122YWtW7dy8skn85GPfISTTz457WJVrK7DR0t8MN73POBVYBxwyiDKKQUQ8lV32stnJJmW0rpLg3fxxRczdepUJk+ezIQJE5gzZ07aRUpUpWsNHQbc4e6vAJckWB7JkZCXhk5z+YykV3XVukuDd+WVV6ZdhIaqtEXwHuD3ZnaTmZ1oZrUuVicFkvZVd3/SXD4j6Rx+qCuO9pS1VQ2yopbPtdI7lH3YzJqA44hGD11rZne7+8eqfkcpjJAXrUtzkbpGBMhQVxztMnz4cF588UUtRV1nXfcjKB2CWomKr+zdvdPM7iRahbQZmEM0jFSkrNBTFGlVliEHyMGqdHDA2LFjaW9vZ/369SmUMt+67lBWjUonlB0HnAHMAh4gGj10enXFk6LR0tDlhRIg6z2iq5q+j6ampqruoCXJqnT46I+ARcCd7v7XxEvVj0YMHw11yKPkR9rnWM9KG6JgNJi+hDwPyU1Dvc+R/oaPVtpHcFbJi53o7rfXXJrAZek+vZJdaefwkxjRFfLggKxpdD1U8cziEpfWvRQB0azMcOR1CYgQJFFpa/5C/TS6HqolEOS6i19XNWEIeTJaHiRRaWftjnYha3Q9VFEgMLPhZvYPZvZj4CUz+4yZDTg+ycxmm9lKM1tlZgv62e8UM3MzK5u/aiRd1YRBLbNkJVFpZ2X+QhY0uh6qdPjovwMbiZaiBng/cBNwWl9/YGZDiVYoPQZoB5aa2WJ3X9FjvxHAp4BfV1f0ZIQyoqPo6nlFlHbHbIiSGtGVdt9HXjS6Hqo0EEx294NKHt9vZgPd5Xk6sMrdnwEws5uBk4Cef3cZ8I/A/ArLkigNeQxDvcbaq/O/b6q0w9XoeqjSQLDczA53918BmNnfAAON4RwDrC553A78TekOZnYoMM7d7zCzVAOBrhrDUq8ropDXOxLpTyMDdTWLzj1sZn+OH+8LrDSzJwB397dW+8ZmNgT4OjC3gn3nAfMA9t1332rfakC6agxPva6I1PkvMrBKA0EtN6HpIFquusvY+LkuI4DJwAPxWiN7AYvN7L3u3q214e7XA9dDNKGshrL0S1eNYarHFVGel3MQqZdKJ5TVcqP6pcABZjaBKACcSdTJ3PWaLwMjux6b2QPABT2DQCPoqrEx0ki/qfNfZGCJLSft7lvN7DzgLqKb3d/g7k+a2aVAm7svTuq9qxXqVWOe+i3SSr+p819kYBWtNRSSJNYaSmLdlTyWaTC0Do1IJK0LvEGvNZR3IV415q3fQuk3CU0aFXKoA1MUCGKhjanOW8UZavpNiimtCjnUC7xa1hqSBsjbUhdah0ZCktYSJqFe4CkQBCpvFafWoZGQpFUhh3qBp9RQwmrNQ4bYbzFYoaXfpLjSSlWGOpy5EIEgrV76weYhVXGKJCOtCjnUC7zcB4I0e+lD7RiS+srTfI+iSLNCDvECL/eBIM3KONSOodBkuSINdTigDCzECjktue8sTrMybkTHUNZv55j1O5HpBjqSB7kPBGn20ic98ifrlShkvyJVq0/yIPeBIM1hmEkPmcx6JQrhV6QDtbhCHQ4ojZP1VjkUoI8g7V76JPOQoVeilQh5xnEl+f9QhwNKY+Sljyj3gQDy2ykUciVaqZAr0koGGqR9oSHpysvIwEIEgrwKuRKtVMgVaaUtrrxeaEC2R3Q1Qh5a5aBAkGkhV6LVGExFmmRFlYcW12DkJe2RpLycIwoEGZfnq9GBJF1R5aHFNRh5SXskKS/niAKBZFbSFVVeWly1SjLtkZeUU3/nSOkx7tbchBls2NQZ5PEqEEhmNSI/W+QWV1Jpj7ylnMqdIz2PccPmzu3bQjze3M8jkPzSGP5kJTUHJw/zXwZS7hhLhXa8CgSSWXm7Z0NokpoQmYeRNgNNIqvkWEI6XqWGJLOKnsNvhCRSY1kfaVNJaquvYywV0vEqEEimFTmHn1VZH2lTySCFcsdYKrTjVSAQkYbKekuuktRWz2PUqCERkR6y3JKrNLWVpWNUZ7EA+VhBUaQR8jhIwdw97TJUpbW11dva2hr6nnmZ/NKXnp1fAAY40UiRvB2vyGBlsU4ws2Xu3lpum1JDA8jb5JdyynV+dV0e5PF4RQYrS2mfSig1NIAiTH4ZaDxz3o5XRLpTIBhAHia/DKSS8cx5Ol4R6U6BYABFWMagXOdXT3k6XhHpLtFAYGazzWylma0yswVltv+Dma0ws9+Y2b1mtl+S5alFHkcI9FS6lABEHcWl8na8ItJdYp3FZjYUuBY4BmgHlprZYndfUbLbo0Cru28ys78D/gk4I6ky1SLrk18qVdr5lcURESJSu8SGj5rZ24CL3f3Y+PGFAO5+eR/7HwJc4+4z+3vdNIaPiohkXX/DR5NMDY0BVpc8bo+f68tHgTvLbTCzeWbWZmZt69evr2MRRUQkiM5iMzsbaAWuKLfd3a9391Z3bx01alRjCyciknNJTijrAMaVPB4bP9eNmR0NfAF4l7v/NcHyiIhIGUkGgqXAAWY2gSgAnAm8v3SHuF/g28Bsd1+XYFlEJMM0gCFZiQUCd99qZucBdwFDgRvc/UkzuxRoc/fFRKmgXYD/NDOAP7v7e5Mqk4hkTxGWeUlbomsNufsSYEmP575U8vvRSb5/yHSFI1KZSm4EI4OjRedSoCucYlLwr00RlnlJWxCjhoqmCAvZSXddwb9jw2acN4K/7vswsHou86L7bpSnQJACXeEUj4J/7eq1zIuCcd8UCFJQhIXspDsF/9qVroVlRDdLuvx9B1edVlMw7pv6CFIw/9hJve4IpoXd8qmrX6CvhVwU/CtTjxvBKBj3TYEgBUVZyK7oyt0CtJSCf2NVetP5IlIgSEnebnUnvZVLRXTRvaAbTy3xvikQiCSkr5SDAQ8tOLKxhRG1xPuhQCCSEKUiwqOWeHkaNSSSkCLc3U7yQS0CkYSEmIoIfXZz6OXLKwUCkQSFlIoIdWmTrsq/Y8NmDLYPtQ2lfEWg1JBIQYQ4oap0ti/Qa75F2uUrCrUIRDKgHimTECdU9TfEtosmfCVPLQKRwNVrjZwQlzappJLXKKvkKRDklFZZzI96pXRCHMU0UCWfdvmKQqmhBmrUiIhQOwWlNvVK6YQ4iqncbN+uDmPNvm4cBYIGaWTlrDs65Us9J6aFNIoJwgxORaRA0CCNrJxD7BSU2uV9jZzQglMRFS4QpDVhpZGVc96XNijapCNdNUvSChUI0sydN7JyzvMVZFH7P3TVLEkq1KihNCfUNHLERr3u6BSiECdFiWRd9loEK1fCrFndnzv9dDj3XNi0CY4/vvffzJ0Lc+eyec1z3Hzb5b02/+CQ44EjYfVq+OAHe//9Zz8L73lP9N4f/3jv7RddxNJfP8WEhV9ij7+8ROfQJraM249d94sr3q99jTkzZrD7423scumXeW3rNnbcYSj77t7MyAeHwdVXw9SpcM898JWv9H79b38bJk2C//ov+Od/7r39pptg3DhYtAiuuw6AOfEPALfcAiNHwve/H/30tGQJ7LQT/Mu/wH/8R+/tDzwQ/XvllXD77d23NTfDnXdGv192Gdx7b/fte+wBt94a/X7hhfDLX3bfPnYs/OAH0e+f/jQ89lj37RMnwvXXs/Sr1/C/L53PL17bwmtDm/jzbnvy4s4trBi9P5cePS9KsZ19NrS3d//7t70NLo+/81NOgRdf7L79qKPgi1+Mfj/uONjco9V24olwwQXR7z3PO6j43OOFF+DUU3tv/7u/gzPOGNS5x9FHR5/bpz/de/vXvgYzZsDDD8PnP997ewLnXjc5OPcAmDcPfv/77tunTo0+P8jmuVcie4FgEPbarXwaZveddxzU66784W1M/uF3aO78KwA7butk6J+e5hV4IxgA75w4CvZtGdR7FdHSr17D5Esu2P75DtvWyf7/031exD4tzdBe7q+7e+VPHQxf/SeatnXSObSJdTvswtgvJlHqyi396jXs908XM+qVF3tfRIg0grtn6uewww7zWv1kebsfeNGdvt/nbt/+c+BFd/pPlrfX/Jru7mtbRrtDr5+1LaMH9boS6evzXb3rqKq+w0e+8i3f1DSs22tsahrmj3zlWw04iurL9MhXvuVrW0b7NszXtowuW85K9hFxdwfavI961aLt2dHa2uptbW01/30SI05etyEMKXN78tcxhvjrg3pt6f/zfcfl91T8HT73pj3Za8O63s+3jGavl56vS1mr1VeZXmoewfCtr21vBQFsbhrGb798JdO+cB7Qu6VUbh+RLma2zN1by24rWiBIQj0rmKVfvYZxV17G6A3rWdcyitUXfLHw/6nr9fmGGLD7KpMTzbDtqfSYQwxsEq7+AkGhRg0lZfUFX2Rz07Buz21uGsbqC6pLPndd4e21YR1DcPbasI7Jl1zA0q9eU8/iZk69Pt91LaOqer4Rqn3v0RvWl/29r31EKqFAUAfTvnAev/3ylTzXMprXMZ5rGV1T83zclZd1a+YDNHf+lXFXXlbP4mZOvT7fegWUeuqrTBuady27f2ngCDGwSUb11XkQ6s9gOotDtw0r2ym6DUu7aLkRYudquTJV0rEdSud3iJ+p9EY/ncWJVtrAbGAlsApYUGb7MGBRvP3XwPiBXjPPgSD00Uf6D99YWRg1FEowkoGlEgiAocDTwP7AjsDjwEE99jkX+Nf49zOBRQO9bp4DQcj/qUIum6Qn9IsXeUN/gSDJPoLpwCp3f8bdXwNuBk7qsc9JwI3x77cAR5lZucEShVCvXHgS1H8h5ajDOh+SDARjgNUlj9vj58ru4+5bgZeBPXq+kJnNM7M2M2tbvz7fJ9i0L5zHXi89zxB/nb1eej6IIAD6Dy/lqcM6HzIxasjdr3f3VndvHTVKJ1ga9B9eyglxJJZUL8lA0AGMK3k8Nn6u7D5mtgOwG9BjZSYJgf7DSzkhpzOlconNLI4r9t8DRxFV+EuB97v7kyX7/D1wsLt/wszOBN7n7qf397ohziwuCs16Fsmu1JaYMLPjgauJRhDd4O5fNbNLiXqvF5vZcOAm4BDgf4Az3f2Z/l5TgUBEpHr9BYJEl8tsBEwAAAcjSURBVKF29yXAkh7Pfank9y3AaUmWQURE+peJzmIREUmOAoGISMEpEIiIFJwCgYhIwSkQiIgUnAKBiEjBKRCIiBRc5u5ZbGbrgT/V+OcjgRfqWJws0DEXg465GAZzzPu5e9nFwTIXCAbDzNr6mlmXVzrmYtAxF0NSx6zUkIhIwSkQiIgUXNECwfVpFyAFOuZi0DEXQyLHXKg+AhER6a1oLQIREelBgUBEpOAKEwjMbLaZrTSzVWa2IO3yJMHMxpnZ/Wa2wsyeNLNPxc/vbmZ3m9kf4n/flHZZ68nMhprZo2Z2e/x4gpn9Ov6uF5nZjmmXsZ7MrMXMbjGzp8zsd2b2tgJ8x5+Jz+nfmtmPzGx43r5nM7vBzNaZ2W9Lniv7vVrkm/Gx/8bMDh3MexciEJjZUOBa4DjgIOAsMzso3VIlYivwWXc/CDgc+Pv4OBcA97r7AcC98eM8+RTwu5LH/whc5e7/C3gJ+GgqpUrON4CfuvuBwBSiY8/td2xmY4DzgVZ3n0x0x8Mzyd/3/H1gdo/n+vpejwMOiH/mAdcN5o0LEQiA6cAqd3/G3V8DbgZOSrlMdefua919efz7RqIKYgzRsd4Y73YjMCedEtafmY0FTgC+Ez824EjglniXvB3vbsA7ge8CuPtr7r6BHH/HsR2A5vhe6DsBa8nZ9+zuPye6ZW+pvr7Xk4B/98ivgBYz27vW9y5KIBgDrC553B4/l1tmNp7oXtC/BvZ097XxpueAPVMqVhKuBv4P8Hr8eA9gg7tvjR/n7bueAKwHvhenw75jZjuT4+/Y3TuAK4E/EwWAl4Fl5Pt77tLX91rXOq0ogaBQzGwX4Fbg0+7+Suk2j8YL52LMsJmdCKxz92Vpl6WBdgAOBa5z90OAV+mRBsrTdwwQ58VPIgqC+wA70zuFkntJfq9FCQQdwLiSx2Pj53LHzJqIgsAP3f3H8dPPdzUb43/XpVW+OpsJvNfMniVK9x1JlD9viVMIkL/vuh1od/dfx49vIQoMef2OAY4G/uju6929E/gx0Xef5++5S1/fa13rtKIEgqXAAfEogx2JOpoWp1ymuovz498FfufuXy/ZtBj4UPz7h4D/1+iyJcHdL3T3se4+nug7vc/dPwDcD5wa75ab4wVw9+eA1WY2KX7qKGAFOf2OY38GDjezneJzvOuYc/s9l+jre10M/G08euhw4OWSFFL13L0QP8DxwO+Bp4EvpF2ehI7x7URNx98Aj8U/xxPlze8F/gDcA+yedlkTOPZZwO3x7/sDjwCrgP8EhqVdvjof61SgLf6ebwPelPfvGLgEeAr4LXATMCxv3zPwI6I+kE6ilt9H+/peASMaCfk08ATRiKqa31tLTIiIFFxRUkMiItIHBQIRkYJTIBARKTgFAhGRglMgEBEpOAUCyS0z28vMbjazp81smZktMbOJCb/nwwNsbzGzcwfx+h83sw+b2VQz+3atryNSSoFAcimeePQT4AF3f7O7HwZcSEJr8HTNcHX3GQPs2gLUHAiAdwA/B94V/ysyaAoEkldHAJ3u/q9dT7j74+7+YDwb84p4bfsnzOwMgLj1cELX/mb2fTM71czGm9mDZrY8/pkRb58VP7+YaKYrZvaX+N9dzOzeeP8nzKxrtduFwJvN7DEzuyLed76ZLY3Xlb+k3MHE6/E/BpxMtITIJcAXzOxfy+0vUg1NKJNcMrPzgQnu/pky204BPkG0cNlIoiVI/oboHg5z3P1D8VIkTwMTiWZxvu7uW8zsAOBH7t5qZrOAO4DJ7v7H+LX/4u67dC2X7O6vmNlI4FdEa8fvRzQDenK8/7uJlkn4ePw+i4F/8mhJ4p7lHgXc4O7vMbNH3H16nT4uKbgdBt5FJHfeTlSZbyNa1OtnwDTgTuAbZjaMKEj83N03x/cAuMbMpgLbiIJDl0e6gkAPBnzNzN5JtET2GMqnpd4d/zwaP96FKGCUS/scCjxuZrsCG6o6YpF+KBBIXj3JGwuSVSS+4n8AOBY4g2hFU4DPAM8T3Q1sCLCl5M9e7ePlPgCMAg5z9854hdThZfYz4HJ377Pj18xGA/8NjI7f+0xgRJwqOsXdn67oAEX6oD4Cyav7gGFmNq/rCTN7q5m9A3gQOMOiex2PIrrj1yPxbouADxN1yv40fm43YK27vw58kOhWiQPZjeheCZ1mdgRRSghgIzCiZL+7gI/E95DAzMbEFf927r7O3acCy4nutvcD4MPuPlVBQOpBgUByyaPOr5OBo+Pho08ClxPd5eknRCt3Pk4UMP6PR8s7Q3Tl/S7gHo9uawrwL8CHzOxx4ED6bgWU+iHQamZPAH9LtHIm7v4i8FDcUX2Fu/838H+BX8b73kL3QAFsv+/2Hu7+AjAD+EV1n4hI39RZLCJScGoRiIgUnAKBiEjBKRCIiBScAoGISMEpEIiIFJwCgYhIwSkQiIgU3P8HIQhD9AqDYDwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = sm.OLS(y, sm.add_constant(X.T)).fit()\n",
    "plt.scatter(range(n_metrics+1), model.pvalues )\n",
    "sigs = model.pvalues[model.pvalues < 0.05]\n",
    "sig_idxs = np.where(model.pvalues < 0.05)[0]\n",
    "plt.scatter(sig_idxs, sigs,color=\"r\" )\n",
    "plt.axhline(0.05, color=\"r\", ls=\"--\", label=\"significant\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Covariate #\")\n",
    "plt.ylabel(\"p-value\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How likely are we see a false positive signal given we look at 20 metrics?\n",
    "To do this we'll fit a regression on 20 metrics without an effect and test if we see 1 or more pvalues less than alpha = 0.05. We'll repeat this simulation 1000 times to estimate the probability it could occur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "n_metrics = 20\n",
    "N = 10000\n",
    "\n",
    "res = []\n",
    "for _ in range(1000):\n",
    "    X = np.array([np.random.normal(0,1,size=N) for _ in range(n_metrics)])\n",
    "    y = np.random.normal(0,1, size=N) \n",
    "\n",
    "    model = sm.OLS(y, sm.add_constant(X.T)).fit()\n",
    "    pvals = model.pvalues[1:] # ignore the coefficient on the intercept\n",
    "    res.append( (pvals < 0.05).sum() >= 1 )  \n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.642"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6415140775914581"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - 0.95**(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How about 2 metrics?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09750000000000003"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 - 0.95**(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.107"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "n_metrics = 2\n",
    "N = 10000\n",
    "\n",
    "res = []\n",
    "for _ in range(1000):\n",
    "    X = np.array([np.random.normal(0,1,size=N) for _ in range(n_metrics)])\n",
    "    y = np.random.normal(0,1, size=N) \n",
    "\n",
    "    model = sm.OLS(y, sm.add_constant(X.T)).fit()\n",
    "    pvals = model.pvalues[1:] # ignore the coefficient on the intercept\n",
    "    res.append( (pvals < 0.05).sum() >= 1 )  \n",
    "\n",
    "np.mean(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Can we fix this?\n",
    "\n",
    "A common method to control the false positive rate is to use the Benjamini-Hochberg procedure. You may have heard of a bonferroni correction, which is more popular, however it tends to reduce the power to detect an effect far more than the Benjamini-Hochberg procedure.\n",
    "\n",
    "The code below shows how the Benjamini-Hochberg procedure controls the false positive rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.634, 0.051)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_metrics = 20\n",
    "N = 10000\n",
    "\n",
    "res = []\n",
    "adj_res = []\n",
    "for _ in range(1000):\n",
    "    X = np.array([np.random.normal(0,1,size=N) for _ in range(n_metrics)])\n",
    "    y = np.random.normal(0,1, size=N) \n",
    "\n",
    "    model = sm.OLS(y, sm.add_constant(X.T)).fit()\n",
    "    \n",
    "    pvals = model.pvalues[1:] # ignore the coefficient on the intercept\n",
    "    _,adj_pvals,_,_ = sm.stats.multipletests(pvals, method=\"fdr_bh\")\n",
    "    \n",
    "    res.append( (pvals < 0.05).sum() >= 1 )  \n",
    "    adj_res.append( (adj_pvals < 0.05).sum() >= 1 )  \n",
    "\n",
    "np.mean(res), np.mean(adj_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Amazing, we reduced the false positive rate from 0.634 back down to 0.05. It's important to note that this comes at a price - the price you're paying is a reduced power to detect a true effect that may exist. But it does a fantastic job at controlling false positives\n",
    "\n",
    "[TODO: add more details on benjamini hochberg correction later]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What happens when we add in a few real effects? Can the Benjamini Hochberg procedure detect them? And how does it compare to the Bonferroni correction "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}