goodness_of_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Goodness of Fit Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook explores a few techniques for checking whether a statistical model is \"adequate,\" i.e., whether it fits the data so poorly that one can reject the hypothesis that the model generated the data.\n",
    "\n",
    "We will look at two kinds of models: models with no fitted parameters and models with free parameters fitted to the data. One should expect that a model that is adjusted to fit the data by picking the values of free parameters will fit better than a model in which the values of those parameters are pre-specified without looking at the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's very common to fit models to data (e.g., normal, log-normal, linear regression, logistic regression, Poisson regression) without worrying about where the model comes from or even whether it is plausible.\n",
    "\n",
    "In general, if you test whether a model fits, and then report the fitted model, $P$-values, confidence intervals, etc., only if it does fit, inferences will be \"selective inferences.\" \n",
    "The same thing goes for procedures that select which variables to include in a model, select transformations of the variables, select the functional form of the model, etc.: The $P$-values, confidence intervals, and so on will be misleading if you do not adjust appropriately for the selection process. \n",
    "\n",
    "For examples, see: \n",
    "\n",
    "+ Benjamini, Y. and D. Yekutieli, 2005. False Discovery Rate-Adjusted Multiple Confidence Intervals for Selected Parameters, _Journal of the American Statistical Association, Theory\n",
    "and Methods_, _100_(469), DOI 10.1198/016214504000001907\n",
    "\n",
    "+ Berk, R., L. Brown, and L. Zhao, 2009. Statistical Inference After Model Selection,\n",
    "_J. Quantitative Criminology_, DOI: 10.1007/s10940-009-9077-7 \n",
    "\n",
    "+ Berk, R., L. Brown,  A. Buja, K. Zhang, and L. Zhao, 2013.\n",
    "Valid Post-Selection Inference, _Annals of Statistics, 41_, 802-837\n",
    "\n",
    "+ Olshen, R., 1974. The Conditional Level of the $F$ test, _Journal of the American\n",
    "Statistical Association, Theory and Methods, 68_(343), 692-698.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models with no Free Parameters "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kolmogorov-Smirnov Test\n",
    "\n",
    "One particularly useful test for models for real-valued data is the Kolmogorov-Smirnov test, one of many goodness-of-fit tests based on the cumulative distribution function.\n",
    "\n",
    "Suppose $\\{X_j \\}_{j=1}^n$ are IID with a continuous distribution.\n",
    "Define $1_A$ = {1, if $A$; 0, otherwise}.\n",
    "The theoretical cdf of $X_i$ is\n",
    "\\begin{equation*} \n",
    "\tF(x) \\equiv \\mathbb P \\{ X_i \\le x \\} .\n",
    "\\end{equation*}\n",
    "\n",
    "Define the empirical cumulative distribution function\n",
    "\\begin{equation*} \n",
    "    \\hat{F}_n (x) \\equiv \\frac{1}{n} \\sum_{i=1}^n 1_{x \\ge X_i}\n",
    "\\end{equation*}\n",
    "\n",
    "Consider the one-sided Kolmogorov-Smirnov statistics\n",
    "\\begin{equation*} \n",
    "\tD_n^+ \\equiv \\sup_x (\\hat{F}_n(x) - F(x) )\n",
    "\\end{equation*}\n",
    "and\n",
    "\\begin{equation*} \n",
    "    D_n^- \\equiv \\sup_x (F(x) - \\hat{F}_n(x)),\n",
    "\\end{equation*}\n",
    "and the two-sided Kolmogorov-Smirnov statistic\n",
    "\\begin{equation*} \n",
    "\tD_n \\equiv \\sup_x |\\hat{F}_n(x) - F(x)|.\n",
    "\\end{equation*}\n",
    "Kolmogorov (1933) and Smirnov (1944) seem to have been the first to study these\n",
    "statistics, which have the same distribution&mdash;a distribution that does not depend on $F$ if $F$ is continuous.\n",
    "\n",
    "Below is a simulation: we take a random sample of size $n$ from a Uniform distribution and plot the ecdf and the true cdf. The ecdf is the stair-step function. In places it is above the true cdf, and in places it is below. $D^-$ measures how far $F$ ever gets above $\\hat{F}_n$. Note that as $n$ increases, the maximum distance between $F$ and $\\hat{F}_n$ tends to be smaller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division, print_function\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "from scipy.stats import binom\n",
    "import pandas as pd\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4b0d78c9e6c24a06a248d793173cd5d4"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plotUniformKolmogorov(n)>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def ecdf(x):\n",
    "    '''\n",
    "       calculates the empirical cdf of data x\n",
    "       returns the unique values of x in ascending order and the cumulative probabity at those values\n",
    "       NOTE: This is not an efficient algorithm: it is O(n^2), where n is the length of x. \n",
    "       A better algorithm would rely on the Collections package or something similar and could work\n",
    "       in O(n log n)\n",
    "    '''\n",
    "    theVals = sorted(np.unique(x))\n",
    "    theProbs = np.array([sum(x <= v) for v in theVals])/float(len(x))\n",
    "    if (theVals[0] > 0.0):\n",
    "        theVals = np.append(0., theVals)\n",
    "        theProbs = np.append(0., theProbs)\n",
    "    return theVals, theProbs\n",
    "\n",
    "\n",
    "def plotUniformKolmogorov(n):\n",
    "    sam = np.random.uniform(size=n)\n",
    "    v, pr = ecdf(sam)\n",
    "    fig, ax = plt.subplots(nrows=1, ncols=1)\n",
    "    ax.plot([0, 1], [0, 1], linestyle='--', color='r', label='true cdf')\n",
    "    ax.step(v, pr, color='b', where='post', label='ecdf')\n",
    "    ax.legend(loc='best')\n",
    "    dLoc = np.argmax(v-pr)\n",
    "    dMinus = (v-pr)[dLoc]\n",
    "    ax.axvline(x=v[dLoc], ymin=pr[dLoc], ymax=v[dLoc], color='g', linewidth='3')\n",
    "    ax.text(0.5, 0.1, r'$D_n^-=$' + str(round(dMinus, 3)), color='g', weight='heavy')\n",
    "    plt.show()\n",
    "\n",
    "interact(plotUniformKolmogorov, n=widgets.IntSlider(min=3, max=300, step=1, value=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Massart-Dvoretzky-Kiefer-Wolfowitz (MDKW) Inequality \n",
    "\n",
    "Dvoretsky, Kiefer, and Wolfowitz (1956, DKW) showed that\n",
    "\\begin{equation*}  \n",
    "\t\\mathbb P \\{D_n^- > t\\} \\le C \\exp(- 2 n t^2)\n",
    "\\end{equation*}\n",
    "for some constant $C$.\n",
    "\n",
    "Massart (1990) showed that $C = 1$ is the sharp constant in the DKW inequality, provided \n",
    "$\\exp(-2 n t^2) \\le \\frac{1}{2}$.\n",
    "The distribution of $D_n^-$ is stochastically larger when \n",
    "$F$ is continuous than when $F$ has jumps (Massart 1990); thus the inequality conservative for IID sampling from discrete or \"mixed\" distributions.\n",
    "\n",
    "Moreover, $D_n^-$ is stochastically larger for sampling with replacement than for sampling \n",
    "without replacement, so the inequality is conservative for sampling from a finite population\n",
    "without replacement as well.\n",
    "\n",
    "Let's call the inequality with $C=1$ the Massart-Dvoretsky-Kiefer-Wolfowitz (MDKW) inequality.\n",
    "It follows from the MDKW inequality that\n",
    "\\begin{equation*}  \n",
    "    \\mathbb P \\left \\{\\sup_x (F(x) - \\hat{F}_n(x)) > \\sqrt{-\\frac{\\ln \\alpha}{2n}} \\right \\} \\le \\alpha\n",
    "\\end{equation*}\n",
    "and (by symmetry) that\n",
    "\\begin{equation*}  \n",
    "    \\mathbb P \\left \\{\\sup_x |F(x) - \\hat{F}_n(x)| > 2\\sqrt{-\\frac{\\ln \\alpha}{2n}} \\right \\} \\le \\alpha.\n",
    "\\end{equation*}\n",
    "\n",
    "This lets us construct goodness-of-fit tests for models for real-valued data (in contrast, for instance, to models for categorical data or models for multidimensional data--although there are related methods for such models)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Testing whether data come from a uniform distributution on $[a, b]$\n",
    "\n",
    "To find a $P$-value for the hypothesis that $\\{X_j \\}_{j=1}^n$ are IID $U[a,b]$, calculate\n",
    "\n",
    "\\begin{equation*}  D_n = \\sup_{x \\in [a, b]} \\left | \\frac{x-a}{b-a} - \\hat{F}_n(x) \\right | .\\end{equation*}\n",
    "\n",
    "(Note that the maximum will be attained at or infinitesimally before some $X_j$.)\n",
    "\n",
    "The $P$-value can be found by solving\n",
    "\n",
    "\\begin{equation*} \n",
    " D_n = 2\\sqrt{-\\frac{\\ln P}{2n}}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*} \n",
    " -2n \\left ( D_n/2  \\right )^2 = \\ln P\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*} \n",
    "  P = \\exp(-2n \\left ( D_n/2 \\right )^2)\n",
    "\\end{equation*}\n",
    "\n",
    "(Of course, if $\\min_j X_j < a$ or $\\max_j X_j > b$, we know with certainty that\n",
    "$\\{X_j\\}$ are not IID $U[a,b]$.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Testing whether data come from a Poisson process\n",
    "\n",
    "This model has a parameter, but that will turn out not to matter, because our inference will be conditional on the observed number of events.\n",
    "\n",
    "We observe $X$ on the interval $[0, T]$. \n",
    "If $X$ has a Poisson process with rate $\\lambda$, then, conditional on the number of events that occur between $0$ and $T$, the times of the events have the same distribution as the order statistics of that many IID uniform random variables on $[0, T]$. \n",
    "That is, if $X(T) = n$, then the times of the $n$ events are like\n",
    "an IID sample of size $n$ from a $U[0, T]$ distribution.\n",
    "Thus, to test whether $X$ is indeed Poisson, we can use the previous test, taking $a=0$ and $b=T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Testing whether data come from a Poisson distribution with known rate $\\lambda$\n",
    "\n",
    "The CDF of the Poisson distribution with mean $\\lambda$ is \n",
    "\\begin{equation*}   e^{-\\lambda} \\sum_{j=0}^{\\lfloor x \\rfloor} \\lambda^j/j! \\end{equation*}\n",
    "\n",
    "This can be used with the KS test to test whether data come from a Poisson distribution with mean $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multinomial Models and the chi-square Goodness-of-Fit Test\n",
    "\n",
    "A very common test for goodness-of-fit is based on the probability that the data fall in different categories.\n",
    "\n",
    "The most natural setting is multinomial models, but (typically by binning the data), the test can be used in many situations--although in its usual form, it is an approximate test, while the KS test is conservative.\n",
    "\n",
    "Recall that the joint distribution of the number of values in each of $k \\ge 2$ categories\n",
    "for $n$ IID draws with probability $\\pi_j$ of selecting value $j$ in each draw is _multinomial_. \n",
    "A random vector $(X_1, \\ldots, X_k)$ has a multinomial joint distribution with parameters $n$\n",
    "and $\\{\\pi_j\\}_{j=1}^k$ iff\n",
    "\\begin{equation*}  \\Pr \\{X_j = x_j \\} = \\prod_{j=1}^k \\pi_j^{x_j} \\frac{n!}{x_1!x_2! \\cdots x_j!}, \\;\\; x_j \\ge 0,\\;\\; \\sum_{j=1}^k x_j = n.\\end{equation*}\n",
    "\n",
    "The chi-square statistic for categorical data is\n",
    "\n",
    "\\begin{equation*} \n",
    "  \\chi^2 \\equiv \\sum_{j=1}^k \\frac{(X_j - E_j)^2}{E_j},\n",
    "\\end{equation*}\n",
    "\n",
    "where $E_j = n \\pi_j$ is the expected value of $X_j$.\n",
    "\n",
    "If $(X_1, \\ldots, X_k)$ are jointly multinomial with category probabilities $\\{ \\pi_j \\}_{j=1}^k$, the probability distribution of\n",
    "$\\chi^2$ is asymptotically the chi-square distribution with $k-1$ degrees of freedom, as $E_j$ approach infinity.\n",
    "The approximation improves as $E_j$ grows; see [SticiGui: The Multinomial Distribution and the Chi-Squared Test for Goodness of Fit](https://www.stat.berkeley.edu/~stark/SticiGui/Text/chiSquare.htm).\n",
    "\n",
    "When $\\{E_j\\}$ are not all large (and even when they are), the null distribution of $\\chi^2$ can be approximated by simulation; the theoretical approximation can be pretty bad.\n",
    "\n",
    "The chi-square test is often applied to artificial categories created by binning continuous\n",
    "data. \n",
    "That makes the test results depend on the (arbitrary) choice of bins.\n",
    "The $P$-values will generally be misleading if the bins are chosen after looking at the data: they should be specified in advance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Testing whether data come from a Poisson distribution with known rate $\\lambda$\n",
    "\n",
    "A different approach to testing whether data come from a Poisson distribution is to bin the possible outcomes $\\{0, 1, 2, \\ldots\\}$ into a finite number $k$ of categories, then rely on the fact that the joint distribution of the number of observations that fall into each of the bins has a multinomial distribution. \n",
    "\n",
    "For instance, one might define four bins, with the first bin containing the possible outcomes $\\{0, 1, 2\\}$, the second containing $\\{3, 4, 5, 6\\}$, \n",
    "the third containing $\\{7, 8, 9\\}$, and the fourth containing $\\{10, 11, \\ldots \\}$. Then $\\pi_1 = e^{-\\lambda} \\sum_{\\ell = 0}^2 \\lambda^\\ell/\\ell!$, $\\pi_2 = e^{-\\lambda} \\sum_{\\ell = 3}^6 \\lambda^\\ell/\\ell!$, $\\pi_3 = e^{-\\lambda} \\sum_{\\ell = 7}^9 \\lambda^\\ell/\\ell!$, and $\\pi_4 = 1 - e^{-\\lambda} \\sum_{\\ell = 0}^9 \\lambda^\\ell/\\ell!$.\n",
    "\n",
    "If you are doing this and intend to use the chi-square approximation to the null distribution of the chi-square statistic, you should define the bins in a way that ensures the expected count in every bin is large. If the sample size and/or some $\\pi_j$ is too small to have $E_j$ large for all $j$, don't use the chi-square approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chi-square GOF when there are estimated parameters\n",
    "\n",
    "If a model has $p$ free parameters fitted to the data, the distribution of the chi-square statistic is still asymptotically the chi-square distribution, but with $k-p-1$ degrees of freedom instead of $k-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Poisson Regression\n",
    "\n",
    "TBD"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}