When working on data science or research teams, it often helps to have a workflow that makes it easy for teammates to review your work and add additional components to the data cleaning or model structure. Additionally, it ensures that the steps in the process are clear so that debugging is easy.

In R, the {tidymodels} package offers one such workflow and in python, folks seem to prefer scikit learn. This week, I’m going to walk through a full workflow in tidymodels using the K-nearest neighbor algorithm. Some of the things I’ll cover include:

• Splitting the data in test/train sets and cross validation folds
• Setting up the model structure (what tidymodels refers to as a recipe)
• Creating preprocessing steps
• Compiling the preprocessing and model structure together into a single workflow
• Tuning the KNN model
• Identifying the best tuned model
• Evaluating the model on the test set
• Saving the entire worflow and model to be used later on with new data
• Making predictions with the saved workflow and model on new data

I’ve provided a number of tidymodels tutorials on this blog so feel free to search the blog for those. Additionally, all of my tidymodels tutorials and templates are available in a GitHub Repo. Alternatively, if Python is more you jam, I did do a previous blog comparing the workflow in tidymodels to Scikit Learn, which can be found HERE.

This entire tutorial and data are available on my GitHub page if you’d like to code along.

Load & Clean Data

For this tutorial I’m going to use the 2021 FIFA Soccer Ratings, which are available on Kaggle. There are several years of ratings there, but we will concentrate on 2021 with the goal being to estimate a player’s contract value based on the various FIFA ratings provided and the position they play.

First, I load the tidyverse and tidymodels libraries and read in the data. I’m going to drop goalies so that we only focus on players who play the field and, since the data has a large number of columns, I’m going to get only the columns we care about: Information about the player (playerID, name, height, weight, club, league, position, and contract value) and then the various FIFA ratings.

The only column with missing data that impacts our analysis is the team position column, as this is a feature in the data set. Additionally, there are 207 players with a contract value of 0 (all 195 players missing a team position have a 0 contract value). So, perhaps these are players who weren’t on a club at the time the ratings were assigned.

I’m going to remove these players from our analysis data set, but I’m going to place them in their own data set because we will use them later to make estimates of their contract value.

Visualizing the Data

Let’s take a look at a visual of the contract value for all of the players in our data set. Since this is severely right skewed, we will plot it on the log scale.

We also notice there are a variety of positions (way more than I would have guessed in soccer!).

Finally, we want to explore how playing position might influence contract value.

tidymodels set up

Now that we have our data organized we want to set up a tidymodels workflow to estimate the contract value of every player using a K-nearest neighbor model.

First, we split the data in train and test splits and then further split the training data into 5 cross validation folds so that we can tune the model to find the best number of K-neighbors.

Next, we specify that we want to use a KNN model and the outcome variable as regression, since it is continuous (as opposed to classification).

Now that the model is specified we need to set up our preprocessing steps. This is one thing that tidymodels is super useful for. Normally, we’d need to do all this preprocessing before fitting the model and, as in the case of KNN and most machine learning models, we’d need to standardize the variables in the training set and then store the means and standard deviations of those variables so that we can use them to standardize the test set or future data. In tidymodels, we set up our preprocessing workflow and store it and it contains all that information for us! We can simply load it and use it for any downstream analysis we want to do. There are three preprocessing steps we want to add to our workflow:

1. log transform the dependent variable (contract value)
2. Dummy code the positions
3. Standardize (z-score) all off the numeric variable in the model (the ratings)

We can view the preprocessing steps using the prep() function.

 

With the model specified and the preprocessing recipe set up we are ready to compile everything into a single workflow.

We can look at the workflow to make sure everything makes sense before we start fitting the model.

Next, we tune the model using the cross validated folds that we set up on our training data set. We are trying to find the optimal number of k-neighbors that minimizes the RMSE of the outcome variable (the log of contract value).

We see the results are stored in a list for every one of our 5 cross validation splits. We can view the model metrics for each of the folds.

The model with 10 neighbors appears to have the lowest RMSE. So, we want to pull that value directly from this table so that we can store it and use it later when we go to fit the final model.

Fitting the Final Model Version 1: Using finalize_workflow()

