Guiding Ideas

• How to use k-nearest neighbors for regression through the use of the knnreg() function from the caret package
• How to use decision trees for regression through the use of the rpart() function from the rpart package.
• How “making predictions” can be thought of as estimating the regression function, that is, the conditional mean of the response given values of the features.
• The difference between parametric and nonparametric methods.
• The difference between model parameters and tuning parameters methods.
• How these nonparametric methods deal with categorical variables and interactions.
• What is model flexibility?
• What is overfitting and how do we avoid it?

R Setup

library(tibble)     # data frame printing
library(dplyr)      # data manipulation

library(caret)      # fitting knn
library(rpart)      # fitting trees
library(rpart.plot) # plotting trees

library(knitr)      # creating tables
library(kableExtra) # styling tables

Mathematical Setup

Let’s return to the setup we defined in the previous lectures. Consider a random variable \(Y\) which represents a response variable, and \(p\) feature variables \(\boldsymbol{X} = (X_1, X_2, \ldots, X_p)\). We assume that the response variable \(Y\) is some function of the features, plus some random noise.

\[ Y = f(\boldsymbol{X}) + \epsilon \]

Our goal is to find some \(f\) such that \(f(\boldsymbol{X})\) is close to \(Y\). More specifically we want to minimize the risk under squared error loss.

\[ \mathbb{E}_{\boldsymbol{X}, Y} \left[ (Y - f(\boldsymbol{X})) ^ 2 \right] = \mathbb{E}_{\boldsymbol{X}} \mathbb{E}_{Y \mid \boldsymbol{X}} \left[ ( Y - f(\boldsymbol{X}) ) ^ 2 \mid \boldsymbol{X} = \boldsymbol{x} \right] \]

We saw previously (see the slides from last two content days) that this risk is minimized by the conditional mean of \(Y\) given \(\boldsymbol{X}\),

\[ \mu(\boldsymbol{x}) \triangleq \mathbb{E}[Y \mid \boldsymbol{X} = \boldsymbol{x}] \]

which we call the regression function.

Our goal then is to estimate this regression function. Let’s use an example where we know the true probability model:

\[ Y = 1 - 2x - 3x ^ 2 + 5x ^ 3 + \epsilon \]

where \(\epsilon \sim \text{N}(0, \sigma^2)\).

Recall that this implies that the regression function is

\[ \mu(x) = \mathbb{E}[Y \mid \boldsymbol{X} = \boldsymbol{x}] = 1 - 2x - 3x ^ 2 + 5x ^ 3 \]

Let’s also pretend that we do not actually know this information, but instead have some data, \((x_i, y_i)\) for \(i = 1, 2, \ldots, n\).

We simulate enough data to make the “pattern” clear-ish to recognize.

When we use a linear model, we first need to make an assumption about the form of the regression function.

For example, we could assume that

\[ \mu(x) = \mathbb{E}[Y \mid \boldsymbol{X} = \boldsymbol{x}] = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 \]

which is fit in R using the lm() function

lm(y ~ x + I(x ^ 2) + I(x ^ 3), data = sim_slr_data)

## 
## Call:
## lm(formula = y ~ x + I(x^2) + I(x^3), data = sim_slr_data)
## 
## Coefficients:
## (Intercept)            x       I(x^2)       I(x^3)  
##      0.8397      -2.7257      -2.3752       6.0906

Notice that what is returned are (maximum likelihood or least squares) estimates of the unknown \(\beta\) coefficients. That is, the “learning” that takes place with a linear models is “learning” the values of the coefficients.

For this reason, we call linear regression models parametric models. They have unknown model parameters, in this case the \(\beta\) coefficients that must be learned from the data. The form of the regression function is assumed.

What if we don’t want to make an assumption about the form of the regression function? While in this case, you might look at the plot and arrive at a reasonable guess of assuming a third order polynomial, what if it isn’t so clear? What if you have 100 features? Making strong assumptions might not work well.

Enter nonparametric models. We will consider two examples: k-nearest neighbors and decision trees.

k-Nearest Neighbors

We’ll start with k-nearest neighbors which is possibly a more intuitive procedure than linear models.¹

If our goal is to estimate the mean function,

\[ \mu(x) = \mathbb{E}[Y \mid \boldsymbol{X} = \boldsymbol{x}] \]

the most natural approach would be to use

\[ \text{average}(\{ y_i : x_i = x \}). \]

That is, to estimate the conditional mean at \(x\), average the \(y_i\) values for each data point where \(x_i = x\).

While this sounds nice, it has an obvious flaw. For most values of \(x\) there will not be any \(x_i\) in the data where \(x_i = x\)!

So what’s the next best thing? Pick values of \(x_i\) that are “close” to \(x\).

