Task description

In this task, we are going to use various supervised learning methods to try to predict the antibody response to influenza vaccination using early gene expression measurements. The main goal of the practical is to understand the concepts behind cross-validation : why we need it, how it works, and how it is used to tune the hyperparameters of a model.

To save us some time, I have preprocessed the data for you. The data is an extract of open-access data from the Human Immune Project Consortium, with detail of the study found here. The relevant information for this exercise is as follows :

Antibodies were measured at days 0 and 28

Gene expression of around 10,000 genes was measured at days 0 and 7. To reduce the computation time, I kept only the 1,000 most variable genes

Flu is an interesting virus in this context, since it is a common disease, and many people may have already been infected with a given flu strain prior to vaccination. There is evidence in the literature that prior infection can modify the response to vaccination. Therefore, to take into account baseline (day 0) information, I have transformed both the response and the predictors into (log2) fold-change values. That is, the ratio between the post-vaccination and pre-vaccination measurements.

Our task will be to see if the day 7 fold-change of expression in some genes is predictive of the day 28 fold-change in antibody levels. That is, we seek to find an early gene expression signature which informs us of the later antibody response.

Setup

0a) Set up your project

Why do we need a project?

When you analyse data in R, you need to tell R where files are. For example, you might load data into your environment by defining a path "C:/Users/me/Documents/data/file.xlsx". However, this is problematic, since it is specific to your own computer setup. The code will then break as soon as you move the files, share them with a colleague, or re-run the analysis on a different machine.

To avoid this, we will introduce some basic organisational skills for an R project. Firstly, we will contain the work in a single R project folder (the root), and refer to every file relative to that folder. This makes specifying paths for loading data and saving outputs more robust.

Using a project folder also helps to organise files in a clean way.

raw data : the original data files, which you should never modify;

analysis : the code (such as this .Rmd file) that reads the raw data and produces results;

output : everything the code creates (figures, tables, processed data), which should be able to be perfectly regenerated by re-running the code.

Create the project folder

If you already created a project for a previous practical, you can reuse it : open it and go to “Add the files” below. Otherwise :

In RStudio, click File > New Project... > New Directory > New Project.

Give the project a name, for example omics_practicals, and choose where to save it on your computer. Click Create Project.

RStudio restarts inside your new project. You can check this in the top-right corner of RStudio, which now shows the project name, and in the Files pane, which shows a single file called e.g. omics_practicals.Rproj.

Next, create three sub-folders inside the project. In the Files pane of RStudio, click New Folder and create data-raw, analysis, and output. Alternatively, run the following code once in the R console.

dir.create("data-raw")
dir.create("analysis")
dir.create("output")

Your project should now look like this:

omics_practicals/
├── omics_practicals.Rproj    <- marks the "root" of the project
├── data-raw/                 <- the raw data files
├── analysis/                 <- the practicals (.Rmd files)
└── output/                   <- figures and results that you save

Add the files

Copy the data file I sent you, SDY1119_day7_GE_foldChange_top1000.RData, into the data-raw folder.

Copy this practical (the .Rmd file) into the analysis folder.

In the Files pane of RStudio, open the analysis folder and click on the .Rmd file to open it.

Important: always start your work by opening the project (double-click the .Rproj file, or use File > Open Project), and open the .Rmd file from within RStudio’s Files pane. If you open the .Rmd file directly from your Downloads folder, it will not be part of the project and the paths below will not work.

How do we refer to files?

By default, R Markdown runs code from the folder where the .Rmd file is saved (here, analysis/), which would make paths confusing. Instead, we will use the here package. The function here() always returns the path to the root of your project (the folder containing the .Rproj file), wherever the code is run from. We then build paths from this root. For example, here("data-raw", "my_file.RData") points to my_file.RData inside the data-raw folder, and here("output", "my_plot.pdf") points to a file in the output folder. The same code works on any computer and operating system, and you never have to type a full path.

