The Goal

We have been learning a lot of code in our quest to explore cross validation techniques. Today, we are going to dig more deeply into the code behind 10-fold CV.

The Data

We are going to stick with our penguins data set from last lab. To load the data, put the following lines of code inside a chunk and press play:

library(palmerpenguins)
data("penguins")
penguins <- na.omit(penguins)
# Remove the feature year: we don't want to use it today
penguins <- data.frame(penguins)[,-8]

Our goal for today is to predict \(Y\) = the body mass of the penguin. As we are provided with only one data set, we will need to determine a way to assess the predictive accuracy of any technique we choose. We will do this today using 10-fold CV. 10-fold CV involves dividing our original data set into 10 different subsets called folds. Each fold will contain roughly the same number of rows, and each row from the original data set will be assigned to exactly one fold.

In class, we discussed the folds in terms of tables at a restaurant. We have 10 folds (so 10 tables) and now we need to think about how many penguins need to be assigned to each fold (table).

Once we have determined how many rows (penguins) need to go in to each fold (table), we need to complete the process of actually assigning each row to a fold (penguin to a table). This involves sampling without replacement. Before we try this on the real data set, let’s try a smaller example.

Small Data Example

Let’s create a smaller data set to work with to illustrate the code.

tickets <- rep(1:5, 4)

tickets

tickets_Q4 <- 

# Set the random seed
set.seed(363663)

# Draw tickets from the box
foldsSmall <- sample(tickets, 20, replace = FALSE)

which(foldsSmall == 1)

fold1 <- penguinsSmall[which(foldsSmall==1), ]

10-fold CV: Running the Loop

Now that we know how to create the folds, let’s move on to the process of running the loop. We need to:

1. Set f = 1
2. Set fold f as validation data and the rest of the folds as training data
3. Build a linear regression model to predict body mass using all other columns as features.
4. Store the predictions for the rows in the validation data.
5. And repeat!

This should sound like a for loop, because it is one! We also need one step before we run the loop, which is creating our storage space:

# Create storage space 
predictions <- data.frame("body_mass_g" = rep(NA,n))

One of the tricks to making a for loop is to write the inside first. In other words, write the code that would give the results for f = 1. Once we have done that, it is easier to figure out what parts of the code need to change with f and create the loop. Let’s try it.

Building the inside first makes it easier to catch any errors. This in essence shows us what each iteration of the loop will do. Once it does what we want, we can write the loop!

compute_RMSE <- function(truth, predictions){

  # Part 1
  part1 <- truth - predictions
  
  #Part 2
  part2 <- part1^2
  
  #Part 3: RSS 
  part3 <- sum(part2)
  
  #Part4: MSE
  part4 <- 1/length(predictions)*part3
  
  # Part 5:RMSE
  sqrt(part4)
}

KNN

Now that we have seen how the process words, let’s try another technique to predict body mass: KNN.

Look at Row 1 in the penguins data and Row 53 in the penguins data.

We want to use the Gower’s Distance to find the distance between these two rows. Since we have two variables that are categorical with more than two levels, our code will actually “see” a slightly different version of the two rows.

We are going to use 10-fold CV to estimate the test RMSE of 5-NN (5-nearest neighbors). Before running KNN, we need to load the StatMatch library and teach R the function to compute Gower’s Distance.

suppressMessages(library(StatMatch))

knnGower <- function(trainX, testX, trainY, K){
  # Convert the data 
  trainX <- model.matrix(trainY~ ., data = trainX)[,-1]
  holder <- 1:nrow(testX)
  testX <- model.matrix(holder~ ., data = testX)[,-1]
  
  # Find the Gower Distance
  gowerHold <- StatMatch::gower.dist( testX, trainX)
  # For each row, find the k smallest
  neighbors <- apply(gowerHold, 1, function(x) which(x %in% sort(x)[1:K]))
  
  if(class(neighbors)[1]=="integer"){
    preds <- trainY[neighbors]
  }
  
  # Take the mean to get the prediction
  if(class(neighbors)[1]=="matrix"){
    preds <- apply(neighbors, 2,function(x) mean(trainY[x]))
  }
  
  if(class(neighbors)[1]=="list"){
    preds <- lapply(neighbors, function(x) mean(trainY[x]))
    preds <- unlist(preds)
  }
  # Return the predictions
  unlist(preds)
}

Recall that to use this function, the structure is:

knnGower( trainingFeatures, testFeatures, trainingResponse, K = numberOfNeighbors)

We want to use 10-fold CV to assess the predictive accuracy of this approach. Our steps for the 10-fold CV with KNN are:

1. Set f = 1
2. Set fold f as validation data and the rest of the folds as training data
3. Use 5-NN to make predictions for the rows in the validation data.
4. Store the predictions
5. And repeat!

Considering K

In the work we have done so far, we have used 5-NN. However, is that the best choice for the number of neighbors? Could we end up with a smaller validation RMSE if we choose a different number of neighbors \(K\)?

To determine this, we actually need to try out different options of \(K\). This means not only do we need a for loop that loops over the folds, we also need one that allows us to explore different choices of \(K\). The code structure we need for this is called a Nested for-loop. You will work with this for your Project 1.