\[ \text{average}( \{ y_i : x_i \text{ equal to (or very close to) x} \} ). \]

This is the main idea behind many nonparametric approaches. The details often just amount to very specifically defining what “close” means.

In the case of k-nearest neighbors we use

\[ \hat{\mu}_k(x) = \frac{1}{k} \sum_{ \{i \ : \ x_i \in \mathcal{N}_k(x, \mathcal{D}) \} } y_i \]

as our estimate of the regression function at \(x\). While this looks complicated, it is actually very simple. Here, we are using an average of the \(y_i\) values of for the \(k\) nearest neighbors to \(x\).

The \(k\) “nearest” neighbors are the \(k\) data points \((x_i, y_i)\) that have \(x_i\) values that are nearest to \(x\). We can define “nearest” using any distance we like, but unless otherwise noted, we are referring to euclidean distance.² We are using the notation \(\{i \ : \ x_i \in \mathcal{N}_k(x, \mathcal{D}) \}\) to define the \(k\) observations that have \(x_i\) values that are nearest to the value \(x\) in a dataset \(\mathcal{D}\), in other words, the \(k\) nearest neighbors.

The plots below begin to illustrate this idea.

• In the left plot, to estimate the mean of \(Y\) at \(x = -0.5\) we use the three nearest neighbors, which are highlighted with green. Our estimate is the average of the \(y_i\) values of these three points indicated by the black x.
• In the middle plot, to estimate the mean of \(Y\) at \(x = 0\) we use the five nearest neighbors, which are highlighted with green. Our estimate is the average of the \(y_i\) values of these five points indicated by the black x.
• In the right plot, to estimate the mean of \(Y\) at \(x = 0.75\) we use the nine nearest neighbors, which are highlighted with green. Our estimate is the average of the \(y_i\) values of these nine points indicated by the black x.

You might begin to notice a bit of an issue here. We have to do a new calculation each time we want to estimate the regression function at a different value of \(x\)! For this reason, k-nearest neighbors is often said to be “fast to train” and “slow to predict.” Training, is instant. You just memorize the data! Prediction involves finding the distance between the \(x\) considered and all \(x_i\) in the data!³

So, how then, do we choose the value of the tuning parameter \(k\)? We validate!

First, let’s take a look at what happens with this data if we consider three different values of \(k\).

For each plot, the black dashed curve is the true mean function.

• In the left plot we use \(k = 25\). The red “curve” is the estimate of the mean function for each \(x\) shown in the plot.
• In the left plot we use \(k = 5\). The blue “curve” is the estimate of the mean function for each \(x\) shown in the plot.
• In the left plot we use \(k = 1\). The green “curve” is the estimate of the mean function for each \(x\) shown in the plot.

Some things to notice here:

• The left plot with \(k = 25\) is performing poorly. The estimated “curve” does not “move” enough. This is an example of an inflexible model.
• The right plot with \(k = 1\) might not perform too well. The estimated “curve” seems to “move” too much. (Notice, that it goes through each point. We’ve fit to the noise.) This is an example of a flexible model.

While the middle plot with \(k = 5\) is not “perfect” it seems to roughly capture the “motion” of the true regression function. We can begin to see that if we generated new data, this estimated regression function would perform better than the other two.

But remember, in practice, we won’t know the true regression function, so we will need to determine how our model performs using only the available data!

This \(k\), the number of neighbors, is an example of a tuning parameter. Instead of being learned from the data, like model parameters such as the \(\beta\) coefficients in linear regression, a tuning parameter tells us how to learn from data. It is user-specified. To determine the value of \(k\) that should be used, many models are fit to the estimation data, then evaluated on the validation. Using the information from the validation data, a value of \(k\) is chosen. (More on this in a bit.)

• Model parameters are “learned” using the same data that was used to fit the model.
• Tuning parameters are “chosen” using data not used to fit the model.

This tuning parameter \(k\) also defines the flexibility of the model and plays a similar role as complexity, or the number of parameters, did in our linear regressions. In KNN, a small value of \(k\) is a flexible model, while a large value of \(k\) is inflexible.⁴

Before moving to an example of tuning a KNN model, we will first introduce decision trees.

Decision Trees

Decision trees are similar to k-nearest neighbors but instead of looking for neighbors, decision trees create neighborhoods. We won’t explore the full details of trees, but just start to understand the basic concepts, as well as learn to fit them in R.

Neighborhoods are created via recursive binary partitions. In simpler terms, pick a feature and a possible cutoff value. Data that have a value less than the cutoff for the selected feature are in one neighborhood (the left) and data that have a value greater than the cutoff are in another (the right). Within these two neighborhoods, repeat this procedure until a stopping rule is satisfied. To make a prediction, check which neighborhood a new piece of data would belong to and predict the average of the \(y_i\) values of data in that neighborhood.

