Notebooks/Resampling/Resampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resampling Methods\n",
    "### by [Richard W. Evans](https://sites.google.com/site/rickecon/), February 2020\n",
    "The code in this Jupyter notebook was written using Python 3.7. It uses data files [`Titanic dataset`](https://raw.githubusercontent.com/BigDataGal/Python-for-Data-Science/master/titanic-train.csv). For the code to run properly, you will either need to have access to the internet or you should have the data file in the same folder as the Jupyter notebook file. Otherwise, you will have to change the respective lines of the code that read in the data to reflect the location of that data.\n",
    "\n",
    "Resampling methods are a way to test the sensitivity of statistical results to estimation using a different sample. It is often too difficult or too expensive to draw a new sample from the population. Resampling methods take advantage of the training-set test-set paradigm to evaluate the sensitivity of estimates to sample variance. The two main classes of resampling methods are:\n",
    "\n",
    "1. Cross validation\n",
    "2. Bootstrapping\n",
    "\n",
    "In choosing models to predict or match data or to infer relationships between variables, James, et al (2013) decompose the process into *model assessment* and *model selection*. Model assessment is treated in this notebook. It is the process and various means of evaluating the fit or accuracy of a given model. Model selection is the process of adjusting parameters, variables, or functional relationships between variables to better fit the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Cross validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Validation set approach\n",
    "This is an approach that we have already studied in the [classifiers notebook](https://github.com/UC-MACSS/persp-model-econ_W19/blob/master/Notebooks/Classification/LogitKNN.ipynb).\n",
    "\n",
    "1. Partition the data into a training set and a test set.\n",
    "2. Estimate the model using the training set.\n",
    "3. Evaluate the fit or predictive accuracy on the test set.\n",
    "\n",
    "The primary measure of fit is the mean squared error (MSE) of the estimated model on the test set. Let the test set have $N$ observations. The MSE of the test set is the sum of squared deviations of the actual dependent variable values minus the predicted values.\n",
    "\n",
    "$$ MSE = \\frac{1}{N}\\sum_{i=1}^N\\left(y_i - \\hat{y}_i\\right)^2 $$\n",
    "\n",
    "In classification problems, researchers sometimes use the ROC curve (receiver operating characteristics) and the area under the ROC curve (AUROC) because it captures both type 1 and type 2 errors in a single visualization (false positive rate versus true positive rate). You want the false positive rate to be low for any given true positive rate. The area under the ROC curve (AUROC). The AUROC is the opposite of a measure of error. The bigger the AUROC, the more accurate predictor is the model.\n",
    "\n",
    "[Insert figure]\n",
    "\n",
    "Let's calculate the MSE from our logistic regression of the titanic example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import needed stuff\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandas import Series, DataFrame\n",
    "\n",
    "import sklearn\n",
    "from sklearn import preprocessing\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import train_test_split, LeaveOneOut, KFold\n",
    "from sklearn import metrics \n",
    "from sklearn.metrics import classification_report, mean_squared_error\n",
    "from pylab import rcParams\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize'] = 10, 8\n",
    "sb.set_style('whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in Titanic data\n",
    "url = ('https://raw.githubusercontent.com/BigDataGal/Python-for-Data-Science/' +\n",
    "      'master/titanic-train.csv')\n",
    "titanic = pd.read_csv(url)\n",
    "titanic.columns = ['PassengerId','Survived','Pclass','Name','Sex','Age',\n",
    "                   'SibSp','Parch','Ticket','Fare','Cabin','Embarked']\n",
    "\n",
    "# Get rid of columns we don't use\n",
    "titanic = titanic.drop(['PassengerId','Name','Ticket','Cabin'], 1)\n",
    "\n",
    "# Impute missing age values\n",
    "def age_approx(cols):\n",
    "    Age = cols[0]\n",
    "    Pclass = cols[1]\n",
    "    \n",
    "    if pd.isnull(Age):\n",
    "        if Pclass == 1:\n",
    "            return 37\n",
    "        elif Pclass == 2:\n",
    "            return 29\n",
    "        else:\n",
    "            return 24\n",
    "    else:\n",
    "        return Age\n",
    "\n",
    "titanic['Age'] = \\\n",
    "    titanic[['Age', 'Pclass']].apply(age_approx, axis=1)\n",
    "    \n",
    "# Drop any observations with missing values\n",
    "titanic.dropna(inplace=True)\n",
    "\n",
    "# Make gender dummies and embark dummies and get rid of\n",
    "# original variables\n",
    "gender = pd.get_dummies(titanic['Sex'], drop_first=True)\n",
    "embark_location = pd.get_dummies(titanic['Embarked'],\n",
    "                                 drop_first=True)\n",
    "titanic.drop(['Sex', 'Embarked'], axis=1, inplace=True)\n",
    "titanic = pd.concat([titanic, gender, embark_location], axis=1)\n",
    "\n",
    "# Drop Pclass variable due to excessive correlation with Fare\n",
    "titanic.drop(['Pclass'], axis=1, inplace=True)\n",
    "\n",
    "titanic.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now partition the data into the same 60% training set sample that we did in the [logistic regression notebook](https://github.com/UC-MACSS/persp-model_W18/blob/master/Notebooks/Classfcn1/KKNlogitLDA.ipynb) and estimate the logistic regression with all the variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = titanic[['Age', 'SibSp', 'Parch', 'Fare', 'male', 'Q', 'S']]\n",
    "y = titanic['Survived']\n",
    "# This function train_test_split is from sklearn.cross_validation\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4,\n",
    "                                                    random_state=25)\n",
    "LogReg = LogisticRegression(max_iter=200)\n",
    "LogReg.fit(X_train, y_train)\n",
    "y_pred = LogReg.predict(X_test)\n",
    "# Note that the squared doesn't matter in a Logistic model\n",
    "\n",
    "# You can code the MSE yourself\n",
    "MSE_vs = ((y_test - y_pred) ** 2).sum() / y_pred.shape[0]\n",
    "print('Validation set MSE = ', MSE_vs)\n",
    "\n",
    "# Or you can use scikit-learn's method\n",
    "print('Validation set MSE = ', mean_squared_error(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the estimation of the logistic regression model on all the observations in the dataset, not just the training set observations. What advantages are there to estimating on a subsample of the data (training set) rather than the full data set?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LogReg1 = LogisticRegression(max_iter=200)\n",
    "LogReg1.fit(X, y)\n",
    "y1_pred = LogReg1.predict(X)\n",
    "\n",
    "# You can code the MSE yourself\n",
    "MSE1_vs = ((y - y1_pred) ** 2).sum() / y1_pred.shape[0]\n",
    "print('Validation set MSE = ', MSE1_vs)\n",
    "\n",
    "# Or you can use scikit-learn's method\n",
    "print('Validation set MSE = ', mean_squared_error(y, y1_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Leave-one-out cross validation\n",
    "Leave-one-out cross validation (LOOCV) is an approach in which the model is assessed using $N$ different training sets and test sets of a specific size. Let the data have $N$ observations. LOOCV is to choose a training set with $N-1$ observations, such that the test set only has one observation $y_i$. Repeat this $N$ with a slightly different training set such that each data point is the test set in exactly one of these sebsets.\n",
    "\n",
    "In this case, the mean squared error MSE has no summation because there is only one observation in the test set.\n",
    "\n",
    "$$ MSE_i = (y_i - \\hat{y}_i)^2 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LOOCV estimate for the test MSE is the average of these $N$ test error estimates.\n",
    "\n",
    "$$ CV_{loo} = \\frac{1}{N}\\sum_{i=1}^N MSE_i $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define loo as a leave-one-out object, then\n",
    "# split it into N different partitions\n",
    "\n",
    "# Note that the LeaveOneOut() function does not work\n",
    "# well with pandas DataFrames\n",
    "Xvars = titanic[['Age', 'SibSp', 'Parch', 'Fare',\n",
    "                 'male', 'Q', 'S']].values\n",
    "yvars = titanic['Survived'].values\n",
    "N_loo = Xvars.shape[0]\n",
    "loo = LeaveOneOut()\n",
    "loo.get_n_splits(Xvars)\n",
    "MSE_vec = np.zeros(N_loo)\n",
    "\n",
    "# This loop will take 20 or 30 seconds\n",
    "for train_index, test_index in loo.split(Xvars):\n",
    "    X_train, X_test = Xvars[train_index], Xvars[test_index]\n",
    "    y_train, y_test = yvars[train_index], yvars[test_index]\n",
    "    LogReg = LogisticRegression(max_iter=200)\n",
    "    LogReg.fit(X_train, y_train)\n",
    "    y_pred = LogReg.predict(X_test)\n",
    "    MSE_vec[test_index] = (y_test - y_pred) ** 2\n",
    "    print('MSE for test set', test_index, ' is', MSE_vec[test_index])\n",
    "\n",
    "MSE_loo = MSE_vec.mean()\n",
    "MSE_loo_std = MSE_vec.std()\n",
    "print('test estimate MSE loocv=', MSE_loo,\n",
    "      ', test estimate MSE standard err=', MSE_loo_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3. k-fold cross validation\n",
    "$k$-fold cross validation is a method in which the dataset is randomly divided into $k$ groups (folds). Define a test set of the model as the $k$th fold. For each test set $k$, the model is estimated on the data from the other $k-1$ folds. Let the number of observations in the $k$th fold be $N_k$, and let $\\mathcal{K}$ be the set of observations in the $k$th fold. The $MSE_k$ of the $k$th fold is:\n",
    "\n",
    "$$ MSE_k = \\frac{1}{N_k}\\sum_{i\\in\\mathcal{K}}(y_i - \\hat{y}_i)^2 $$\n",
    "\n",
    "Then the $k$-fold estimate for the test MSE is the average of these $k$ test error estimates.\n",
    "\n",
    "$$ CV_{kf} = \\frac{1}{k}\\sum_{j=1}^k MSE_j $$\n",
    "\n",
    "LOOCV is a special case of $k$-fold cross validation in which $k=N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the Titanic data again and test our logit model performance with a $k$-fold cross validation with $k=6$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 800\n",
    "kf = KFold(n_splits=k, random_state=10, shuffle=True)\n",
    "kf.get_n_splits(Xvars)\n",
    "\n",
    "MSE_vec_kf = np.zeros(k)\n",
    "\n",
    "k_ind = int(0)\n",
    "for train_index, test_index in kf.split(Xvars):\n",
    "    # print(\"TRAIN:\", train_index, \"TEST:\", test_index)\n",
    "    # print('k index=', k_ind)\n",
    "    X_train, X_test = Xvars[train_index], Xvars[test_index]\n",
    "    y_train, y_test = yvars[train_index], yvars[test_index]\n",
    "    LogReg = LogisticRegression(max_iter=300)\n",
    "    LogReg.fit(X_train, y_train)\n",
    "    y_pred = LogReg.predict(X_test)\n",
    "    MSE_vec_kf[k_ind] = ((y_test - y_pred) ** 2).mean()\n",
    "    # print('MSE for test set', k_ind, ' is', MSE_vec_kf[k_ind])\n",
    "    k_ind += 1\n",
    "\n",
    "MSE_kf = MSE_vec_kf.mean()\n",
    "MSE_kf_std = MSE_vec_kf.std()\n",
    "print('test estimate MSE k-fold=', MSE_kf,\n",
    "      'test estimate MSE standard err=', MSE_kf_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4. Bias versus variance\n",
    "Recall the test estimate MSE from the LOOCV of approximately 0.2115 and the MSE(LOOCV) standard error of about 0.4084. What happens to the estimated MSE and MSE standard error in the $k$-fold cross validation above as $k$ increases? Try values of $k=2, 10, 50, 100, 800$.\n",
    "\n",
    "Note that the LOOCV method has low bias (estimated on large number of data) but high variance (errors are based on one draw). In contrast, the $k$-fold method has more bias (estimated with less data) but lower variance. Each test set has more observations.\n",
    "\n",
    "* $k$-fold cross validation can often provide more accurate estimates of the test error rate.\n",
    "* $k$-fold is less computationally intensive\n",
    "* LOOCV has the least bias\n",
    "* LOOCV is the most computationally expensive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Bootstrapping\n",
    "This name comes from the expression \"to pull oneself up by ones own bootstraps.\" In a way similar to the cross validation methods of the last section, we can use *the bootstrap* to quantify the undertainty associated with a given estimator, learning model, or method. In the econometrics and statistics literature, this often shows up as \"bootstrapped standard errors\". Bootstrapping is valuable because it is so widely applicable to a range of models.\n",
    "\n",
    "1. Randomly draw $S$ datasets of size $N_S$ with replacement. Define each training set of observations as $\\mathcal{K}_s$ and each corresponding test set as $\\mathcal{-K}_{s}$.\n",
    "2. Calculate the MSE for each test set $\\mathcal{-K}_{s}$\n",
    "\n",
    "The bootstrap estimate for the test MSE is the average MSE from each random test set.\n",
    "\n",
    "$$ CV_{boot} = \\frac{1}{S}\\sum_{s=1}^S MSE_s $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_bs = 800\n",
    "\n",
    "MSE_vec_bs = np.zeros(N_bs)\n",
    "\n",
    "for bs_ind in range(N_bs):\n",
    "    X_train, X_test, y_train, y_test = \\\n",
    "        train_test_split(X, y, test_size=0.4)\n",
    "    LogReg = LogisticRegression(max_iter=200)\n",
    "    LogReg.fit(X_train, y_train)\n",
    "    y_pred = LogReg.predict(X_test)\n",
    "    MSE_vec_bs[bs_ind] = ((y_test - y_pred) ** 2).mean()\n",
    "    print('MSE for test set', bs_ind, ' is', MSE_vec_bs[bs_ind])\n",
    "\n",
    "MSE_bs = MSE_vec_bs.mean()\n",
    "MSE_bs_std = MSE_vec_bs.std()\n",
    "print('test estimate MSE bootstrap=', MSE_bs,\n",
    "      'test estimate MSE standard err=', MSE_bs_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. How to use cross validation for model assessment and selection\n",
    "*Model assessment* is the process of evaluating the performance of a particular model estimated on training data on its prediction accuracy on test data. There are many criteria for model assessment. The most common measure of model accuracy on test data is the mean squared error $MSE$ or root mean squared error $rMSE$. However, we have seen that the measure $MSE$ varies depending on which cross validation method is used.\n",
    "\n",
    "[JWHT13] define *Model selection* as the process of \"selecting the proper level of flexibility for a model. That is, they define model selection as a process of tuning a particular family or class of models to maximize accuracy on test set prediction. So model selection involves model assessment. However, one can expand this definition of model selection to include testing multiple families or classes of model in terms of accuracy--a horse race.\n",
    "\n",
    "The narrower [JWHT13] definition of model selection is analogous to maximizing the efficiency of a particular horse and rider in a horse race in which the horse is racing against itself. My broader definition of model selection is analogous running a number of horses in a race, with each horse being optimized for efficiency. In this broader definition, many variables are at play and the data must be sampled many times (think of cross-validation optimization on each horse [model]).\n",
    "\n",
    "This broader definition of model selection is computationally intensive. But that is where TensorFlow shines. [TensorFlow](https://www.tensorflow.org/) is an open source software library developed by the Google Brain AI group. Broadly, TensorFlow is a system of libraries that facilitate efficient, parallel, and scalable use of available processors (CPUs and GPUs) as well as memory management. For statistical learning, TensorFlow is optimized to efficiently run model assessment and model selection algorithms.\n",
    "\n",
    "Many good empirical papers run a horse race on tuned statistical learning models to maximize predictive accuracy. It is hard to know *ex ante* which model will be the most accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Gopalan, \"Predicting Infant Mortality: Minimizing False Negatives\"\n",
    "[G18] is able to maintain overall accuracy, and decrease false negatives from 74% to 7%. She uses a regularization method on one of her variables (class), called \"Tomek links\", to make the variable more informative. She then tests five different predictive models (random forest, AdaBost, XGBoost, decision tree, logistic regression). She also tried a couple of additional data transformations. In her case, the random forest model had the best performance.\n",
    "\n",
    "**[G13, Table 6] Comparison of Classifiers: Tomek Links.** FNR=false negative rate, FPR=False positive rate, AUROC=area under the ROC curve.\n",
    "\n",
    "| Model | Train-FNR | Train-FPR  | Train-AUROC | Test-FNR  | Test-FPR  | Test-AUROC |\n",
    "| --- | --- | --- | --- | --- | --- | --- |\n",
    "| Random Forest | 0.06 | 0.06 | 0.62  | 0.07 | 0.06 | 0.62  |\n",
    "| AdaBoost      | 0.40 | 0.07 | 0.51  | 0.45 | 0.07 | 0.50  |\n",
    "| XGBoost       | 0.40 | 0.07 | 0.51  | 0.43 | 0.07 | 0.51  |\n",
    "| Decision Tree | 0.03 | 0.06 | 0.62  | 0.07 | 0.06 | 0.62  |\n",
    "| Logistic Reg  | 0.41 | 0.07 | 0.51  | 0.43 | 0.07 | 0.65  |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Problem with full-sample estimation\n",
    "Overfitting. Maximizing accuracy in training set can pick incorporate noise from the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. References\n",
    "* Gopalan, Sushmita V., \"Predicting Infant Mortality: Minimizing False Negatives,\" MACSS Thesis, University of Chicago (April 2018).\n",
    "* James, Gareth, Deaniela Witten, Trevor Hastie, and Robert Tibshirani, [*An Introduction to Statistical Learning with Applications in R*](http://link.springer.com.proxy.uchicago.edu/book/10.1007%2F978-1-4614-7138-7), New York, Springer (2013)."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}