Ch05-ISLR-resample-lab.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "id": "b4fa2bec",
   "metadata": {},
   "source": [
    "---\n",
    "title: \"Ch5-resample-lab\"\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67f91564",
   "metadata": {},
   "source": [
    "\n",
    "# Lab: Cross-Validation and the Bootstrap\n",
    "\n",
    "In this lab, we explore the resampling techniques covered in this chapter. Some of the commands in this lab may take a while to run on your computer.\n",
    "\n",
    "## The Validation Set Approach\n",
    "\n",
    "We explore the use of the validation set approach in order to estimate the test error rates that result from fitting various linear models on the `Auto` data set.\n",
    "\n",
    "Before we begin, we use the `set.seed()` function in order to set a  for `R`'s random number generator, so that the reader of this book will\n",
    "obtain precisely the same results as those shown below. It is generally\n",
    "a good idea to set a random seed when performing an analysis such as cross-validation that contains an\n",
    "element of randomness, so that the results obtained can be reproduced precisely at a later time.\n",
    "\n",
    "We begin by using the\n",
    "`sample()` function to split the set of observations into two halves, by selecting a random subset of $196$ observations out of the original $392$ observations. We  refer\n",
    "to these observations as the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eb07b78",
   "metadata": {
    "name": "chunk1"
   },
   "outputs": [],
   "source": [
    "library(ISLR2)\n",
    "set.seed(1)\n",
    "train <- sample(392, 196)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ead673",
   "metadata": {},
   "source": [
    "(Here we use a shortcut in the sample command; see `?sample` for details.)\n",
    "We then use the `subset` option in `lm()` to fit a linear regression using only the observations corresponding to the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eb5f80aa",
   "metadata": {
    "name": "chunk2"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>1</li><li>3</li><li>2</li><li>5</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1\n",
       "\\item 3\n",
       "\\item 2\n",
       "\\item 5\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1\n",
       "2. 3\n",
       "3. 2\n",
       "4. 5\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 1 3 2 5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lm.fit <- lm(mpg ~ horsepower, data = Auto, subset = train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ae38b1",
   "metadata": {},
   "source": [
    "We now  use\n",
    " the `predict()` function to estimate the response for all $392$ observations,  and\n",
    " we  use\n",
    "  the `mean()` function to calculate the MSE of the $196$ observations in the validation set. Note that the `-train` index below selects  only the observations that are not in the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0bc8b67",
   "metadata": {
    "name": "chunk3"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>1</li><li>6</li><li>2</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1\n",
       "\\item 6\n",
       "\\item 2\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1\n",
       "2. 6\n",
       "3. 2\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 1 6 2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "attach(Auto)\n",
    "mean((mpg - predict(lm.fit, Auto))[-train]^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0cd124e",
   "metadata": {},
   "source": [
    "Therefore, the estimated test MSE for the linear regression fit is $23.27$. We can use the `poly()` function to estimate the test error for the quadratic and cubic regressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "39555396",
   "metadata": {
    "name": "chunk4"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "3"
      ],
      "text/latex": [
       "3"
      ],
      "text/markdown": [
       "3"
      ],
      "text/plain": [
       "[1] 3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "3"
      ],
      "text/latex": [
       "3"
      ],
      "text/markdown": [
       "3"
      ],
      "text/plain": [
       "[1] 3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>2</li><li>10</li><li>5</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 2\n",
       "\\item 10\n",
       "\\item 5\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 2\n",
       "2. 10\n",
       "3. 5\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]  2 10  5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto, \n",
    "    subset = train)\n",
    "mean((mpg - predict(lm.fit2, Auto))[-train]^2)\n",
    "lm.fit3 <- lm(mpg ~ poly(horsepower, 3), data = Auto, \n",
    "    subset = train)\n",
    "mean((mpg - predict(lm.fit3, Auto))[-train]^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d0549e0",
   "metadata": {},
   "source": [
    "These error rates are $18.72$ and $18.79$, respectively.\n",
    "If we choose a different training set instead, then we will obtain somewhat different errors on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "259e0074",
   "metadata": {
    "name": "chunk5"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'x'</li><li>'y'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'x'\n",
       "\\item 'y'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'x'\n",
       "2. 'y'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"x\" \"y\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [],
      "text/latex": [],
      "text/markdown": [],
      "text/plain": [
       "character(0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "set.seed(2)\n",
    "train <- sample(392, 196)\n",
    "lm.fit <- lm(mpg ~ horsepower, subset = train)\n",
    "mean((mpg - predict(lm.fit, Auto))[-train]^2)\n",
    "lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto, \n",
    "    subset = train)\n",
    "mean((mpg - predict(lm.fit2, Auto))[-train]^2)\n",
    "lm.fit3 <- lm(mpg ~ poly(horsepower, 3), data = Auto, \n",
    "    subset = train)\n",
    "mean((mpg - predict(lm.fit3, Auto))[-train]^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1b1497c",
   "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 $25.73$, $20.43$, and $20.39$, respectively.\n",
    "\n",
    "These results are consistent with our previous findings: a model that predicts `mpg` using a quadratic function of `horsepower` performs better than a model that involves only a linear function of `horsepower`, and there is little evidence in favor of a model that uses a cubic function of `horsepower`.\n",
    "\n",
    "## Leave-One-Out Cross-Validation\n",
    "\n",
    "The LOOCV estimate can be automatically computed for any generalized linear model   using the `glm()` and `cv.glm()` functions.  In the lab for Chapter 4, we used the `glm()` function to perform logistic regression by passing in  the `family = \"binomial\"` argument.\n",
    " But if we use `glm()` to fit a model without passing in the `family` argument, then it  performs linear regression, just like the `lm()` function.\n",
    "So for instance,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "99306afa",
   "metadata": {
    "name": "chunk6"
   },
   "outputs": [],
   "source": [
    "glm.fit <- glm(mpg ~ horsepower, data = Auto)\n",
    "coef(glm.fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87c22e0e",
   "metadata": {},
   "source": [
    "and"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e6035c0",
   "metadata": {
    "name": "chunk7"
   },
   "outputs": [],
   "source": [
    "lm.fit <- lm(mpg ~ horsepower, data = Auto)\n",
    "coef(lm.fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf3f7a7",
   "metadata": {},
   "source": [
    " yield identical linear regression models. In this lab, we will  perform linear regression using\n",
    " the `glm()` function rather than the `lm()` function because the former can be used together with\n",
    "`cv.glm()`. The `cv.glm()` function is part of the `boot` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8584a3d3",
   "metadata": {
    "name": "chunk8"
   },
   "outputs": [],
   "source": [
    "library(boot)\n",
    "glm.fit <- glm(mpg ~ horsepower, data = Auto)\n",
    "cv.err <- cv.glm(Auto, glm.fit)\n",
    "cv.err$delta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f64a4086",
   "metadata": {},
   "source": [
    "The `cv.glm()` function produces a list with several components.  The two numbers in the `delta` vector contain the cross-validation results. In this case the numbers are identical (up to two decimal places) and correspond to the LOOCV statistic given in ( 5.1). Below, we discuss a situation in which the two numbers differ. Our cross-validation estimate for the test error is approximately $24.23$.\n",
    "\n",
    "We can repeat this procedure for increasingly complex polynomial fits.\n",
    " To automate the process, we use the  `for()` function to initiate a  which iteratively fits polynomial regressions for polynomials of order $i=1$ to $i=10$, computes the associated cross-validation error, and stores it in the $i$th element of the vector `cv.error`.\n",
    " We begin by initializing the vector. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f353e50",
   "metadata": {
    "name": "chunk9"
   },
   "outputs": [],
   "source": [
    "cv.error <- rep(0, 10)\n",
    "for (i in 1:10) {\n",
    "  glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)\n",
    "  cv.error[i] <- cv.glm(Auto, glm.fit)$delta[1]\n",
    "}\n",
    "cv.error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67634a57",
   "metadata": {},
   "source": [
    "As in Figure 5.4, we see a sharp drop in the estimated test MSE between the linear and quadratic fits, but then no clear improvement from using higher-order polynomials.\n",
    "\n",
    "## $k$-Fold Cross-Validation\n",
    "\n",
    "The `cv.glm()` function can also be used to implement $k$-fold CV. Below we use $k=10$, a common choice for $k$, on the `Auto` data set.\n",
    "We once again set a random seed and initialize a vector in which we will store the CV errors corresponding to the polynomial fits of orders one to ten."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12806f8a",
   "metadata": {
    "name": "chunk10"
   },
   "outputs": [],
   "source": [
    "set.seed(17)\n",
    "cv.error.10 <- rep(0, 10)\n",
    "for (i in 1:10) {\n",
    "  glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)\n",
    "  cv.error.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]\n",
    "}\n",
    "cv.error.10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c6887a5",
   "metadata": {},
   "source": [
    "Notice that the computation time is shorter than that of LOOCV.\n",
    "(In principle, the computation time for LOOCV for a least squares linear model should be faster than for $k$-fold CV, due to the availability\n",
    "of the formula ( 5.2) for LOOCV; however, unfortunately the `cv.glm()` function does not make use of this formula.)\n",
    "We still see little evidence that using cubic or higher-order polynomial terms leads to lower test error than simply using a quadratic fit.\n",
    "\n",
    "We saw in Section 5.3.2 that the two numbers associated with `delta` are essentially the same when LOOCV is performed.\n",
    "When we instead perform $k$-fold CV, then the two numbers associated with `delta` differ slightly. The first is the standard $k$-fold CV estimate,\n",
    "as in ( 5.3). The second is a bias-corrected version. On this data set, the two estimates are very similar to each other.\n",
    "\n",
    "## The Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6cbb6ad",
   "metadata": {},
   "source": [
    "We illustrate the use of the bootstrap in the simple example of Section 5.2, as well as on an example involving estimating the\n",
    "accuracy of the linear regression model on the `Auto` data set.\n",
    "\n",
    "\\subsubsection{Estimating the Accuracy of a Statistic of Interest}\n",
    "\n",
    "One of the great advantages of the bootstrap approach is that it can be\n",
    "applied in almost all situations. No complicated mathematical calculations\n",
    "are required. Performing a bootstrap analysis in `R` entails only two\n",
    "steps. First, we must create a function that computes the statistic of\n",
    "interest. Second, we use the `boot()` function, which is part of the `boot` library, to perform the bootstrap by repeatedly\n",
    "sampling observations from the data set with replacement.\n",
    "\n",
    "The `Portfolio` data set in the `ISLR2` package is simulated data of $100$ pairs of returns, generated in the fashion described in Section 5.2.\n",
    "To illustrate the use of the bootstrap on this data, we must first\n",
    "create a function, `alpha.fn()`, which takes as input the $(X,Y)$ data\n",
    "as well as a vector indicating which observations should be used to\n",
    "estimate $\\alpha$. The function then outputs the estimate for $\\alpha$\n",
    "based on the selected observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "996fc134",
   "metadata": {
    "name": "chunk11"
   },
   "outputs": [],
   "source": [
    "alpha.fn <- function(data, index) {\n",
    "  X <- data$X[index]\n",
    "  Y <- data$Y[index]\n",
    "  (var(Y) - cov(X, Y)) / (var(X) + var(Y) - 2 * cov(X, Y))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76724d5e",
   "metadata": {},
   "source": [
    "This function *returns*, or outputs, an  estimate for $\\alpha$ based on applying ( 5.7) to the observations indexed by the argument `index`.\n",
    "For instance, the following command tells `R` to estimate $\\alpha$ using\n",
    "all $100$ observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06100bce",
   "metadata": {
    "name": "chunk12"
   },
   "outputs": [],
   "source": [
    "alpha.fn(Portfolio, 1:100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75cda089",
   "metadata": {},
   "source": [
    "The next command  uses the `sample()` function to randomly select\n",
    "$100$ observations from the range $1$ to $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": null,
   "id": "e9166f4d",
   "metadata": {
    "name": "chunk13"
   },
   "outputs": [],
   "source": [
    "set.seed(7)\n",
    "alpha.fn(Portfolio, sample(100, 100, replace = T))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33d913d4",
   "metadata": {},
   "source": [
    "We can implement a bootstrap analysis by performing this command many times, recording all of\n",
    "the corresponding estimates for $\\alpha$, and computing the resulting\n",
    "standard deviation.\n",
    "However, the `boot()` function automates this approach. Below we produce $R=1,000$ bootstrap estimates for $\\alpha$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4d57757",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk14"
   },
   "outputs": [],
   "source": [
    "boot(Portfolio, alpha.fn, R = 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58823115",
   "metadata": {},
   "source": [
    "The final output shows that using the original data, $\\hat{\\alpha}=0.5758$,\n",
    "and that the bootstrap estimate for ${\\rm SE}(\\hat{\\alpha})$ is $0.0897$.\n",
    "\n",
    "\\subsubsection{Estimating the Accuracy of a Linear Regression Model}\n",
    "\n",
    "The bootstrap approach can be used  to assess the\n",
    "variability of the coefficient estimates and predictions from a statistical learning method. Here we use the bootstrap approach in order to assess the variability of\n",
    "the estimates for $\\beta_0$ and $\\beta_1$, the intercept and slope terms for the linear regression model\n",
    "that uses  `horsepower` to predict `mpg` in the `Auto` data set. We will compare the estimates obtained using the bootstrap to those obtained using the formulas\n",
    "for ${\\rm SE}(\\hat{\\beta}_0)$ and ${\\rm SE}(\\hat{\\beta}_1)$ described\n",
    "in Section 3.1.2.\n",
    "\n",
    "We first create a simple function, `boot.fn()`, which takes in the\n",
    "`Auto` data set as well as a set of indices for the observations, and\n",
    "returns the intercept and slope estimates for the linear regression model. We then apply this function\n",
    "to the full set of $392$ observations in order to compute the estimates of $\\beta_0$ and $\\beta_1$ on the entire data set using the usual linear regression coefficient estimate\n",
    "formulas from Chapter 3. Note that we do not need the `{` and `}` at the beginning and end of the function because it is only one line long."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "556ab840",
   "metadata": {
    "name": "chunk15"
   },
   "outputs": [],
   "source": [
    "boot.fn <- function(data, index)\n",
    "  coef(lm(mpg ~ horsepower, data = data, subset = index))\n",
    "boot.fn(Auto, 1:392)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b2c04cf",
   "metadata": {},
   "source": [
    " The `boot.fn()` function can also be used in order to create\n",
    " bootstrap estimates for the intercept and slope terms by randomly sampling from among the observations with replacement. Here we give two examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f62400a",
   "metadata": {
    "name": "chunk16"
   },
   "outputs": [],
   "source": [
    "set.seed(1)\n",
    "boot.fn(Auto, sample(392, 392, replace = T))\n",
    "boot.fn(Auto, sample(392, 392, replace = T))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45b3bb50",
   "metadata": {},
   "source": [
    "Next, we use the `boot()` function to compute the standard errors of 1,000 bootstrap estimates for the intercept and slope terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "becb7122",
   "metadata": {
    "name": "chunk17"
   },
   "outputs": [],
   "source": [
    "boot(Auto, boot.fn, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8b27a40",
   "metadata": {},
   "source": [
    "This indicates that the bootstrap estimate for ${\\rm SE}(\\hat{\\beta}_0)$ is $0.84$, and that the bootstrap estimate for ${\\rm SE}(\\hat{\\beta}_1)$ is $0.0073$.\n",
    "As discussed in Section 3.1.2, standard formulas can be used to compute the standard errors for the regression coefficients in a linear model. These can be obtained using the  `summary()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cec25af1",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk18"
   },
   "outputs": [],
   "source": [
    "summary(lm(mpg ~ horsepower, data = Auto))$coef"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a4d5a4",
   "metadata": {},
   "source": [
    "The standard error estimates for $\\hat{\\beta}_0$ and\n",
    "$\\hat{\\beta}_1$ obtained using the formulas from\n",
    "Section 3.1.2 are $0.717$ for the intercept and $0.0064$\n",
    "for the slope. Interestingly, these are somewhat different from the\n",
    "estimates obtained using the bootstrap.  Does this indicate a problem\n",
    "with the bootstrap? In fact, it suggests the opposite.  Recall that\n",
    "the standard formulas given in Equation 3.8 on page66 rely on certain assumptions. For example, they depend\n",
    "on the unknown parameter $\\sigma^2$, the noise variance. We then estimate $\\sigma^2$\n",
    "using the RSS. Now although the formulas for the standard errors do not rely on the linear model\n",
    "being correct, the estimate for $\\sigma^2$ does.\n",
    "We see in\n",
    "Figure 3.8 on page91 that there is a non-linear relationship in\n",
    "the data, and so the residuals from a linear fit will be inflated, and so will $\\hat{\\sigma}^2$.\n",
    "Secondly, the standard formulas assume (somewhat unrealistically) that the $x_i$ are fixed, and all the variability comes from the variation in the errors $\\epsilon_i$.\n",
    " The bootstrap approach does not rely on any of these assumptions, and so it is\n",
    "likely giving a more accurate estimate of the standard errors of\n",
    "$\\hat{\\beta}_0$ and $\\hat{\\beta}_1$ than is the `summary()`\n",
    "function.\n",
    "\n",
    "Below we compute the bootstrap standard error estimates and the standard\n",
    "linear regression estimates that result from fitting the quadratic model to the data. Since this model provides a good fit to the data (Figure 3.8), there is now a better correspondence between the bootstrap estimates and the standard estimates of ${\\rm SE}(\\hat{\\beta}_0)$, ${\\rm SE}(\\hat{\\beta}_1)$ and ${\\rm SE}(\\hat{\\beta}_2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d5ca6e1",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "chunk19"
   },
   "outputs": [],
   "source": [
    "boot.fn <- function(data, index)\n",
    "  coef(\n",
    "      lm(mpg ~ horsepower + I(horsepower^2), \n",
    "        data = data, subset = index)\n",
    "    )\n",
    "set.seed(1)\n",
    "boot(Auto, boot.fn, 1000)\n",
    "summary(\n",
    "    lm(mpg ~ horsepower + I(horsepower^2), data = Auto)\n",
    "  )$coef"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cfca3e8",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "jupytext": {
   "formats": "ipynb,Rmd"
  },
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.5"
  },
  "substitutions": {
   "ISLPmod": "`ISLP`",
   "Rlang": "`R`",
   "mpl": "`matplotlib`",
   "numpy": "`numpy`",
   "pandas": "`pandas`",
   "pylang": "`python`",
   "smlib": "`statsmodels`"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}