With the data above, which has a single feature \(x\), consider three possible cutoffs: -0.5, 0.0, and 0.75.

For each plot, the black vertical line defines the neighborhoods. The green horizontal lines are the average of the \(y_i\) values for the points in the left neighborhood. The red horizontal lines are the average of the \(y_i\) values for the points in the right neighborhood.

What makes a cutoff good? Large differences in the average \(y_i\) between the two neighborhoods. More formally we want to find a cutoff value that minimizes

\[ \sum_{i \in N_L} \left( y_i - \hat{\mu}_{N_L} \right) ^ 2 + \sum_{i \in N_R} \left(y_i - \hat{\mu}_{N_R} \right) ^ 2 \]

where

• \(N_L\) are the data in the left neighborhood, that is \(x < c\)
• \(N_R\) are the data in the right neighborhood, that is \(x > c\)
• \(\hat{\mu}_{N_L}\) is the mean of the \(y_i\) for data in the left neighborhood
• \(\hat{\mu}_{N_R}\) is the mean of the \(y_i\) for data in the right neighborhood

This quantity is the sum of two sum of squared errors, one for the left neighborhood, and one for the right neighborhood.

The table above summarizes the results of the three potential splits. We see that (of the splits considered, which are not exhaustive⁵) the split based on a cutoff of \(x = -0.50\) creates the best partitioning of the space.

Now let’s consider building a full tree.

In the plot above, the true regression function is the dashed black curve, and the solid orange curve is the estimated regression function using a decision tree. We see that there are two splits, which we can visualize as a tree.

The above “tree”⁶ shows the splits that were made. It informs us of the variable used, the cutoff value, and some summary of the resulting neighborhood. In “tree” terminology the resulting neighborhoods are “terminal nodes” of the tree. In contrast, “internal nodes” are neighborhoods that are created, but then further split.

The “root node” is the neighborhood contains all observations, before any splitting, and can be seen at the top of the image above. We see that this node represents 100% of the data. The other number, 0.21, is the mean of the response variable, in this case, \(y_i\).

Looking at a terminal node, for example the bottom left node, we see that 23% of the data is in this node. The average value of the \(y_i\) in this node is -1, which can be seen in the plot above.

We also see that the first split is based on the \(x\) variable, and a cutoff of \(x = -0.52\). Note that because there is only one variable here, all splits are based on \(x\), but in the future, we will have multiple features that can be split and neighborhoods will no longer be one-dimensional. However, this is hard to plot.

Let’s build a bigger, more flexible tree.

There are two tuning parameters at play here which we will call by their names in R which we will see soon:

• cp or the “complexity parameter” as it is called.⁷ This parameter determines which splits are accepted. A split must improve the performance of the tree by more than cp in order to be used. When we get to R, we will see that the default value is 0.1.
• minsplit, the minimum number of observations in a node (neighborhood) in order to consider splitting within a neighborhood.

There are actually many more possible tuning parameters for trees, possibly differing depending on who wrote the code you’re using. We will limit discussion to these two.⁸ Note that they effect each other, and they effect other parameters which we are not discussing. The main takeaway should be how they effect model flexibility.

First let’s look at what happens for a fixed minsplit by variable cp.

We see that as cp decreases, model flexibility increases. We see more splits, because the increase in performance needed to accept a split is smaller as cp is reduced.

Now the reverse, fix cp and vary minsplit.

We see that as minsplit decreases, model flexibility increases. By allowing splits of neighborhoods with fewer observations, we obtain more splits, which results in a more flexible model.

Example: Credit Card Data

# load data, coerce to tibble
crdt = as_tibble(ISLR::Credit)

# data prep
crdt = crdt %>%
  select(-ID) %>%
  select(-Rating, everything())

# test-train split
set.seed(1) 
crdt_trn_idx = sample(nrow(crdt), size = 0.8 * nrow(crdt))
crdt_trn = crdt[crdt_trn_idx, ]
crdt_tst = crdt[-crdt_trn_idx, ]

# check data
head(crdt_trn, n = 10)

## # A tibble: 10 × 11
##    Income Limit Cards   Age Education Gender   Student Married Ethnicity Balance
##     <dbl> <int> <int> <int>     <int> <fct>    <fct>   <fct>   <fct>       <int>
##  1  183.  13913     4    98        17 " Male"  No      Yes     Caucasian    1999
##  2   35.7  2880     2    35        15 " Male"  No      No      African …       0
##  3  123.   8376     2    89        17 " Male"  Yes     No      African …    1259
##  4   20.8  2672     1    70        18 "Female" No      No      African …       0
##  5   39.1  5565     4    48        18 "Female" No      Yes     Caucasian     772
##  6   36.5  3806     2    52        13 " Male"  No      No      African …     188
##  7   45.1  3762     3    80         8 " Male"  No      Yes     Caucasian      70
##  8   43.5  2906     4    69        11 " Male"  No      No      Caucasian       0
##  9   23.1  3476     2    50        15 "Female" No      No      Caucasian     209
## 10   53.2  4943     2    46        16 "Female" No      Yes     Asian         382
## # ℹ 1 more variable: Rating <int>

