# 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](https://immunespace.org/query/study/SDY1119#study-citations). 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 : > 1. In RStudio, click `File > New Project... > New Directory > New Project`. > 2. Give the project a name, for example `omics_practicals`, and choose where to save it on your computer. Click `Create Project`. > 3. 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. ```{r, eval = F} 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 > 1. Copy the data file I sent you, `SDY1119_day7_GE_foldChange_top1000.RData`, into the `data-raw` folder. > 2. Copy this practical (the `.Rmd` file) into the `analysis` folder. > 3. 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`. ```{r, echo = F, results = 'hide', message = F} # This sets the options for the outputs when knitting the rmd file knitr::opts_chunk$set(warning = F, message = T, fig.width = 8, fig.height = 8) ``` ```{r, eval = F, results = 'hide', message = FALSE, warning = FALSE} 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. ```{r, echo = F, results = 'hide', message = F} library(dplyr) library(caret) library(glmnet) library(ggplot2) library(pls) library(ranger) library(here) # Default look of all the ggplot2 figures theme_set(theme_minimal(base_size = 14)) # Set a seed for reproducibility set.seed(123) ``` ## 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. ```{r} # Path to the root of the project here() # Check that the data file exists file.exists(here("data-raw", "SDY1119_day7_GE_foldChange_top1000.RData")) ``` > 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. ```{r} 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. ```{r} df_all[, 1:8] %>% head() ``` > 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. ```{r} 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. ```{r} # The response, as a numeric vector y_vec <- as.numeric(scale(df_all$response)) # The genes (all the columns after the response), as a numeric matrix gene_cols <- (which(colnames(df_all) == "response") + 1):ncol(df_all) X_mat <- df_all[, gene_cols] %>% as.matrix() %>% scale() # Checks : the means should be 0, the standard deviations 1, and there should be no missing values c(mean = mean(y_vec), sd = sd(y_vec)) summary(colMeans(X_mat)) summary(apply(X_mat, 2, sd)) sum(is.na(X_mat)) ``` > 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. ```{r} # 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. ```{r} # Fit the ridge regression to all the participants fit_all <- glmnet(x = X_mat, y = y_vec, alpha = 0, lambda = 1) ``` > Now use the `predict()` function to make predictions of the response, using the same data (`X_mat`) as an input to the `newx` argument. ```{r} # Predict the response of the same participants pred_train <- as.numeric(predict(fit_all, newx = X_mat)) ``` > 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. ```{r} 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. ```{r} metrics_train <- calculate_metrics(y_vec, predicted = pred_train) metrics_train ``` > **Answer:** The sRMSE is `r sprintf("%.3f", metrics_train$sRMSE)` and R2 is `r sprintf("%.3f", metrics_train$R2)`, indicating very high performance. > 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. ```{r} 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? ```{r} # Plot and save the figure plot_train <- plot_predictions(y_vec, pred_train, "Ridge regression : predictions on the training data") ggsave(here("output", "predictions_ridge_train.pdf"), plot_train, width = 8, height = 8) plot_train ``` > **Answer:** The predictions are effectively perfect. This is an indication of overfitting: the model exploits the complexity of the data to build a perfect relationship between the predictors and the response. However, these relationships must be tested on new data to be truly evaluated, as this is likely over-optimistic. ## 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? ```{r} n <- nrow(X_mat) # Randomly select 70% of the participants for training, and keep the others for testing train_idx <- sample(n, size = round(0.7 * n)) test_idx <- setdiff(seq_len(n), train_idx) # Fit the model on the training set only fit_split <- glmnet(x = X_mat[train_idx, ], y = y_vec[train_idx], alpha = 0, lambda = 1) # Predict the testing set pred_test <- as.numeric(predict(fit_split, newx = X_mat[test_idx, ])) # Evaluation metrics, on the testing set metrics_test = calculate_metrics(true = y_vec[test_idx], predicted = pred_test) metrics_test # Plot and save the figure plot_test <- plot_predictions(y_vec[test_idx], pred_test, "Ridge regression : predictions on a testing set") ggsave(here("output", "predictions_ridge_split.pdf"), plot_test, width = 8, height = 8) plot_test ``` > **Answer:** Now, the metrics are sRMSE = `r sprintf("%.3f", metrics_test$sRMSE)` and R2 = `r sprintf("%.3f", metrics_test$R2)` (this will depend on your data split). While these indicate some reasonable level of predictive performance, they are not perfect, as before. This is more likely to be an accurate indication of the predictive signal in the data. # 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). ```{r} # Create the 10 folds : a list containing, for each fold, the indices of the participants in it folds <- createFolds(y_vec, k = 10, list = TRUE, returnTrain = FALSE) ``` > Now, you will write a basic implementation of K-fold CV using a `for` loop over the 10 folds. In the `i`th iteration of the loop, form the training data by removing the observations belonging to the `i`th 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. ```{r} # Vector which will contain the prediction of each participant pred_cv <- rep(NA_real_, n) for (i in seq_along(folds)) { # The testing set is fold i, and the training set is all the other participants test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) # Fit the model on the training set, and predict the testing set fit_fold <- glmnet(x = X_mat[train_idx, ], y = y_vec[train_idx], alpha = 0, lambda = 1) pred_cv[test_idx] <- as.numeric(predict(fit_fold, newx = X_mat[test_idx, ])) } # Evaluation metrics, on the predictions of all the participants metrics_cv <- calculate_metrics(true = y_vec, predicted = pred_cv) metrics_cv # Plot and save the figure plot_cv <- plot_predictions(y_vec, pred_cv, "Ridge regression (lambda = 1) : cross-validated predictions") ggsave(here("output", "predictions_ridge_cv.pdf"), plot_cv, width = 8, height = 8) plot_cv ``` > Compute the same two metrics and plot the predictions. How do the results compare with those of 1b? > **Answer:** We have sRMSE = `r sprintf("%.3f", metrics_cv$sRMSE)` and R2 = `r sprintf("%.3f", metrics_cv$R2)`. These are in the same order of magnitude as the train test split, however using CV we use all of the data and reduce the impact of the random data split on the results. # 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](https://topepo.github.io/caret/available-models.html). 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). ```{r} 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? ```{r} pred_ridge_tuned <- rep(NA_real_, n) best_lambda_ridge_tuned <- rep(NA_real_, length(folds)) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) # Tune lambda on the training set only, with a 5-fold inner cross-validation fit_ridge_tuned <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "glmnet", tuneGrid = ridge_grid, trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) # Predict the testing set with the model fitted using the best lambda pred_ridge_tuned[test_idx] <- predict(fit_ridge_tuned, newdata = X_mat[test_idx, ]) best_lambda_ridge_tuned[i] <- fit_ridge_tuned$bestTune$lambda } best_lambda_ridge_tuned metrics_ridge_tuned <- calculate_metrics(true = y_vec, predicted = pred_ridge_tuned) metrics_ridge_tuned plot_ridge_tuned <- plot_predictions(y_vec, pred_ridge_tuned, "Ridge regression : nested cross-validation") ggsave(here("output", "predictions_ridge_tuned.pdf"), plot_ridge_tuned, width = 8, height = 8) plot_ridge_tuned ``` > **Answer:** The metrics are sRMSE = `r sprintf("%.3f", metrics_ridge_tuned$sRMSE)` and R2 = `r sprintf("%.3f", metrics_ridge_tuned$R2)`. These are effectively the same as the predictions using the fixed hyperparameter approach above, suggesting that the impact of the choice of hyperparameters may not be of great importance. `r if (length(unique(best_lambda_ridge_tuned)) == 1) paste0("The optimal hyperparameter value was fixed at ", unique(best_lambda_ridge_tuned), " across iterations, suggesting a stable hyperparameter selection procedure.") else paste0("The optimal hyperparameter value varied between ", min(best_lambda_ridge_tuned), " and ", max(best_lambda_ridge_tuned), " across iterations, suggesting that its selection is sensitive to the training data.")` # 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))`. ```{r} pred_lasso <- rep(NA_real_, n) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) fit_lasso <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "glmnet", tuneGrid = expand.grid(alpha = 1, lambda = lambda_grid), trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) pred_lasso[test_idx] <- predict(fit_lasso, newdata = X_mat[test_idx, ]) } metrics_lasso <- calculate_metrics(true = y_vec, predicted = pred_lasso) metrics_lasso plot_lasso <- plot_predictions(y_vec, pred_lasso, "Lasso regression : nested cross-validation") ggsave(here("output", "predictions_lasso.pdf"), plot_lasso, width = 8, height = 8) plot_lasso ``` ## 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`. ```{r} pred_enet <- rep(NA_real_, n) # Both alpha and lambda are tuned : the grid contains all their combinations enet_grid <- expand.grid(alpha = seq(0, 1, 0.25), lambda = lambda_grid) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) fit_enet <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "glmnet", tuneGrid = enet_grid, trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) pred_enet[test_idx] <- predict(fit_enet, newdata = X_mat[test_idx, ]) } metrics_enet <- calculate_metrics(true = y_vec, predicted = pred_enet) metrics_enet plot_enet <- plot_predictions(y_vec, pred_enet, "Elastic net regression : nested cross-validation") ggsave(here("output", "predictions_enet.pdf"), plot_enet, width = 8, height = 8) plot_enet ``` ## 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))`. ```{r} pred_pcr <- rep(NA_real_, n) # The numbers of components to try ncomp_grid <- data.frame(ncomp = c(1, 2, 5, 10, 20)) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) fit_pcr <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "pcr", tuneGrid = ncomp_grid, trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) pred_pcr[test_idx] <- predict(fit_pcr, newdata = X_mat[test_idx, ]) } metrics_pcr <- calculate_metrics(true = y_vec, predicted = pred_pcr) metrics_pcr plot_pcr <- plot_predictions(y_vec, pred_pcr, "Principal component regression : nested cross-validation") ggsave(here("output", "predictions_pcr.pdf"), plot_pcr, width = 8, height = 8) plot_pcr ``` ## 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))`. ```{r} pred_pls <- rep(NA_real_, n) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) fit_pls <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "pls", tuneGrid = ncomp_grid, trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) pred_pls[test_idx] <- predict(fit_pls, newdata = X_mat[test_idx, ]) } metrics_pls <- calculate_metrics(true = y_vec, predicted = pred_pls) metrics_pls plot_pls <- plot_predictions(y_vec, pred_pls, "Partial least squares regression : nested cross-validation") ggsave(here("output", "predictions_pls.pdf"), plot_pls, width = 8, height = 8) plot_pls ``` ## 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))`. ```{r} pred_rf <- rep(NA_real_, n) rf_grid <- expand.grid( mtry = c(100, 300), splitrule = "variance", min.node.size = c(5, 10) ) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) fit_rf <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "ranger", tuneGrid = rf_grid, trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) pred_rf[test_idx] <- predict(fit_rf, newdata = X_mat[test_idx, ]) } metrics_rf <- calculate_metrics(true = y_vec, predicted = pred_rf) metrics_rf plot_rf <- plot_predictions(y_vec, pred_rf, "Random forest : nested cross-validation") ggsave(here("output", "predictions_rf.pdf"), plot_rf, width = 8, height = 8) plot_rf ``` # 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. ```{r} model_metrics <- bind_rows( list( "Ridge" = metrics_ridge_tuned, "Lasso" = metrics_lasso, "Elastic net" = metrics_enet, "PCR" = metrics_pcr, "PLS" = metrics_pls, "Random forest" = metrics_rf ), .id = "model" ) write.csv(model_metrics, here("output", "model_metrics.csv"), row.names = FALSE) model_metrics ``` > 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? > **Answer:** The models generally produced comparable results. The best results were obtained with the `r model_metrics$model[which.min(model_metrics$sRMSE)]` approach, producing sRMSE = `r sprintf("%.3f", model_metrics$sRMSE[which.min(model_metrics$sRMSE)])` and R2 = `r sprintf("%.3f", model_metrics$R2[which.min(model_metrics$sRMSE)])`. However, these gains were marginal, suggesting that the choice of predictive model is less important than the inherent signal present in the data. > 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. > **Answer:** Choosing the model with the best cross-validated performance, and then reporting that performance as our estimate of how well it would predict new participants, is over-optimistic. In our case, the models have very similar performances, so the apparent advantage of the best model could easily be due to chance. To obtain an unbiased estimate of the performance of the whole procedure, including the choice of the model, we need data which played no part in this choice, e.g. a separate testing set. # 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? ```{r} # Fit and tune the lasso model on the whole dataset lasso_full <- train( x = X_mat, y = y_vec, method = "glmnet", tuneGrid = expand.grid(alpha = 1, lambda = lambda_grid), trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) # The coefficients of the final model at the best lambda, as a matrix lasso_coef_mat <- as.matrix(coef(lasso_full$finalModel, s = lasso_full$bestTune$lambda)) # One row per gene with a non-zero coefficient (excluding the intercept), sorted by decreasing size lasso_coefs <- data.frame( gene = rownames(lasso_coef_mat), coefficient = lasso_coef_mat[, 1] ) %>% filter(coefficient != 0, gene != "(Intercept)") %>% arrange(desc(abs(coefficient))) message(nrow(lasso_coefs), " genes were selected.") head(lasso_coefs, 10) ``` > **Answer:** `r nrow(lasso_coefs)` genes (out of `r ncol(X_mat)`) have non-zero coefficients. > 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? ```{r} plot_lasso_coefs <- head(lasso_coefs, 30) %>% ggplot(aes(x = reorder(gene, coefficient), y = coefficient)) + geom_col() + coord_flip() + labs(x = "Gene", y = "Coefficient (standardised)", title = "Lasso regression : largest coefficients") ggsave(here("output", "lasso_coefficients.pdf"), plot_lasso_coefs, width = 8, height = 8) plot_lasso_coefs ``` > **Answer:** A positive coefficient means that a higher day-7 fold-change in the expression of the gene is associated with a higher antibody response, all the other genes being equal, and a negative coefficient means the opposite. > Can you think of a way to interpret the names of the selected genes, to see whether particular biological processes are over-represented? ```{r} paste(lasso_coefs$gene, collapse = ", ") ``` ## 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. ```{r} # List which will contain, for each fold, the names of the selected genes lasso_selected <- vector("list", length(folds)) for (i in seq_along(folds)) { train_idx <- setdiff(seq_len(n), folds[[i]]) # Tune the lasso on the training set, with a 5-fold inner cross-validation fit_lasso_fold <- train( x = X_mat[train_idx, ], y = y_vec[train_idx], method = "glmnet", tuneGrid = expand.grid(alpha = 1, lambda = lambda_grid), trControl = trainControl(method = "cv", number = 5), metric = "RMSE" ) # The genes with non-zero coefficients in the final model of this fold coef_mat_fold <- as.matrix(coef(fit_lasso_fold$finalModel, s = fit_lasso_fold$bestTune$lambda)) lasso_selected[[i]] <- setdiff(rownames(coef_mat_fold)[coef_mat_fold[, 1] != 0], "(Intercept)") } # The proportion of folds in which each gene was selected, sorted from the most to the least frequent selection_freq <- table(unlist(lasso_selected)) / length(folds) selection_df <- data.frame( gene = names(selection_freq), prop = as.numeric(selection_freq) ) %>% arrange(desc(prop)) head(selection_df, 10) ``` > 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. ```{r} plot_selection <- selection_df %>% filter(prop >= 0.5) %>% ggplot(aes(x = reorder(gene, prop), y = prop)) + geom_col() + coord_flip() + labs(x = "Gene", y = "Proportion of folds in which the gene is selected", title = "Selection frequency of the genes in cross-validation (lasso)") ggsave(here("output", "lasso_selection_frequency.pdf"), plot_selection, width = 8, height = 8) plot_selection ``` # 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. ```{r, eval = F} 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](https://mixomics.org/parameter-tuning/), 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](https://mixomics-users.discourse.group/t/error-in-check-entry-pls-x-y-ncomp-keepx-keepy-mode-mode-scale-scale-unmapped-y-contains-samples-with-no-associated-class-may-be-caused-by-nas-in-input-y-vector/1202/7), 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. ```{r} # The response as a matrix, with the workaround for the known bug y_spls <- as.matrix(y_vec) + 0.000001 spls_tuning <- mixOmics::tune.spls( X = X_mat, Y = y_spls, ncomp = 3, test.keepX = c(5, 10, 20, 50, 100), validation = "Mfold", folds = 5, measure = "MSE", mode = "regression", nrepeat = 3 ) # Plot the tuning results, and save the plot as a PDF plot(spls_tuning) pdf(here("output", "spls_tuning.pdf"), width = 8, height = 6) print(plot(spls_tuning)) invisible(dev.off()) # The best number of components, and the best number of genes to keep for each component spls_best_ncomp <- spls_tuning$choice.ncomp[[1]] spls_best_keepX <- spls_tuning$choice.keepX[seq_len(spls_best_ncomp)] spls_best_ncomp spls_best_keepX ``` > **Answer:** The best number of components is `r spls_best_ncomp`, keeping `r paste(spls_best_keepX, collapse = ", ")` gene(s) for the component(s). ## 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. ```{r} pred_spls <- rep(NA_real_, n) for (i in seq_along(folds)) { test_idx <- folds[[i]] train_idx <- setdiff(seq_len(n), test_idx) # Tune ncomp and keepX on the training participants only tuning_fold <- mixOmics::tune.spls( X = X_mat[train_idx, ], Y = y_spls[train_idx, , drop = FALSE], ncomp = 3, test.keepX = c(5, 10, 20, 50, 100), validation = "Mfold", folds = 5, measure = "MSE", mode = "regression", nrepeat = 3 ) ncomp_fold <- tuning_fold$choice.ncomp[[1]] keepX_fold <- tuning_fold$choice.keepX[seq_len(ncomp_fold)] # Fit the model with the best values on the training participants fit_spls <- mixOmics::spls( X_mat[train_idx, ], y_spls[train_idx, , drop = FALSE], keepX = keepX_fold, ncomp = ncomp_fold ) # Predict the testing participants : the first response, with all the components of the model pred_fold <- predict(fit_spls, newdata = X_mat[test_idx, ]) pred_spls[test_idx] <- pred_fold$predict[, 1, ncomp_fold] } metrics_spls <- calculate_metrics(true = y_vec, predicted = pred_spls) metrics_spls plot_spls <- plot_predictions(y_vec, pred_spls, "Sparse PLS regression : nested cross-validation") ggsave(here("output", "predictions_spls.pdf"), plot_spls, width = 8, height = 8) plot_spls ``` > 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? ```{r} model_metrics_with_spls <- bind_rows( model_metrics, data.frame(model = "sPLS", sRMSE = metrics_spls$sRMSE, R2 = metrics_spls$R2) ) write.csv(model_metrics_with_spls, here("output", "model_metrics_with_spls.csv"), row.names = FALSE) model_metrics_with_spls ``` > **Answer:** The sPLS model has sRMSE = `r sprintf("%.3f", metrics_spls$sRMSE)` and R2 = `r sprintf("%.3f", metrics_spls$R2)`, compared with sRMSE = `r sprintf("%.3f", min(model_metrics$sRMSE))` for the best of the other models (`r model_metrics$model[which.min(model_metrics$sRMSE)]`).