0b) Install and load the required packages

If you ran the installation script I sent you, all the packages are already installed. Otherwise, the following chunk checks which of the packages needed for this practical you are missing and installs them : dplyr, ggplot2, caret, glmnet, pls, ranger and here.

cran_packages <- c(
  "dplyr",
  "ggplot2",
  "caret",
  "glmnet",
  "pls",
  "ranger",
  "here"
)

# Install missing CRAN packages
install.packages(setdiff(cran_packages, installed.packages()[, "Package"]))

Now load the packages.

0c) Check your project setup

Run the following chunk to check that everything is in the right place. The first line should print the path to your project folder (the one containing the .Rproj file), and the data file should be found.

# Path to the root of the project
here()

## [1] "/home/ah3/Desktop/Work/Teaching/PHDS_omics_data_course"

# Check that the data file exists
file.exists(here("data-raw", "SDY1119_day7_GE_foldChange_top1000.RData"))

## [1] TRUE

If the result is FALSE, or if the path printed above is not your project folder, go back to the previous steps : check that the file is in data-raw, and that you opened the project and then this .Rmd file from the Files pane.

0d) Load and explore the data

Load the data into the global environment, using here() to build the path to the file. The function load() reads the object saved in the file, and get() retrieves it.

df_all <- get(load(here("data-raw", "SDY1119_day7_GE_foldChange_top1000.RData")))

We can explore the data first to see what we’re dealing with. Let’s print the first few columns and rows.

df_all[, 1:8] %>%
  head()

##   participant_id    sex age response         a2m       abca1       abcb8
## 1 SUB179655.1119 Female  66      0.5  0.07524736 -0.54463090  0.41377154
## 2 SUB179657.1119   Male  75     -0.5  0.67610584  1.32433852 -0.26258410
## 3 SUB179659.1119   Male  70      1.0 -0.57017270 -3.46176587  0.09333269
## 4 SUB179662.1119   Male  65      1.0  0.84066094 -0.07048454  0.52547832
## 5 SUB179665.1119 Female  66      4.0  0.43971367 -1.55241659 -0.03439272
## 6 SUB179673.1119   Male  70      0.5  0.86480922 -2.32612497  0.05644670
##         abcc8
## 1  0.39287681
## 2 -0.31099800
## 3  0.15474707
## 4  0.04068224
## 5  0.51448642
## 6  0.58639101

We see that the first few columns give us some clinical information (participant id, sex, age), then we have the response value to be predicted, then the columns represent each gene, sorted alphabetically. Let’s plot the response as a histogram, and save the figure in the output folder.

response_hist <- ggplot(df_all, aes(x = response)) +
  geom_histogram(bins = 20) +
  labs(x = "Response (log2 fold-change in antibodies)", y = "Number of participants")

ggsave(here("output", "response_histogram.pdf"), response_hist, width = 8, height = 6)

response_hist

We have a distribution with a relatively long tail. Some participants have negative values. These are participants whose antibody levels decreased after vaccination! When we log transform a fold-change value between 0 and 1 (i.e. less antibodies at day 28 compared to day 0), we get a negative.

0e) Standardise the data

Many models are sensitive to the scale of the variables. Standardise the response and every gene, so that each has mean 0 and standard deviation 1, using the scale() function. Store the response as a numeric vector called y_vec, and the genes as a numeric matrix called X_mat (columns are genes and rows are participants; the genes are all the columns after response). Before proceeding, check explicitly that the column means are 0 and the standard deviations are 1 across your dataset.

# Your code here

For the rest of the tasks in this practical, you should use the standardised predictors and response as your inputs.

Now that we have finished with the preprocessing of the data, we can remove the unnecessary objects from the global environment with the rm() function. This is good practice with high-dimensional data, as memory can quickly be filled by large, unnecessary objects.

# List the objects to keep
objects_to_keep <- c("X_mat", "y_vec")

