Get the quiz sheet for this module from the general Assessments page. Fill it in, then submit to the online grading system before the deadline.
For this exercise, we will implement some of the machine learning models covered in this module, and some of the general approaches we covered previously. This is a continuation of the data analysis exercise you’ve been working on over the last few weeks, analyzing flu data. The goal is that by the end of this exercise, you will have an example of a nice full analysis. So let’s get started!
Note: This exercise is likely time-consuming and you might get stuck. Please plan accordingly. Start early and ask for help in Slack.
This is a solo-exercise. Make sure your previous repo for this exercise is fully up-to-date and synced.
The goal is to have both a complete and well-organized data analysis example at the end of the exercise. So as you go through this exercise, in addition to doing the indicated tasks and writing the code, keep organizing things.
You want to remove anything that’s not part of the exercise (any files/folders from the template). You also want to make sure you have readme files or other documentation that briefly explain what files are where and in which order to run them. Make sure your R/Rmd scripts are easy to understand and that they all run.
As you re-organize, you can decide the mix of R and Rmd files and what you want to combine and what to split. However you do it, make sure it’s documented and easy to understand.
We previously covered pre-processing, but haven’t specifically looked at that yet in an exercise. So let’s add that part. This can be added in various places to the code. You can for instance make it part of the cleaning code, or include it in the recipe part of your tidymodels workflow (or a mix).
You probably realized that some of the variables contain the same information in slightly different ways. Specifically, Weakness, Cough and Myalgia exist on both a severity score and as Yes/No. Further, there are 2 variables for cough yes/no. These variables are strongly correlated and thus at a minimum, don’t help when you model, and are actually more likely to confuse your model and lead to errors/warnings/poor model performance (e.g. “predictions from a rank-deficient fit may be misleading). So let’s fix that. For those symptoms where you have both multiple levels and yes/no, remove all the yes/no versions. That should remove 4 variables.
We also want to code the 3 symptom severity factors as ordinal, so make sure they are coded as ordered factors. The order should of course be None < Mild < Moderate < Severe.
If you look at your data, you’ll see that some predictors are fairly unbalanced, with most patients reporting No and only a few Yes. If almost everyone is in one category and almost nobody in others, that often (but not always) means those predictor variables are not very helpful in fitting/predicting the outcome. Furthermore, if you do cross-validation and one of your samples happen to not include one predictor level at all, things might go wrong when applied to the holdout-sample. (There are methods to deal with this specific problem, but they are more complicated than what we will do here.) Thus, it is worth considering if we want to remove them. The recipes package in tidymodels has the function step_nzv() which can do that for you automatically. If you have lots of predictors, that might be good to use. But it’s often better to decide manually for each variable based on your scientific expertise if you want to remove it or not. We’ll take that approach here. After looking at the data, we decide to remove those binary predictors that have <50 entries in one category (there are 2). Write code to remove them.
You should end up with a data frame that has 730 observations and 26 variables. This is the dataset we’ll use for modeling.
I suggest you put all the analysis code into one R Markdown file. Note that for this exercise, some of your code might take long to run. In such situations, it is often good to have an alternative setup, where some of the heavy computations are done by separate R scripts/functions, and the results saved, pulled into and displayed in the main R Markdown file. How you want to do it here is up to you, as long as the R Markdown file shows all the main results and your code/file(s) is/are well documented.
You will likely be able to recycle bits of our earlier code, but you are also welcome to start fresh.
To keep things a bit simpler, for this exercise we focus on a single outcome, the continuous, numerical value of Body temperature. Thus, we are fitting regression models here. Once you finished the whole workflow for that outcome, you are welcome to (optionally) do it again and run a classification for the categorical outcome of Nausea. Since we are doing a regression, and we don’t have any specific expert knowledge that tells us how we should build our performance metric, we’ll go with one of the usual ones, namely RMSE. But remember to always think about your performance metric and don’t just use the default without at least considering other options.
Start by setting the random seed to 123. This should make everything reproducible and everyone should get the same results.
Split the dataset into 70% training, 30% testing. Also use the outcome BodyTemp as stratification. This allows for more balanced outcome values in the train and test sets. See e.g., section 3 of the Get Started tutorial.
We want to do 5-fold cross-validation, 5 times repeated. (There’s no specific reason to do this 5x5 pattern, other than to show you that there are different ways to pick the sample, and that I want you to not use the default.) For the CV folds, we also want to stratify on BodyTemp, as we did for the main train/test split. Use the vfold_cv() function to create a resample object for the training data with these specifications.
Create a recipe for the data and fitting. You won’t need to do much, just make sure you code the categorical variables as dummy variables, otherwise things might not work smoothly. For that, you want to use the step_dummy() function and pick all nominal predictor variables (which are actually all predictor variables here, since the only continuous variable is our outcome).
Write some code to compute the performance of a null model, i.e. a “model” that doesn’t use any predictor information. For a continuous outcome and RMSE as our metric, a null model is one that always predicts the mean of the outcome. Compute the RMSE for both training and test data for such a “model”. We’ll use that later to compare it to the performance of our real models. Of course, we expect/hope our real models that use predictor information to be better. If they aren’t that means they are no good.
We’ll fit a tree, a LASSO model, and a random forest. I chose those because they are used in the tutorial on the tidymodels website. You can of course add further models. For the tree, see the Tune model parameters section of the Get Started tutorial. For LASSO and the random forest, check out the Case Study section of the Get Started tutorial. Note that you will need to adjust the code for our scenario since we have a continuous outcome.
If you follow the tutorial, you’ll likely use the packages rpart, glmnet and ranger to fit those 3 models. Make sure they are installed and loaded.
I suggest you write code for each model separately. A lot of the code will look similar, so once you got the first one set up, the other two should be easier. They mainly differ in the commands specifying the tuning parameters and the tuning grid.
Each of these models requires some tuning. For the choices regarding the tuning parameters, you can follow the examples. Most of the models have more things that can be tuned, but for now you can stick to what they show in the tutorial. Follow the examples by setting up a workflow, set a tuning grid, and then use the tune_grid() function to tune the model using cross-validation.
Note that the tuning part, i.e., calling tune_grid() might take a good bit of time to run (possibly minutes).
The steps (block of code) you should have here are 1) model specification, 2) workflow definition, 3) tuning grid specification and 4) tuning using cross-validation and the tune_grid() function.
Once you have done the tuning, you can take a look at some diagnostics by sending your object returned from the tune_grid() function to autoplot(). For instance if you tuned the tree and saved the result as tree_tune_res, you can run tree_tune_res %>% autoplot(). Depending on the model, the plot will be different, but in general it shows you what happened during the tuning process.
Next, you want to get the model that the tuning process has determined is the best. You can get the best-fit model with select_best() and finalize_workflow() and then do one more fit to the training data with this final workflow using the fit() function. Follow the examples in the tutorial.
To evaluate the final fit for each model, do the following.
Make two plots, one that shows model predictions from the tuned model versus actual outcomes, and one that plots residuals. The actual outcomes you get straight from the data, the predicted outcomes you can get by applying the predict() function to the final fit.
Look at/print the model performance and compare it with the null model (still only on training data). Here, we want the performance of the tuned, best-fitting model on the CV dataset (we are not yet touching the test data). You can get that for instance with the show_best() function, which gives you the mean cross-validated performance for the best models. It also shows the standard deviation for the performance. Compare that model performance with the null model.
The mean and standard deviation of the performance give you a measure of overall performance and variability in that measure. The plots show you if there are any systematic deviations between model and data. Taken together, these can be compared for the different models and based on those (and as wanted, other considerations) a final model can be chosen.
Implement the model tuning/fitting and evaluating steps for all 3 models.
Once you have implemented above steps for the 3 models, you should a “best” fit for each one based on the tuning process. For each best model you should have performance, uncertainty around the performance measure, and some diagnostic plots. While for any real research project, you likely want to look deeper (e.g. at uncertainty in predictions instead of just overall performance), for now this is enough. Pick one of the three models. Explain why you pick it. There is no single answer that’s correct, I just want you to learn to reason/argue for why you are doing something: in this case justify why you are picking the model you do.
Once you picked your final model, you are allowed to once – and only once – fit it to the test data and check how well it performs on that data. This gives you a somewhat honest estimate of how the model might perform for new, unseen data. You can do that using the last_fit() function applied to the model you end up choosing. For the final model applied to the test set, report performance and the diagnostic plots as above.
And that concludes what is likely a fairly long exercise. The code itself is not that long, but it will take you time to cobble it together from the tidymodel tutorial and possibly other sources.
Make sure your analysis and all results are nicely documented and when knitted, your main analysis Rmd file produces a nice-looking HTML page. In a future exercise, you will add this to your portfolio. No need to post about the exercise, I know where to find it.
Let’s use this week’s discussion to talk about projects. Post a short summary of your project, very briefly describing your data, your question(s), your (planned) methods and (expected) results. Sort of like an abstract for a paper. This way, everyone can get a quick glimpse as to what everyone else is doing. Also provide a link to the repo for easy access.
Then, also post any specific questions/struggles/concerns you might have regarding your project. It’s quite likely that there is some overlap in approaches and questions among you and your classmates. Hopefully, through this discussion you can provide each other with some help and input.
Post this by Wednesday. Then read the summaries of your classmates’ projects. If you see a question/topic that you think you might be able to help with, either by answering a specific question, or by providing some general feedback, do so.
And of course, as you look at each others projects, it’s a 2-way street. You can provide help/feedback, but you are also welcome to take inspiration from what others are doing and integrate some of that into your own project.
Since you get feedback from me at each of your submission points, I plan to stay out of this week’s discussion 😁. So this is all you helping each other. Each of you will be reviewing 2 projects in a lot of depth at the end, but hopefully with the discussion this week you get a bit of an idea of what everyone else is doing, and you can provide feedback to others and earlier.
Comments for specific models
Here are some more suggestions and hints. Most of this is optional but worth trying.
For the tree model, if you want to plot the tree, you can use the
rpart.plotpackage and run this commandrpart.plot(extract_fit_parsnip(best_tree_fit)$fit)(assuming your result from the final workflow fit is calledbest_tree_fit). You might get a warning message, but the tree will show. You will likely find when you look at the actual/predicted plot or the residual plot that the tree model does not perform very well, and the model only predicts a few discrete outcome values. That’s also noticeable when you compare RMSE for the tree model and the null model, they are very similar.For the lasso model, you will likely find that it performs a bit better than the tree, but not a lot. If you want to see a plot for how the number of predictors included in the LASSO model changes with the tuning parameter, you can extract the model from your final fit (say it’s called
best_lasso_fit) withx <- best_lasso_fit$fit$fit$fitand thenplot(x, "lambda"). I know, this is awful code having to dig that deep into thebest_lasso_fitobject. Maybe there’s a better way of doing this with sometidymodelscommand (potentially some combination of$fitandextract_fit_parnsip()). If you know one, use it.For the random forest model, you will likely again find that it performs a bit better than the tree, but not a lot. The tuning setup might require some fiddling around, I had a few initial tries where the whole tuning failed. For a model like random forest, all variables stay in the model. There are ways to look at the variables that are most important. If you want to do that, you again need to pull out the fit object. Say it’s called
best_rf_fit, you can do that withx <- best_rf_fit$fit$fit$fitand then use thevip()function from thevippackage to plot the importance of the variables.