Version 1 that I’ll show for fitting the final model uses the finalize_workflow() function. I like this approach if I’m building my model in a local session and then want to see how well it performs on the test session (which is what this function does). This isn’t the approach I take when I need to save everything for downstream use (which we will cover in Version 2).

First, fit the best model using the finalize_workflow() function.

Then, we get the model predictions for this final workflow on the the test data set.

We can plot the predicted contract value compared to the actual contract value and calculate the test set RMSE.

Fitting the Final Model Version 2: Re-sepcify the model on the training set, save it, and use it later

In this second version of fitting the final model, we will take the best neighbors from the tuning phase, re-specify the model to the training set, reset the workflow, and then save these components so that we can use them later.

First, re-specify the model and notice I set the neighbors argument to best_neighbors, which we saved above when we tuned the model.

We then use this finalized/re-specified model to reset the workflow. Notice it is added to the add_model() function however the recipe, with our preprocessing steps, does not change.

With the workflow set up we now need to do the final two sets:

1. Fit the model to the entire data set and extract the recipe using extract_recipe() since we will need to save this to preprocess any new data before making predictions.
2. Fit the model to the final data and extract the model itself, using extract_fit_parsnip() so we can use the model to make future predictions.

Now we save the recipe and model as .rda files. We can then load these two structures and use them on new data!

We have 207 players with 0 contract value, of which 195 had no team position. So, that leaves us with 12 players that we can estimate contract value for with our saved model.

First, we use the bake() function to apply the saved recipe to the new data so that we can preprocess everything appropriately.

With the new data preprocessed we can now predict the log of contract value with our saved model.

We can add these predictions back to the new data set, exponentiate the predicted contract value to get it back to the normal scale and then create a visual of our estimates for the 12 players.

Wrapping Up

That’s a quick walk through on how to set up a complete tdymodels workflow. This approach would work for any type of model you want to build, not just KNN! As always, the code and data are available on my GitHub page.

For the next installment of our simulation blog series we will use simulation to look at the Multicollinearity assumption in linear regression.

For a refresher, here are the previous 5 articles in this series:

• Part 1 discussed the basic functions for simulating and sampling data in R.
• Part 2 walked us through how to perform bootstrap resampling and then simulate bivariate and multivariate distributions.
• Part 3 we worked making group comparisons by simulating thousands of t-tests
• Part 4 building simulations for linear regression
• Part 5 using simulation to investigate the homoskedasticity assumption in regression.

The entire code for this series is accessible on my GITHUB page.

Multicollinearity

Multicollinearity occurs when two independent variables in a regression model are highly correlated with each other. Such a situation can produce problems with interpretation of the beta coefficients of the model, may increase standard errors in the model, and can lead to over fitting of the data. We can simulate this issue in order to get a better understanding of how multicollinearity can influence a regression model.

Constructing the simulation

We will use the mvnorm package to help us construct a simulation where the two independent variables share a certain level of correlation between each other.

First, we will create the true model parameters: an intercept of 2, a beta1 of 5, and a beta2 of 10. We also create a vector of correlation coefficients from 0 to 0.99 and a few data frames to store the results of our model. We will also specify that at each correlation coefficient we want 200 random draws from the multivariate normal distribution.

## load packages
library(tidymodels)
library(patchwork)
library(mvtnorm)

set.seed(999)

# create the true model parameters
intercept <- 2
beta1 <- 5
beta2 <- 10

## Number of draws from a multivariate normal distribution
n <- 200

## Create a data frame to store model results
sim_params <- data.frame(intercept = NA,
                      intercept_se = NA,
                      beta1 = NA,
                      beta1_se = NA,
                      beta2 = NA,
                      beta2_se = NA,
                      model_rse = NA)

## create levels of multicollinearity between the two independent variables
cor_coefs <- c(seq(from = 0, to = 0.9, by = 0.1), 0.99)

# data frame to store the average beta coefficient and their standard deviations form the simulation
mean_betas <- data.frame(beta1 = NA,
                       sd_beta1 = NA,
                       beta2 = NA,
                       sd_beta2 = NA)