# Remove everything else from the global environment
rm(list = setdiff(ls(), objects_to_keep))

# Garbage collection to free unused memory
invisible(gc())

1) Building a predictive model

In this section, we are going to demonstrate how to fit a predictive model to a dataset and use it to make predictions. We will demonstrate how fitting a model to a dataset, and making predictions on the same samples used to fit the model, results in over-optimistic evaluation of the prediction error.

1a) Fitting a ridge regression model

As a first example, we are going to fit a ridge regression model to the entire dataset. Ridge is a regularised linear regression with a penalty term controlling the size of the coefficients. It is implemented in the glmnet() function of the glmnet package, with the argument alpha = 0. The strength of the penalty is controlled by the argument lambda. For now, we will fix lambda = 1. The glmnet() function takes the predictors as a matrix x and the response as a vector y, and predictions can be made with the predict() function, giving the predictors as the argument newx.

Fit a ridge regression with lambda = 1 to all of the observations.

# Your code here

Now use the predict() function to make predictions of the response, using the same data (X_mat) as an input to the newx argument.

# Your code here

Now we’ve made these predictions, we can evaluate them with performance metrics. For the purposes of this exercise, we will be using two metrics. The first is the root-mean-squared-error (RMSE) between the true values $y_i$ and the predictions $\hat{y}_i$ :

\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \]

The lower this is, the better. Since we standardised the response, it has a standard deviation of 1, and the RMSE can be read as follows : if RMSE = 1, your predictions are as good on average as a model which simply predicts the average response (i.e. one which does not use any of the predictor information). Anything below 1 indicates a model which predicts better than this naïve model.

The second metric is the R2 statistic, which is the squared Pearson correlation between the response and the predictions :

\[ R^{2} = cor_{\text{pearson}}(y, \hat{y})^{2} \]

The higher this is, the better. It can be interpreted as “the proportion of variability in the true response which is explained by the predictions”.

Below is a helper function which you can use to compute these metrics, by inputting your true and predicted response values.

calculate_metrics <- function(true, predicted) {
  sRMSE = sqrt(mean((true - predicted) ^ 2))/sd(true)
  R2 = cor(true, predicted, method = "pearson") ^ 2
  
  metrics = list("sRMSE" = sRMSE, "R2" = R2)
  
  return(metrics)
}

Compute and report the two metrics for your first ridge regression fit on the entire dataset.

# Your code here

We could also interpret our predictions by visualising the true versus predicted values. The below helper function allows you to do this, with the arguments being the true values, predicted values, and a string to use as a descriptive title for the figure.

plot_predictions <- function(true, predicted, title = "Insert a descriptive plot title here") {
  # The same limits for both axes : the range of the true and predicted values, with some padding
  results_range <- range(c(true, predicted))
  padding <- diff(results_range) * 0.02
  lims <- c(results_range[1] - padding, results_range[2] + padding)

  ggplot(data.frame(true, predicted), aes(x = true, y = predicted)) +
    geom_point(alpha = 0.7, size = 2) +
    geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", linewidth = 0.8) +
    coord_equal(xlim = lims, ylim = lims) +
    labs(title = title, x = "True (observed)", y = "Predicted") +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"))
}

Plot your results using the above helper function, and save the resulting image with a descriptive title, using the ggsave(), with here() to build the path (you should do this for all the figures in the practical). What do you conclude about the quality of the predictions?

# Your code here

1b) A training set and a testing set

The predictions above are made on the data which was used to fit the model, so they say nothing about how the model would perform on new participants. The simplest way to evaluate the model on unseen data is a train-test split. Observations are randomly assigned to either the training set or the testing set. A model is trained on the training set, and used to predict on the unseen testing set. The predictions can then be compared to their true values to obtain an unbiased estimate of the model’s performance.

Randomly select 70% of the participants as the training set (hint : the sample() function), and keep the remaining 30% as the testing set. Fit the same ridge regression model (lambda = 1) on the training set only, and predict the response of the testing set. Compute the same two metrics on the testing set, and plot the predictions.

