Ch05-resample-lab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "85ad9863",
   "metadata": {},
   "source": [
    "\n",
    "# Chapter 5\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac8b08af",
   "metadata": {},
   "source": [
    "# Lab: Cross-Validation and the Bootstrap\n",
    "In this lab, we explore the resampling techniques covered in this\n",
    "chapter. Some of the commands in this lab may take a while to run on\n",
    "your computer.\n",
    "\n",
    "We again begin by placing most of our imports at this top level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e7712cfe",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from ISLP import load_data\n",
    "from ISLP.models import (ModelSpec as MS,\n",
    "                         summarize,\n",
    "                         poly)\n",
    "from sklearn.model_selection import train_test_split\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "784a2ba3",
   "metadata": {},
   "source": [
    "There are several new imports needed for this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "21c2ed4f",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "from sklearn.model_selection import \\\n",
    "     (cross_validate,\n",
    "      KFold,\n",
    "      ShuffleSplit)\n",
    "from sklearn.base import clone\n",
    "from ISLP.models import sklearn_sm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac3acd5",
   "metadata": {},
   "source": [
    "## The Validation Set Approach\n",
    "We explore the use of the validation set approach in order to estimate\n",
    "the test error rates that result from fitting various linear models on\n",
    "the  `Auto`  data set.\n",
    "\n",
    "We use the function `train_test_split()` to split\n",
    "the data into training and validation sets. As there are 392 observations,\n",
    "we split into two equal sets of size 196 using the\n",
    "argument `test_size=196`. It is generally a good idea to set a random seed\n",
    "when performing operations like this that contain an\n",
    "element of randomness, so that the results obtained can be reproduced\n",
    "precisely at a later time. We set the random seed of the splitter\n",
    "with the argument `random_state=0`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8af59641",
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "Auto = load_data('Auto')\n",
    "Auto_train, Auto_valid = train_test_split(Auto,\n",
    "                                         test_size=196,\n",
    "                                         random_state=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76383f0",
   "metadata": {},
   "source": [
    "Now we can fit a linear regression using only the observations corresponding to the training set `Auto_train`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d9b0b7c8",
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "hp_mm = MS(['horsepower'])\n",
    "X_train = hp_mm.fit_transform(Auto_train)\n",
    "y_train = Auto_train['mpg']\n",
    "model = sm.OLS(y_train, X_train)\n",
    "results = model.fit()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d196dd08",
   "metadata": {},
   "source": [
    "We now use the `predict()` method of `results` evaluated on the model matrix for this model\n",
    "created using the validation data set. We also calculate the validation MSE of our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3e77d831",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23.61661706966988"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_valid = hp_mm.transform(Auto_valid)\n",
    "y_valid = Auto_valid['mpg']\n",
    "valid_pred = results.predict(X_valid)\n",
    "np.mean((y_valid - valid_pred)**2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4369ee6",
   "metadata": {},
   "source": [
    "Hence our estimate for the validation MSE of  the linear regression\n",
    "fit is $23.62$.\n",
    "\n",
    "We can also estimate the validation error for\n",
    "higher-degree polynomial regressions. We first provide a function `evalMSE()` that takes a model string as well\n",
    "as a training and test set and returns the MSE on the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0aa4bfcc",
   "metadata": {
    "execution": {}
   },
   "outputs": [],
   "source": [
    "def evalMSE(terms,\n",
    "            response,\n",
    "            train,\n",
    "            test):\n",
    "\n",
    "   mm = MS(terms)\n",
    "   X_train = mm.fit_transform(train)\n",
    "   y_train = train[response]\n",
    "\n",
    "   X_test = mm.transform(test)\n",
    "   y_test = test[response]\n",
    "\n",
    "   results = sm.OLS(y_train, X_train).fit()\n",
    "   test_pred = results.predict(X_test)\n",
    "\n",
    "   return np.mean((y_test - test_pred)**2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0271dc50",
   "metadata": {},
   "source": [
    "Let’s use this function to estimate the validation MSE\n",
    "using linear, quadratic and cubic fits. We use the `enumerate()`  function\n",
    "here, which gives both the values and indices of objects as one iterates\n",
    "over a for loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a0dbd55f",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([23.61661707, 18.76303135, 18.79694163])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE = np.zeros(3)\n",
    "for idx, degree in enumerate(range(1, 4)):\n",
    "    MSE[idx] = evalMSE([poly('horsepower', degree)],\n",
    "                       'mpg',\n",
    "                       Auto_train,\n",
    "                       Auto_valid)\n",
    "MSE\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7401536",
   "metadata": {},
   "source": [
    "These error rates are $23.62, 18.76$, and $18.80$, respectively. If we\n",
    "choose a different training/validation split instead, then we\n",
    "can expect somewhat different errors on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "885136a4",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20.75540796, 16.94510676, 16.97437833])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Auto_train, Auto_valid = train_test_split(Auto,\n",
    "                                          test_size=196,\n",
    "                                          random_state=3)\n",
    "MSE = np.zeros(3)\n",
    "for idx, degree in enumerate(range(1, 4)):\n",
    "    MSE[idx] = evalMSE([poly('horsepower', degree)],\n",
    "                       'mpg',\n",
    "                       Auto_train,\n",
    "                       Auto_valid)\n",
    "MSE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00785402",
   "metadata": {},
   "source": [
    "Using this split of the observations into a training set and a validation set,\n",
    "we find that the validation set error rates for the models with linear, quadratic, and cubic terms are $20.76$, $16.95$, and $16.97$, respectively.\n",
    "\n",
    "These results are consistent with our previous findings: a model that\n",
    "predicts `mpg` using a quadratic function of `horsepower`\n",
    "performs better than a model that involves only a linear function of\n",
    "`horsepower`, and there is no evidence of an improvement in using a cubic function of `horsepower`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c071b8",
   "metadata": {},
   "source": [
    "## Cross-Validation\n",
    "In theory, the cross-validation estimate can be computed for any generalized\n",
    "linear model.  {}\n",
    "In practice, however, the simplest way to cross-validate in\n",
    "Python is to use `sklearn`, which has a different interface or API\n",
    "than `statsmodels`, the code we have been using to fit GLMs.\n",
    "\n",
    "This is a problem which often confronts data scientists: \"I have a function to do task $A$, and need to feed it into something that performs task $B$, so that I can compute $B(A(D))$, where $D$ is my data.\" When $A$ and $B$ don’t naturally speak to each other, this\n",
    "requires the use of a *wrapper*.\n",
    "In the `ISLP` package,\n",
    "we provide \n",
    "a wrapper, `sklearn_sm()`, that enables us to easily use the cross-validation tools of `sklearn` with\n",
    "models fit by `statsmodels`.\n",
    "\n",
    "The class `sklearn_sm()` \n",
    "has  as its first argument\n",
    "a model from `statsmodels`. It can take two additional\n",
    "optional arguments: `model_str` which can be\n",
    "used to specify a formula, and `model_args` which should\n",
    "be a dictionary of additional arguments used when fitting\n",
    "the model. For example, to fit a logistic regression model\n",
    "we have to specify a `family` argument. This\n",
    "is passed as `model_args={'family':sm.families.Binomial()}`.\n",
    "\n",
    "Here is our wrapper in action:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6d957d8c",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24.23151351792924"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp_model = sklearn_sm(sm.OLS,\n",
    "                      MS(['horsepower']))\n",
    "X, Y = Auto.drop(columns=['mpg']), Auto['mpg']\n",
    "cv_results = cross_validate(hp_model,\n",
    "                            X,\n",
    "                            Y,\n",
    "                            cv=Auto.shape[0])\n",
    "cv_err = np.mean(cv_results['test_score'])\n",
    "cv_err\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c17e2bc8",
   "metadata": {},
   "source": [
    "The arguments to `cross_validate()` are as follows: an\n",
    "object with the appropriate `fit()`, `predict()`,\n",
    "and `score()` methods,  an\n",
    "array of features `X` and a response `Y`. \n",
    "We also included an additional argument `cv` to `cross_validate()`; specifying an integer\n",
    "$K$ results in $K$-fold cross-validation. We have provided a value \n",
    "corresponding to the total number of observations, which results in\n",
    "leave-one-out cross-validation (LOOCV). The `cross_validate()`  function produces a dictionary with several components;\n",
    "we simply want the cross-validated test score here (MSE), which is estimated to be 24.23."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c7901f2",
   "metadata": {},
   "source": [
    "We can repeat this procedure for increasingly complex polynomial fits.\n",
    "To automate the process, we again\n",
    "use a for loop which iteratively fits polynomial\n",
    "regressions of degree 1 to 5, computes the\n",
    "associated cross-validation error, and stores it in the $i$th element\n",
    "of the vector `cv_error`. The variable `d` in the for loop\n",
    "corresponds to the degree of the polynomial. We begin by initializing the\n",
    "vector. This command may take a couple of seconds to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e2b5ce95",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([24.23151352, 19.24821312, 19.33498406, 19.42443033, 19.03323827])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cv_error = np.zeros(5)\n",
    "H = np.array(Auto['horsepower'])\n",
    "M = sklearn_sm(sm.OLS)\n",
    "for i, d in enumerate(range(1,6)):\n",
    "    X = np.power.outer(H, np.arange(d+1))\n",
    "    M_CV = cross_validate(M,\n",
    "                          X,\n",
    "                          Y,\n",
    "                          cv=Auto.shape[0])\n",
    "    cv_error[i] = np.mean(M_CV['test_score'])\n",
    "cv_error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03706248",
   "metadata": {},
   "source": [
    "As in Figure 5.4, we see a sharp drop in the estimated test MSE between the linear and\n",
    "quadratic fits, but then no clear improvement from using higher-degree polynomials.\n",
    "\n",
    "Above we introduced the `outer()`  method of the `np.power()`\n",
    "function.  The `outer()` method is applied to an operation\n",
    "that has two arguments, such as `add()`, `min()`, or\n",
    "`power()`.\n",
    "It has two arrays as\n",
    "arguments, and then forms a larger\n",
    "array where the operation is applied to each pair of elements of the\n",
    "two arrays. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1dda1bd7",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 5,  7],\n",
       "       [ 7,  9],\n",
       "       [11, 13]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([3, 5, 9])\n",
    "B = np.array([2, 4])\n",
    "np.add.outer(A, B)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5092f1b",
   "metadata": {},
   "source": [
    "In the CV example above, we used $K=n$, but of course we can also use $K<n$. The code is very similar\n",
    "to the above (and is significantly faster). Here we use `KFold()` to partition the data into $K=10$ random groups. We use `random_state` to set a random seed and initialize a vector `cv_error` in which we will store the CV errors corresponding to the\n",
    "polynomial fits of degrees one to five."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fb25fa70",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([24.20766449, 19.18533142, 19.27626666, 19.47848403, 19.13720581])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cv_error = np.zeros(5)\n",
    "cv = KFold(n_splits=10,\n",
    "           shuffle=True,\n",
    "           random_state=0) # use same splits for each degree\n",
    "for i, d in enumerate(range(1,6)):\n",
    "    X = np.power.outer(H, np.arange(d+1))\n",
    "    M_CV = cross_validate(M,\n",
    "                          X,\n",
    "                          Y,\n",
    "                          cv=cv)\n",
    "    cv_error[i] = np.mean(M_CV['test_score'])\n",
    "cv_error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ec6afb",
   "metadata": {},
   "source": [
    "Notice that the computation time is much shorter than that of LOOCV.\n",
    "(In principle, the computation time for LOOCV for a least squares\n",
    "linear model should be faster than for $K$-fold CV, due to the\n",
    "availability of the formula (5.2)  for LOOCV;\n",
    "however, the generic `cross_validate()`  function does not make\n",
    "use of this formula.)  We still see little evidence that using cubic\n",
    "or higher-degree polynomial terms leads to a lower test error than simply\n",
    "using a quadratic fit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5edf407f",
   "metadata": {},
   "source": [
    "The `cross_validate()` function is flexible and can take\n",
    "different splitting mechanisms as an argument. For instance, one can use the `ShuffleSplit()` funtion to implement\n",
    "the validation set approach just as easily as K-fold cross-validation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d78795cd",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([23.61661707])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "validation = ShuffleSplit(n_splits=1,\n",
    "                          test_size=196,\n",
    "                          random_state=0)\n",
    "results = cross_validate(hp_model,\n",
    "                         Auto.drop(['mpg'], axis=1),\n",
    "                         Auto['mpg'],\n",
    "                         cv=validation);\n",
    "results['test_score']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a081be63",
   "metadata": {},
   "source": [
    "One can estimate the variability in the test error by running the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0407ad56",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(23.802232661034168, 1.421845094109185)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "validation = ShuffleSplit(n_splits=10,\n",
    "                          test_size=196,\n",
    "                          random_state=0)\n",
    "results = cross_validate(hp_model,\n",
    "                         Auto.drop(['mpg'], axis=1),\n",
    "                         Auto['mpg'],\n",
    "                         cv=validation)\n",
    "results['test_score'].mean(), results['test_score'].std()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b66db3cb",
   "metadata": {},
   "source": [
    "Note that this standard deviation is not a valid estimate of the\n",
    "sampling variability of the mean test score or the individual scores, since the randomly-selected training\n",
    "samples overlap and hence introduce correlations. But it does give an\n",
    "idea of the Monte Carlo variation\n",
    "incurred by picking different random folds.\n",
    "\n",
    "## The Bootstrap\n",
    "We illustrate the use of the bootstrap in the simple example\n",
    " {of Section 5.2,}  as well as on an example involving\n",
    "estimating the accuracy of the linear regression model on the  `Auto`\n",
    "data set.\n",
    "### Estimating the Accuracy of a Statistic of Interest\n",
    "One of the great advantages of the bootstrap approach is that it can\n",
    "be applied in almost all situations. No complicated mathematical\n",
    "calculations are required. While there are several implementations\n",
    "of the bootstrap in Python, its use for estimating\n",
    "standard error is simple enough that we write our own function\n",
    "below for the case when our data is stored\n",
    "in a dataframe.\n",
    "\n",
    "To illustrate the bootstrap, we\n",
    "start with a simple example.\n",
    "The  `Portfolio`  data set in the `ISLP` package is described\n",
    "in Section 5.2. The goal is to estimate the\n",
    "sampling variance of the parameter $\\alpha$ given in formula (5.7).  We will\n",
    "create a function\n",
    "`alpha_func()`, which takes as input a dataframe `D` assumed\n",
    "to have columns `X` and `Y`, as well as a\n",
    "vector `idx` indicating which observations should be used to\n",
    "estimate \n",
    "$\\alpha$. The function then outputs the estimate for $\\alpha$ based on\n",
    "the selected observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f04f15bd",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "Portfolio = load_data('Portfolio')\n",
    "def alpha_func(D, idx):\n",
    "   cov_ = np.cov(D[['X','Y']].loc[idx], rowvar=False)\n",
    "   return ((cov_[1,1] - cov_[0,1]) /\n",
    "           (cov_[0,0]+cov_[1,1]-2*cov_[0,1]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88bd6a4",
   "metadata": {},
   "source": [
    "This function returns an estimate for $\\alpha$\n",
    "based on applying the minimum\n",
    "    variance formula (5.7) to the observations indexed by\n",
    "the argument `idx`.  For instance, the following command\n",
    "estimates $\\alpha$ using all 100 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f98c0323",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.57583207459283"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_func(Portfolio, range(100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58a78f00",
   "metadata": {},
   "source": [
    "Next we randomly select\n",
    "100 observations from `range(100)`, with replacement. This is equivalent\n",
    "to constructing a new bootstrap data set and recomputing $\\hat{\\alpha}$\n",
    "based on the new data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bcd40175",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6074452469619004"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "alpha_func(Portfolio,\n",
    "           rng.choice(100,\n",
    "                      100,\n",
    "                      replace=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6058be4",
   "metadata": {},
   "source": [
    "This process can be generalized to create a simple function `boot_SE()` for\n",
    "computing the bootstrap standard error for arbitrary\n",
    "functions that take only a data frame as an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ab6602cd",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "def boot_SE(func,\n",
    "            D,\n",
    "            n=None,\n",
    "            B=1000,\n",
    "            seed=0):\n",
    "    rng = np.random.default_rng(seed)\n",
    "    first_, second_ = 0, 0\n",
    "    n = n or D.shape[0]\n",
    "    for _ in range(B):\n",
    "        idx = rng.choice(D.index,\n",
    "                         n,\n",
    "                         replace=True)\n",
    "        value = func(D, idx)\n",
    "        first_ += value\n",
    "        second_ += value**2\n",
    "    return np.sqrt(second_ / B - (first_ / B)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d94d383e",
   "metadata": {},
   "source": [
    "Notice the use of `_` as a loop variable in `for _ in range(B)`. This is often used if the value of the counter is\n",
    "unimportant and simply makes sure  the loop is executed `B` times.\n",
    "\n",
    "Let’s use our function to evaluate the accuracy of our\n",
    "estimate of $\\alpha$ using $B=1{,}000$ bootstrap replications. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4a323513",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09118176521277699"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_SE = boot_SE(alpha_func,\n",
    "                   Portfolio,\n",
    "                   B=1000,\n",
    "                   seed=0)\n",
    "alpha_SE\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22343f53",
   "metadata": {},
   "source": [
    "The final output shows that the bootstrap estimate for ${\\rm SE}(\\hat{\\alpha})$ is $0.0912$.\n",
    "\n",
    "### Estimating the Accuracy of a Linear Regression Model\n",
    "The bootstrap approach can be used to assess the variability of the\n",
    "coefficient estimates and predictions from a statistical learning\n",
    "method. Here we use the bootstrap approach in order to assess the\n",
    "variability of the estimates for $\\beta_0$ and $\\beta_1$, the\n",
    "intercept and slope terms for the linear regression model that uses\n",
    "`horsepower` to predict `mpg` in the  `Auto`  data set. We\n",
    "will compare the estimates obtained using the bootstrap to those\n",
    "obtained using the formulas for ${\\rm SE}(\\hat{\\beta}_0)$ and\n",
    "${\\rm SE}(\\hat{\\beta}_1)$ described in Section 3.1.2.\n",
    "\n",
    "To use our `boot_SE()` function, we must write a function (its\n",
    "first argument)\n",
    "that takes a data frame `D` and indices `idx`\n",
    "as its only arguments. But here we want to bootstrap a specific\n",
    "regression model, specified by a model formula and data. We show how\n",
    "to do this in a few simple steps.\n",
    "\n",
    "We start by writing a generic\n",
    "function `boot_OLS()` for bootstrapping a regression model that takes a formula to\n",
    "define the corresponding regression. We use the `clone()` function to\n",
    "make a copy of the formula that can be refit to the new dataframe. This means\n",
    "that any derived features such as those defined by `poly()`\n",
    "(which we will see shortly),\n",
    "will be re-fit on the resampled data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0220f3af",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "def boot_OLS(model_matrix, response, D, idx):\n",
    "    D_ = D.loc[idx]\n",
    "    Y_ = D_[response]\n",
    "    X_ = clone(model_matrix).fit_transform(D_)\n",
    "    return sm.OLS(Y_, X_).fit().params"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df0c7f05",
   "metadata": {},
   "source": [
    "This is not quite what is needed as the first argument to\n",
    "`boot_SE()`. The first two arguments which specify the model will not change in the\n",
    "bootstrap process, and we would like to *freeze* them.   The\n",
    "function `partial()` from the `functools` module  does precisely this: it takes a function\n",
    "as an argument, and freezes some of its arguments, starting from the\n",
    "left. We use it to freeze the first two model-formula arguments of `boot_OLS()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "62037dcb",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "hp_func = partial(boot_OLS, MS(['horsepower']), 'mpg')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61fbe248",
   "metadata": {},
   "source": [
    "Typing `hp_func?` will show that it has two arguments `D`\n",
    "and `idx` --- it is a version of `boot_OLS()` with the first\n",
    "two arguments frozen --- and hence is ideal as the first argument for `boot_SE()`.\n",
    "\n",
    "The `hp_func()` function can now be used in order to create\n",
    "bootstrap estimates for the intercept and slope terms by randomly\n",
    "sampling from among the observations with replacement. We first\n",
    "demonstrate its utility on 10 bootstrap samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b8bdb7a4",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[39.88064456, -0.1567849 ],\n",
       "       [38.73298691, -0.14699495],\n",
       "       [38.31734657, -0.14442683],\n",
       "       [39.91446826, -0.15782234],\n",
       "       [39.43349349, -0.15072702],\n",
       "       [40.36629857, -0.15912217],\n",
       "       [39.62334517, -0.15449117],\n",
       "       [39.0580588 , -0.14952908],\n",
       "       [38.66688437, -0.14521037],\n",
       "       [39.64280792, -0.15555698]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "np.array([hp_func(Auto,\n",
    "          rng.choice(392,\n",
    "                     392,\n",
    "                     replace=True)) for _ in range(10)])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a831036",
   "metadata": {},
   "source": [
    "Next, we use the `boot_SE()` {}  function to compute the standard\n",
    "errors of 1,000 bootstrap estimates for the intercept and slope terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "36808258",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept     0.848807\n",
       "horsepower    0.007352\n",
       "dtype: float64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp_se = boot_SE(hp_func,\n",
    "                Auto,\n",
    "                B=1000,\n",
    "                seed=10)\n",
    "hp_se\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c65fbf",
   "metadata": {},
   "source": [
    "This indicates that the bootstrap estimate for ${\\rm SE}(\\hat{\\beta}_0)$ is\n",
    "0.85, and that the bootstrap\n",
    "estimate for ${\\rm SE}(\\hat{\\beta}_1)$ is\n",
    "0.0074.  As discussed in\n",
    "Section 3.1.2, standard formulas can be used to compute\n",
    "the standard errors for the regression coefficients in a linear\n",
    "model. These can be obtained using the `summarize()`  function\n",
    "from `ISLP.sm`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c9aea297",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept     0.717\n",
       "horsepower    0.006\n",
       "Name: std err, dtype: float64"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp_model.fit(Auto, Auto['mpg'])\n",
    "model_se = summarize(hp_model.results_)['std err']\n",
    "model_se\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d870ad6b",
   "metadata": {},
   "source": [
    "The standard error estimates for $\\hat{\\beta}_0$ and $\\hat{\\beta}_1$\n",
    "obtained using the formulas  from Section 3.1.2  are\n",
    "0.717 for the\n",
    "intercept and\n",
    "0.006 for the\n",
    "slope. Interestingly, these are somewhat different from the estimates\n",
    "obtained using the bootstrap.  Does this indicate a problem with the\n",
    "bootstrap? In fact, it suggests the opposite.  Recall that the\n",
    "standard formulas given in\n",
    " {Equation 3.8 on page 82}\n",
    "rely on certain assumptions. For example,\n",
    "they depend on the unknown parameter $\\sigma^2$, the noise\n",
    "variance. We then estimate $\\sigma^2$ using the RSS. Now although the\n",
    "formula for the standard errors do not rely on the linear model being\n",
    "correct, the estimate for $\\sigma^2$ does.  We see\n",
    " {in Figure 3.8 on page 108}  that there is\n",
    "a non-linear relationship in the data, and so the residuals from a\n",
    "linear fit will be inflated, and so will $\\hat{\\sigma}^2$.  Secondly,\n",
    "the standard formulas assume (somewhat unrealistically) that the $x_i$\n",
    "are fixed, and all the variability comes from the variation in the\n",
    "errors $\\epsilon_i$.  The bootstrap approach does not rely on any of\n",
    "these assumptions, and so it is likely giving a more accurate estimate\n",
    "of the standard errors of $\\hat{\\beta}_0$ and $\\hat{\\beta}_1$ than\n",
    "the results from `sm.OLS`.\n",
    "\n",
    "Below we compute the bootstrap standard error estimates and the\n",
    "standard linear regression estimates that result from fitting the\n",
    "quadratic model to the data. Since this model provides a good fit to\n",
    "the data (Figure 3.8), there is now a better\n",
    "correspondence between the bootstrap estimates and the standard\n",
    "estimates of ${\\rm SE}(\\hat{\\beta}_0)$, ${\\rm SE}(\\hat{\\beta}_1)$ and\n",
    "${\\rm SE}(\\hat{\\beta}_2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "79c56529",
   "metadata": {
    "execution": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept                                  2.067840\n",
       "poly(horsepower, degree=2, raw=True)[0]    0.033019\n",
       "poly(horsepower, degree=2, raw=True)[1]    0.000120\n",
       "dtype: float64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quad_model = MS([poly('horsepower', 2, raw=True)])\n",
    "quad_func = partial(boot_OLS,\n",
    "                    quad_model,\n",
    "                    'mpg')\n",
    "boot_SE(quad_func, Auto, B=1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fccbbbd",
   "metadata": {},
   "source": [
    "We  compare the results to the standard errors computed using `sm.OLS()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4d0b4edc",
   "metadata": {
    "execution": {},
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept                                  1.800\n",
       "poly(horsepower, degree=2, raw=True)[0]    0.031\n",
       "poly(horsepower, degree=2, raw=True)[1]    0.000\n",
       "Name: std err, dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = sm.OLS(Auto['mpg'],\n",
    "           quad_model.fit_transform(Auto))\n",
    "summarize(M.fit())['std err']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a86ff6e",
   "metadata": {},
   "source": [
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "formats": "Rmd,ipynb",
   "main_language": "python"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}