Next, we will create a nested for() loop to construct out simulations.

•  The outer part of the loop begins by creating the multivariate normal distribution. We use the rmvnorm(), which means we first specify a correlation matrix using our vector of correlation coefficients we specified above. Once we have the two correlated variables we can put them into the inner loop.
• The inner loop is where we create a regression equation for the given correlation between the two variables. We create 100 regression simulations for each correlation coefficient.
• Once the inner loop is finished running we store the results at the bottom of the outer loop. Instead of storing all of the results, we take the average of the 100 beta coefficients and their respective standard errors for each correlation coefficient.

## loop
for(j in 1:length(cor_coefs)){
  
  ## Create a correlation matrix between beta1 and beta2
  beta_corr <- matrix(c(1, cor_coefs[j], cor_coefs[j], 1), nrow = 2, ncol = 2)
  
  ## create a multivariate normal distribution 
  cor_df <- rmvnorm(n = n, mean = c(0, 0), sigma = beta_corr)
  X1 <- cor_df[, 1]
  X2 <- cor_df[, 2]
  
  ## simulate 100 regression simulations
  for(i in 1:100){
    
    # set up the model
    y_hat <- intercept + beta1*X1 + beta2*X2 + rnorm(n = n, mean = 0, sd = 1)
    
    # construct a regression equation
    model <- lm(y_hat ~ X1 + X2)
    
    # store the variance-covariance matrix
    vcv <- vcov(model)
    
    # estimates for the intercept
    sim_params[i, 1] <- model$coef[1]
  
    # estimates for the beta1
    sim_params[i, 3] <- model$coef[2]
    
    # estimates for beta2
    sim_params[i, 5] <- model$coef[3]
  
    # SE for the intercept
    sim_params[i, 2] <- sqrt(diag(vcv)[1])
    
    # SE for beta1
    sim_params[i, 4] <- sqrt(diag(vcv)[2])
    
    # SE for beta2
    sim_params[i, 6] <- sqrt(diag(vcv)[3])
  
    # model RSE
    sim_params[i, 7] <- summary(model)$sigma
  }

  mean_betas[j, ] <- c(mean(sim_params[, 3]), mean(sim_params[, 4]), mean(sim_params[, 5]), mean(sim_params[, 6]))

}

The final result is a data frame with 10 rows, one for each correlation coefficient and the average of 100 regression simulations for the beta coefficients and their standard errors.

# Add the correlation coefficients to the final data frame
mean_betas$cor_coef <- cor_coefs mean_betas %>%
  knitr::kable()

Finally, we plot each of the results with the correlation coefficients on the x-axis.

# Plot the results
par(mfrow = c(2,2))
plot(x = mean_betas$cor_coef,
     y = mean_betas$beta1,
     main = "Beta1",
     lwd = 3,
     type = "b",
     ylim = c(2, 7),
     xlab = "Correlation betwen X1 and X2",
     ylab = "Beta1")
plot(x = mean_betas$cor_coef,
     y = mean_betas$sd_beta1,
     main = "Beta1 Standard Error",
     lwd = 3,
     type = "b",
     xlab = "Correlation betwen X1 and X2",
     ylab = "SE of B1")
plot(x = mean_betas$cor_coef,
     y = mean_betas$beta2,
     main = "Beta2",
     lwd = 3,
     type = "b",
     ylim = c(7, 13),
     xlab = "Correlation betwen X1 and X2",
     ylab = "Beta 2")
plot(x = mean_betas$cor_coef,
     y = mean_betas$sd_beta2,
     main = "Beta2 Standard Error",
     lwd = 3,
     type = "b",
     xlab = "Correlation betwen X1 and X2",
     ylab = "SE of B2")

The beta coefficients themselves remain relatively unchanged in our simulation across the various correlations levels. However, once the correlation between the two independent variables reaches about 0.7 the standard errors around the beta coefficients begin to increase exponentially, increasing our uncertainty about the true parameter values.

Wrapping Up

We’ve now covered using simulation to investigate two assumptions of linear regression. Our next installment will investigate another linear regression assumption before we proceed on to simulating other types of models.

For the full code, check out my GITHUB page.