set.seed(5039378)
split <- rsample::initial_split(df_dummy, prop = 0.75)
train <- rsample::training(split)
test <- rsample::testing(split)Modelling
TidyModelling
The R packages included in the Tidymodels provide an excellent framework for undertaking the modelling aspect of the work.
A very useful training workshop was held 2023-10-18 as part of the R in Pharma event. Nicola Rennie has made her material is available on-line and it provides a good starting point for applying the various modelling methodologies to prediction of Thyroid cancer.
Whilst this provides an excellent introduction to the Tidy Modelling framework the book Tidy Modeling with R goes deeper into the methods and options available. Another useful reference on which this books builds is R for Data Science (2e) which should serve as a useful reference for learning R and adopting good practices.
It is also recommended to read the documentation that goes with the Tidymodels package, in particular the Get Started page which includes a predictive modelling case study.
Setting up Tidymodels
In the absence of the data set that is to be analysed we simulate some dummy data on which to demonstrate the methods.
Training and Testing
To get going with Tidymodels we first need to split our data into Testing and Training subsets. This is done so that we do not have an over-fitted model as we fit the model to the training subset and then test its predictive accuracy in the subset that we withheld, the test subset.
The allocation of individuals to train or test is performed randomly and so we set a seed to ensure the pseudo-random number generator produces the same split each and every time this script is run (this makes our work reproducible). A decision has to be made about how to split the data, often a slightly larger proportion is used for the training subset, here we chose to use a 3:1 split (i.e. 75% of observations are used in the training set).
Cross Validation
We will still want to make some assessment of the model estimated from the training set and this is achieved by resampling which involves taking subsets of our training data and fitting the models on those subsets and then looking at the distribution of metrics across the multiple model fits. A commonly used approach for this methodology is Cross-Validation. V-fold cross-validation splits the data into V folds, the first fold is excluded from the data set and the model assessed, then this subset is replaced and the second fold is excluded and the model assessed again. This is repeated until each fold has been excluded and the model estimated on the remainder. This method can then be repeated R times where the folds are varied in each repetition to give a better estimate of the model parameters.
The following figure gives an overview of how V-fold cross-validation works (without repetition).
We can set this up by passing the train subset into the rsample::vflod_cv() function and define the number of folds (v) and the number of repeats (repeats).
cv_folds <- rsample::vfold_cv(train, v = 10, repeats = 10)Another method of cross-validation is Leave One Out where one observation is removed from the (training) data set, the model is fitted on the remaining samples and then used to make a prediction on the excluded sample. This is then repeated on all samples. We can set this up using the rsample::loo_cv() function.
cv_loo <- rsample::loo_cv(train)Recipe
In the Tidymodels framework the starting point is to create a Recipe, this sets up the “ingredients” for the model and defines what steps should be taken prior to fitting the model, regardless of what model is being fitted. Steps that can go into making a recipe are…
- Define the model in terms of outcome variable and predictors.
- Creating dummy variables for categorical variables.
- Normalizing data, e.g. log-transformation, rescaling.
thyroid_recipe <- recipes::recipe(final_pathology ~ gender + age + asa_score + smoking_status + nodule_fna_thy, data = train) |>
# recipes::step_num2factor(final_pathology, levels = c("Benign", "Malignant")) |>
recipes::step_dummy(gender, asa_score, smoking_status, nodule_fna_thy) |>
recipes::step_normalize(all_numeric())Workflow
Once a recipe has been defined it can be added to a workflow which will apply this step every time the workflow is run using the different models and post-processing steps that we will add.
thyroid_workflow <- workflows::workflow() |>
workflows::add_recipe(thyroid_recipe)Methods to Consider
There are a wealth of options when it comes to “Machine Learning”, these days even logistic regression is grouped into the term! There are however a large number of more sophisticated methods for analysing the data, and it is often wise to apply a number of methods to ensure they are converging on similar solutions rather than cherry-picking any one method.
Logistic Regression
Logistic regression, even in its most basic form is, still considered a “machine learning” algorithm. However because we do not know which out of an array of variable will be useful predictors and, following Occam’s Razor we would tend to prefer simpler explanatory models over complex ones we need a method of determining what subset of variables gives good prediction. An old approach to this was to use Stepwise regression, perhaps based on univariable analyses to select which to use as a base but these approaches have fallen out favour for various reasons as they are ultimate biased (for an overview see (Steyerberg et al. 2001)).
A popular alternative is the Least Absolute Shrinkage and Selection Operator (LASSO) proposed by (Tibshirani 1996) which performs L1 regularisation and allows the coefficients for variables in a series of fitted models to “shrink” towards zero but not drop out completely. It is similar to Ridge Regression which avoids over-fitting by reducing the sum of squares of the regression coefficients but unlike Ridge Regression it allows variables to be selected as the coefficients can (almost) drop out by virtue of their coefficients shrinking towards zero.
Specify the Model
Whilst we have defined the relationship between variables in the Worfklow above we now need to say what model we wish to use to test the relationship between variables.
We first set up a simple logistic regression model using parsnip. The mixture argument is a value 0 <= mixture <= 1 which determines how much L1 regularisation is used in the model. A value of mixture = 1 is equivalent to full L1 regularisation and a LASSO model whilst a value of mixture = 0 is equivalent to ridge regression.
tune_spec_lasso <- parsnip::logistic_reg(penalty = hardhat::tune(), mixture = 1) |>
parsnip::set_engine("glmnet")Tuning
The next step is to tune the model using the tune::tune_grid() function which will calculate accuracy or Root Mean Square Error (or other metrics) for a recipe with multiple samples. We defined our resamples above in two forms (not that we defined cross-validation as a way of sampling from our training data above).
lasso_grid <- tune::tune_grid(
object = workflows::add_model(thyroid_workflow, tune_spec_lasso),
resamples = cv_folds,
grid = dials::grid_regular(penalty(), levels = 50)
)Fit Final Model
Now select the model with the highest Receiver Operating Characteristic Area Under the Curve (ROC AUC).
lasso_kfold_roc_auc <- lasso_grid |>
tune::select_best(metric = "roc_auc")…and finalize the workflow by adding this value.
final_lasso_kfold <- workflows::finalize_workflow(
workflows::add_model(thyroid_workflow, tune_spec_lasso),
lasso_kfold_roc_auc
)Model Evaluation
Collect metrics using the original data.
tune::last_fit(object = final_lasso_kfold, split = split) |>
tune::collect_metrics()And plot the importance of variables
final_lasso_kfold |>
fit(train) |>
hardhat::extract_fit_parsnip() |>
vip::vi(lambda = highest_roc_auc_lasso$penalty) |>
dplyr::mutate(
Importance = abs(Importance),
Variable = fct_reorder(Variable, Importance)
) |>
ggplot(mapping = aes(x = Importance, y = Variable, fill = Sign)) +
geom_col()Random Forest
Random forests build on the concept of regression trees by building many such trees on random subsets of data and averaging them. They build “deep” trees partitioning each time until all individuals are classified, part of the art is deciding how “deep” to go as too many partitions results in over-fitting and lack of generalisability. It is in this regard an “ensemble” method as it is averaging across many possible methods.
Gradient Boosting
Gradient Boosting is another “ensemble” method but in contrast to Random Forests where deep partitioning trees are formed, gradient boosting uses “shallow” trees with simpler decision rules and layers are built in a stage-wise fashion. It allows a choice of loss functions for optimisation and typically out-performs Random Forests when it comes to prediction modelling.
Support Vector Machine
Model Assessment
When considering the utility of prediction there are a number of metrics that can be used to determine how useful a model is. These focus on a number of terms which have very specific statistical meaning such as Sensitivity and Specificity.
| Test Positive | Test Negative | |
|---|---|---|
| Actual Positive | True Positive (TP) | False Negative (FN) |
| Actual Negative | False Positive (FP) | True Negative (TN) |
There are a wealth of metrics based on combinations of these True/False Positive/Negative including precision, recall and so forth (refer to the above linked article for further details).
Receiver Operating Characteristics
The Receiver Operating Characteristic (ROC) is a useful graphical tool for potting the false positive rate v’s the true positive rate (sensitivity).
References
Steyerberg, Ewout W., Marinus J. C. Eijkemans, Frank E. Harrell, and J. Dik F. Habbema. 2001. “Prognostic Modeling with Logistic Regression Analysis: In Search of a Sensible Strategy in Small Data Sets.” Med. Decis. Making 21 (1): 45–56. https://doi.org/10.1177/0272989X0102100106.
Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society. Series B (Methodological) 58 (1): 267–88. https://doi.org/10.2307/2346178.