How do the results compare with those of 1a? What does this tell us?

# Your code here

2) K-fold cross-validation

The train-test split uses the data inefficiently : only 30% of the participants are predicted, and the result depends on the random split. K-fold cross-validation is a better alternative. Here, we first split the data into K distinct folds (partitions of the data). Cross-validation is a loop, where in each stage, a different fold is chosen to act as the testing set, and a model is trained on all the other folds. Therefore, K-fold cross-validation results in a series of K models, and in a prediction on unseen data for every participant. In this section, we are going to improve on our previous train-test split by implementing 10-fold CV.

First, create the folds with the createFolds() function (use the arguments k = 10 and list = TRUE, which gives a list containing the indices of the participants in each fold).

# Your code here

Now, you will write a basic implementation of K-fold CV using a for loop over the 10 folds. In the ith iteration of the loop, form the training data by removing the observations belonging to the ith fold from the data. Then, fit a ridge regression model to this training data and use it to make predictions, as we did in question 1. Hint: you will need to store the predictions after every iteration. This can be done by creating an empty vector with rep(NA_real_, n) prior to initiating the for loop, and filling it with the predictions using the indices of the test participants within every iteration.

# Your code here

Compute the same two metrics and plot the predictions. How do the results compare with those of 1b?

3) Tuning hyperparameters with `caret`

Some models have hyperparameters to specify, which are parameters which control certain aspects of the learning process. Ridge regression, for example, has the hyperparameter lambda, which controls the strength of the penalisation term. Since we are unlikely to have much reason to prefer any particular values of these hyperparameters, we may prefer to let the data decide for us which are the best values. We can also do this with cross-validation : we try different values of the hyperparameters, and we keep those which give the best predictions (here, the lowest RMSE).

However, we must be careful. If we choose the hyperparameters using the same predictions that we use to evaluate the model, the performance will be over-optimistic (the same problem we evoked in question 1). Instead, we typically do something called nested cross-validation in order to first build a model with optimal hyperparameters, and then use that model to predict on unseen data. A nested cross-validation has two loops. The first, called the outer loop, separates the training data from the testing data (which is what you implemented above). The second, called the inner loop, uses cross-validation on the training data only to select the best values of the model hyperparameters. The testing data is never used to choose the hyperparameters.

Coding the inner loop from scratch is tedious, especially when the model has several hyperparameters. Luckily, the caret package exists to automate hyperparameter tuning. Caret standardises syntax across packages to allow for many predictive models to be easily implemented. To train models and find the best hyperparameter(s), we use the train() function in caret. This function takes the following main arguments :

x : the matrix or data frame of predictors (columns) by samples (rows)

y : the vector of responses

method : the name of the regression method to be used. A full list of methods implemented in caret can be found here. Ridge regression is implemented in glmnet, so the method name is “glmnet”.

tuneGrid : a data frame where each column is a named hyperparameter from the method, and each row is a unique combination of these hyperparameters to be used to fit a model to the training data. The expand.grid() function creates such a data frame from the values to try for each hyperparameter. For the “glmnet” method, the hyperparameters are alpha (0 for ridge regression) and lambda.

trControl : an object created using the trainControl() function which defines how the tuning process should be performed. For example, trainControl(method = "cv", number = 5) implements 5-fold (inner) CV for deciding the hyperparameters.

metric : the metric to be optimised when choosing the best model. We will use the root mean squared error (metric = “RMSE”).

The result of train() is a model which has been fitted on all the training data with the best hyperparameters. It can be used with predict() as before (with the testing predictors as the argument newdata), and the best hyperparameters are stored in its element $bestTune.