crdt_knn_01 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 1)
crdt_knn_10 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 10)
crdt_knn_25 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 25)
crdt_knn_100= knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 75)

head(crdt_knn_10$learn$X)

##   Age GenderFemale StudentYes
## 1  98            0          0
## 2  35            0          0
## 3  89            0          1
## 4  70            1          0
## 5  48            1          0
## 6  52            0          0

dist(head(crdt_knn_10$learn$X))

##           1         2         3         4         5
## 2 63.000000                                        
## 3  9.055385 54.009258                              
## 4 28.017851 35.014283 19.052559                    
## 5 50.009999 13.038405 41.024383 22.000000          
## 6 46.000000 17.000000 37.013511 18.027756  4.123106

sqrt(sum((crdt_knn_10$learn$X[3, ] - crdt_knn_10$learn$X[4, ]) ^ 2))

## [1] 19.05256

predict(crdt_knn_10, crdt_tst[1:5, ])

## [1] 334.9091 346.0000 459.3333 342.0909 450.8000

knn_mod_list = list(
  crdt_knn_01 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 1),
  crdt_knn_10 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 10),
  crdt_knn_25 = knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 25),
  crdt_knn_100= knnreg(Rating ~ Age + Gender + Student, data = crdt_trn, k = 75)
)

knn_tst_pred = lapply(knn_mod_list, predict, crdt_tst)

calc_rmse = function(actual, predicted) {
  sqrt(mean((actual - predicted) ^ 2))
}

sapply(knn_tst_pred, calc_rmse, crdt_tst$Rating)

##  crdt_knn_01  crdt_knn_10  crdt_knn_25 crdt_knn_100 
##     163.9868     142.2242     135.3838     142.6786

set.seed(5)
crdt_trn_idx = sample(nrow(crdt), size = 0.8 * nrow(crdt))
crdt_trn = crdt[crdt_trn_idx, ]
crdt_tst = crdt[-crdt_trn_idx, ]

crdt_tree = rpart(Rating ~ Age + Gender + Student, data = crdt_trn, cp = .0095)

crdt_tree

## n= 320 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 320 7045031.0 350.6250  
##    2) Age< 82.5 306 6586876.0 344.8170  
##      4) Age< 35.5 54  471770.1 311.8704 *
##      5) Age>=35.5 252 6043929.0 351.8770  
##       10) Age>=69.5 68 1571017.0 326.0882  
##         20) Age< 77.5 39  498793.0 287.3590 *
##         21) Age>=77.5 29  935056.1 378.1724  
##           42) Gender= Male 11  445504.2 316.2727 *
##           43) Gender=Female 18  421648.0 416.0000 *
##       11) Age< 69.5 184 4410974.0 361.4076 *
##    3) Age>=82.5 14  222217.4 477.5714 *

rpart.plot(crdt_tree)

crdt_tree_big = rpart(Rating ~ Age + Gender + Student, data = crdt_trn,
                      cp = 0.0, minsplit = 20)

rpart.plot(crdt_tree_big)

tree_mod_list = list(
  crdt_tree_0000 = rpart(Rating ~ Age + Gender + Student, data = crdt_trn, cp = 0.000),
  crdt_tree_0001 = rpart(Rating ~ Age + Gender + Student, data = crdt_trn, cp = 0.001),
  crdt_tree_0010 = rpart(Rating ~ Age + Gender + Student, data = crdt_trn, cp = 0.010),
  crdt_tree_0100 = rpart(Rating ~ Age + Gender + Student, data = crdt_trn, cp = 0.100)
)

tree_val_pred = lapply(tree_mod_list, predict, crdt_tst)

sapply(tree_val_pred, calc_rmse, crdt_tst$Rating)

## crdt_tree_0000 crdt_tree_0001 crdt_tree_0010 crdt_tree_0100 
##       182.0236       180.1045       172.6976       177.2816

crdt_tree_all = rpart(Rating ~ ., data = crdt_trn)

rpart.plot(crdt_tree_all)

calc_rmse(
  actual = crdt_tst$Rating,
  predicted = predict(crdt_tree_all, crdt_tst)
)

## [1] 45.17942

predict(crdt_tree_all, crdt_tst[1:5, ])

##        1        2        3        4        5 
## 520.7812 645.4000 520.7812 440.7368 520.7812