You are going to add 5-fold hyperparameter tuning in the ridge regression CV you implemented above. All this will require is to replace the call to glmnet() above with a call to caret::train(), using the arguments method = "glmnet", tuneGrid = ridge_grid, and trControl = trainControl(method = "cv", number = 5). As before, you may make predictions in the test set with the predict() function. Also, store the tuned values of the hyperparameter lambda, extracted from the bestTune$lambda of the caret::train() object. You may use the grid of hyperparameter options suggested below as inputs (its values are spread over several orders of magnitude, since lambda acts on a multiplicative scale).

lambda_grid <- c(0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 2, 3, 4, 5, 10, 20, 50, 100)

ridge_grid = expand.grid(alpha = 0, lambda = lambda_grid)

As before, compute the sRMSE and R2 of your predictions. Do these differ much from the non-tuned approach in question 2, where we simply fixed the value of lambda at 1?

Is the value of lambda selected the same in every iteration of the outer loop? What does this tell us?

# Your code here

4) Other models

Now that you have a nested cross-validation framework with caret, it is easy to try other models : you only need to change the method and the tuneGrid. For each of the following models, repeat your code from the previous section (with the same outer folds), and choose the hyperparameter values to try yourselves. Compute the metrics, plot the predictions, and save the figure. To compare the models later, store the predictions and metrics of each model in objects with distinct names (e.g. pred_lasso and metrics_lasso).

How should you choose a hyperparameter grid? There is no rule, and a grid is always somewhat arbitrary, but here is some practical advice. (1) Start with a small, coarse grid : every extra value multiplies the computation time (the number of models fitted is the number of rows of the grid, times 5 inner folds, times 10 outer folds). (2) For penalty parameters such as lambda, which act on a multiplicative scale, use values spread over several orders of magnitude (e.g. 0.001, 0.01, 0.1, 1, 10) rather than evenly spaced values. (3) After tuning, look at the selected values (the $bestTune of the train() object) : if a selected value is at the edge of your grid, the best value may lie beyond it, so extend the grid in that direction. Each model below comes with an example grid to start from.

4a) Lasso regression

Lasso regression is another regularised method which has the advantage over ridge regression that it can reduce model coefficients fully to 0. It is implemented in the "glmnet" method with alpha = 1. The hyperparameter to tune is lambda. For example, you could use the grid expand.grid(alpha = 1, lambda = c(0.001, 0.01, 0.1, 1, 10)).

# Your code here

4b) Elastic net regression

Elastic net regression combines the ridge and lasso penalties. It is also implemented in the "glmnet" method, but this time the hyperparameter alpha can take values between 0 (ridge) and 1 (lasso), and is tuned together with lambda. For example, you could use the grid expand.grid(alpha = c(0, 0.5, 1), lambda = c(0.01, 0.1, 1, 10)), which contains all the combinations of the values of alpha and lambda.

# Your code here

4c) Principal component regression

Principal component regression (PCR) first computes the principal components of the predictors, and then regresses the response on the first few components. It is implemented in the "pcr" method (package pls). The hyperparameter to tune is the number of components to keep, ncomp. Note that it cannot be larger than the number of participants in the training set. For example, you could use the grid data.frame(ncomp = c(1, 5, 10, 20)).

# Your code here

4d) Partial least squares regression

Partial least squares (PLS) regression is similar to PCR, but the components are chosen to explain the covariance between the predictors and the response, rather than the variance of the predictors only. It is implemented in the "pls" method. Again, the hyperparameter to tune is the number of components to keep, ncomp. You could use the same example grid as for PCR : data.frame(ncomp = c(1, 5, 10, 20)).

# Your code here

4e) Random forest

A random forest is an ensemble of decision trees, each grown on a random subset of the participants, where only a random subset of the genes is considered at each split. The predictions of the trees are averaged. It is implemented in the "ranger" method (package ranger), which has three hyperparameters. You must give a value for each of them in your tuning grid :

mtry : the number of genes randomly considered at each split. This is the main hyperparameter to tune; it must be at most the number of genes.

splitrule : the criterion to split the nodes. For a regression, use "variance".

min.node.size : the minimum number of participants in a terminal node of a tree (e.g. 5 or 10).

For example, you could use the grid expand.grid(mtry = c(50, 200), splitrule = "variance", min.node.size = c(5, 10)).

# Your code here

5) Summary of predictive performance

We have now tried multiple regression approaches in a nested cross-validation framework in order to make predictions of the day 28-post-vaccination fold-change in antibody values using day-7 fold-change in gene expression levels.

Collect the metrics of all your models in a single table, and save it in your output folder.

# Your code here

In a few sentences summarise your conclusions on the predictive performance of these models. Did any model, or class of models, perform noticeably better than the rest?

Let’s say we choose the best model by that which produced the best predictions in our cross-validation framework. Comment on how this strategy could introduce bias into our results.

6) Interpreting the selected variables (optional)

In the previous exercises, we have focused on code implementation and on maximising the predictive accuracy of the models. However, an important part of predictive analysis is the interpretation of any predictive signatures. Different regression methods may have different natural ways to interpret the importance of the input variables for the response. For example, in multiple linear regression, we can interpret the effect size estimates, confidence intervals, and p-values for a test of non-zero effect size. Random forest has a metric called variable importance, which tells us the loss in predictive accuracy when a given variable’s values are randomly permuted.

In this section, we are going to use the example of a lasso regression. When using L1-regularised methods like lasso (and elastic net), we may choose to interpret the predictors which have non-zero coefficients (i.e. the selected variables). We will do this in two setups : a model fit on the entire dataset, as in question 1, and cross-validation, as in questions 2 to 4.

6a) Variables selected by a model fit on the entire dataset

Fit a tuned lasso regression model on the entire dataset (i.e. ignore cross-validation for now), using caret::train() with method = "glmnet", alpha = 1 and a 5-fold cross-validation to tune lambda. The final model, which corresponds to the best value of lambda, can be accessed from the $finalModel element of the train() object, and the best value of lambda from its $bestTune$lambda element.

Extract the coefficients of the final model at the best value of lambda with the coef() function, using the argument s for the value of lambda. The output of coef() is not particularly easy to use : convert it to a matrix with as.matrix(), and then to a data frame with one row per gene and its coefficient. Keep only the genes with a non-zero coefficient (excluding the intercept), sorted by decreasing absolute value of the coefficient. How many genes have non-zero coefficients?

# Your code here

Since we standardised the predictors, the coefficients of the genes are on the same scale, and their absolute values can be compared. Make a bar plot of the coefficients of the (up to) 30 genes with the largest absolute coefficients, sorted by coefficient, and save the figure. What does the sign of a coefficient mean for the relationship between the gene and the antibody response?

# Your code here

Can you think of a way to interpret the names of the selected genes, to see whether particular biological processes are over-represented?

# Your code here

6b) Variables selected in cross-validation

There is a limitation to the interpretation above : the selected genes depend on the particular data used to fit the model. If we had a slightly different sample of participants, would the same genes be selected?

Cross-validation gives us a way to assess this. In each iteration of the cross-validation, a model is fit on a different part of the data, and may select different genes. With L1-regularisation methods like lasso and elastic net, we can therefore interpret the proportion of folds in which each gene is selected. These approaches can broadly be called stability selection methods, as genes which are more frequently selected are more stable predictors.

Modify your lasso cross-validation from question 4a (with the same outer folds), such that in each iteration you store the names of the genes with non-zero coefficients (as in 6a, using the model fitted by train() on the training data of the iteration). Then, compute the proportion of the 10 folds in which each gene was selected, and sort the genes from the most to the least frequently selected.

# Your code here

Make a plot showing the proportion of selection on one axis and the gene names on the other axis, sorted from the most to the least frequently selected, for the genes selected in at least half of the folds. Save the figure.

# Your code here

7) Sparse partial least squares with `mixOmics` (optional)

Sparse partial least squares (sPLS) is an approach which uses L1 (lasso) regularisation on the weighting vectors in a partial least squares approach. That is, it finds sparse components (made from a few genes only) which are maximally covarying with the response variable. This gives a way to select variables which is different from the lasso : the selected genes are those which contribute to the components.

We are going to implement this model using the mixOmics package, from Bioconductor. Unfortunately, this implementation is not available in caret, so we will have to use mixOmics’s own functions to tune the model. If you have not installed it yet, run the following chunk.

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("mixOmics")

Note : we do not load mixOmics with library() : it replaces some functions of other packages (such as select() from dplyr). We will call its functions with the mixOmics:: prefix.

7a) Tuning the sPLS model

The sPLS model has two hyperparameters. The first is called keepX, which is the number of genes to keep for each component. This is roughly equivalent to tuning the strength of the penalty term, like we have done previously in this practical. The second is called ncomp, which is the number of components to keep.

Read the sPLS (regression) example on the mixOmics site, and use the mixOmics::tune.spls() function on the entire dataset to determine the best values of ncomp and keepX (again, ignore the outer cross-validation for now : the function does its own cross-validation to tune the hyperparameters). Use the standardised genes X_mat as the predictors and the response y_vec. You could try up to 3 components, and test the numbers of genes to keep with, for example, test.keepX = c(5, 10, 20, 50, 100). Use the arguments mode = "regression", validation = "Mfold", folds = 5, measure = "MSE" and nrepeat = 3, and plot the tuning results. Save the plot.

Note : there is a known bug in the sPLS implementation in the regression setting, which may give you an error Unmapped Y contains samples with no associated class. This seems to be a problem when some values of the response are 0 : a workaround is to give the response as a matrix (as.matrix()) and to add a small value like 0.000001 to the response values, and the error should go away.

# Your code here

7b) sPLS in a cross-validation setup

In 7a, we tuned the hyperparameters of the sPLS model on the entire dataset. To evaluate sPLS fairly, and to compare it with the other models, we need to implement a nested cross-validation, as we did for the other models. Since sPLS is not available in caret, you have to code the steps yourself with the mixOmics functions.

Use the same 10 outer folds as before (folds), and write a for loop over the folds. In each iteration :

Tune ncomp and keepX on the training participants only, with mixOmics::tune.spls() as in 7a, and extract the best values.

Fit a sPLS model with these values on the training participants with mixOmics::spls().

Predict the response of the testing participants with the predict() function, giving the testing predictors as the argument newdata. The predictions are in the element $predict of the result, which is an array with three dimensions : participants, responses (here, just one) and components. Take the predictions for your response and for the number of components of the model.

Store the predictions of all the participants in a vector, as before. Then compute the metrics, plot the predictions, and save the figure.

# Your code here

Finally, add the metrics of the sPLS model to the table of the metrics of the other models (model_metrics, from section 5), and save the new table in your output folder. How does the sPLS model compare with the other models?

# Your code here

Task description

Setup

0a) Set up your project

Why do we need a project?

Create the project folder

Add the files

How do we refer to files?

0b) Install and load the required packages

0c) Check your project setup

0d) Load and explore the data

0e) Standardise the data

1) Building a predictive model

1a) Fitting a ridge regression model

1b) A training set and a testing set

2) K-fold cross-validation

3) Tuning hyperparameters with caret

4) Other models

4a) Lasso regression

4b) Elastic net regression

4c) Principal component regression

4d) Partial least squares regression

4e) Random forest

5) Summary of predictive performance

6) Interpreting the selected variables (optional)

6a) Variables selected by a model fit on the entire dataset

6b) Variables selected in cross-validation

7) Sparse partial least squares with mixOmics (optional)

7a) Tuning the sPLS model

7b) sPLS in a cross-validation setup

3) Tuning hyperparameters with `caret`

7) Sparse partial least squares with `mixOmics` (optional)