Prediction

Rules

• What it takes to develop a successful (clinical) prediction model. TEN important things to avoid

Stability

Stability of clinical prediction models developed using statistical or machine learning methods 2023

Feature selection

• https://topepo.github.io/caret/recursive-feature-elimination.html
• Feature Selection Approaches from http://r-statistics.co
• Leave-one-out feature importance (LOFO) + LightGBM

Recursive feature elimination

• Feature Recursive Elimination (FRE) is a feature selection algorithm developed by Isabelle Guyon and her colleagues. It is a recursive algorithm that is used to select a subset of relevant features from a large pool of candidate features in a dataset.
• The basic idea behind FRE is to recursively eliminate features that are least informative or least relevant to the prediction task. At each iteration, the algorithm removes the feature with the lowest contribution to the prediction performance, as measured by a performance metric such as accuracy or F1-score. The algorithm continues to eliminate features until a stopping criterion is met, such as a fixed number of features, a minimum prediction performance, or a user-defined stopping threshold.
• FRE has been applied to a variety of machine learning problems, including classification, regression, and clustering. It is often used in combination with other feature selection algorithms, such as wrapper methods or filter methods, to provide a more comprehensive and robust approach to feature selection.
• One advantage of FRE is that it provides a simple and straightforward way to select features that are most relevant to the prediction task. It also allows for the evaluation of the importance of individual features and provides a way to visualize the relationship between the features and the prediction performance. Another advantage of FRE is that it is computationally efficient and can handle large datasets with many features.
• R packages
  • caret: ?rfe
  • Boruta: ?Boruta
  • mlr: Feature Selection

• In mathematical terms, RFE can be formulated as follows:
  1. Initialize the feature set F to contain all features in the dataset X.
  2. Fit a model, such as a Support Vector Machine (SVM), to the data X using the current feature set F.
  3. Rank the features in F based on their importance, as determined by the model coefficients or other feature importance measures.
  4. Remove the feature with the lowest importance from the feature set F.
  5. Repeat steps 2-4 until a stopping criterion is reached.
  The stopping criterion can be defined as a fixed number of features to be included in the final model, a certain threshold of cross-validation accuracy, or a threshold of classification error, among others.
  At each iteration of the RFE process, the model is refitted with the remaining features, and the importance of each feature is re-evaluated. By removing the least important features at each iteration, the RFE process can identify a subset of the most important features that contribute to the prediction performance of the model.
• Backward elimination is not a special case of Recursive Feature Elimination (RFE). While both methods aim to select the most relevant features for a model, they use different approaches to achieve this goal. Backward elimination is specific to linear regression models and relies on p-values to determine which features to remove, while RFE can be used with any model that assigns weights to features or has a feature importance attribute and removes the weakest features based on their importance.
• When using Recursive Feature Elimination (RFE) with the Random Forest algorithm, RFE starts by training a Random Forest model on the entire set of features and computing the importance of each feature. The least important features are then removed from the current set of features, and the process is repeated on the pruned set until the desired number of features is reached. The importance of each feature can be determined by the Random Forest algorithm’s internal method for measuring feature importance.
• The number of features to select in Recursive Feature Elimination (RFE) is a hyperparameter that can be chosen based on the specific problem and dataset. There are several approaches to determining the optimal number of features:
  • Domain knowledge: If you have prior knowledge about the problem domain, you may have an idea of how many features are relevant and should be included in the model.
  • Cross-validation: You can use cross-validation to evaluate the performance of the model with different numbers of features and choose the number that results in the best performance.
  • RFECV: In scikit-learn, you can use the RFECV class, which performs RFE with cross-validation to automatically find the optimal number of features.

SVM RFE

• Support Vector Machines (SVMs) can be used to perform Recursive Feature Elimination (RFE). RFE is a feature selection method that involves iteratively removing the least important features from a dataset and re-fitting a model until a stopping criterion is reached. The goal of RFE is to identify a subset of the most important features that can be used to build a predictive model with good accuracy.
• SVMs are a type of machine learning algorithm that can be used for classification and regression tasks. They work by finding the hyperplane that maximizes the margin between the classes in a dataset. The hyperplane is defined by a subset of the features, called support vectors, that have the largest influence on the classification decision.
• To perform RFE with SVMs, one can use the support vectors as a measure of feature importance and remove the features with the smallest magnitude of coefficients in the SVM model. At each iteration, the SVM model is refitted with the remaining features and the process is repeated until a stopping criterion is reached.
• In this way, RFE with SVMs can be used to identify a subset of the most important features that contribute to the prediction performance of the SVM model. RFE with SVMs can also be used to handle high-dimensional datasets with many features, as it can help reduce the dimensionality of the data and improve the interpretability of the model.
• Common stopping criteria that are used in RFE with SVMs, including:
  • Number of Features: One common stopping criterion is to specify a fixed number of features to be included in the final model. For example, one can specify that the RFE process should stop when the number of features is reduced to a certain value, such as 10 or 50.
  • Cross-Validation Accuracy: Another common stopping criterion is to use cross-validation accuracy as a measure of performance. The RFE process can be stopped when the cross-validation accuracy reaches a certain threshold, or when it starts to decrease, indicating that further feature elimination is not beneficial for improving performance.
  • Classification Error: A third common stopping criterion is to use the classification error, or misclassification rate, as a measure of performance. The RFE process can be stopped when the classification error reaches a certain threshold, or when it starts to increase, indicating that further feature elimination is not beneficial for improving performance.
  • The choice of stopping criterion will depend on the specific requirements of the user and the characteristics of the dataset. It is important to select a stopping criterion that balances the need for feature reduction with the need to preserve performance and avoid overfitting.

• One example:

  library(caret)
  library(mlbench)
  require(randomForest) # require(e1071)
  data(Sonar)
  X <- Sonar[, 1:60]
  y <- Sonar[, 61]
  ctrl <- rfeControl(functions = rfFuncs, method = "cv", number= 10) # 10-fold CV
      # ctrl <- rfeControl(functions = svmFuncs, method = "cv", number =10)
  rfe_result <- rfe(x = X, y= y, sizes=c(1:10), rfeControl =ctrl) 
      # range of feature subset sizes 
  print(rfe_result)
  plot(rfe_result)
  

• Another example:

  library(caret)
  data(iris)
  
  # Split the data into training and testing sets
  set.seed(123)
  train_index <- createDataPartition(iris$Species, p = 0.8, list = FALSE)
  train_data <- iris[train_index, ]
  test_data <- iris[-train_index, ]
  
  # Preprocess the data
  preProcess_obj <- preProcess(train_data[, -5], method = c("center", "scale"))
  train_data_preprocessed <- predict(preProcess_obj, train_data[, -5])
  train_data_preprocessed$Species <- train_data$Species
  
  # Perform Recursive Feature Elimination with a SVM classifier
  ctrl <- rfeControl(functions = rfFuncs, method = "repeatedcv", repeats = 3, verbose = FALSE)
  svm_rfe_model <- rfe(x = train_data_preprocessed[, -5], 
                       y = train_data_preprocessed$Species, 
                       sizes = c(1:4), 
                       rfeControl = ctrl, 
                       method = "svmLinear")
  print(svm_rfe_model)
  # Recursive feature selection
  #
  # Outer resampling method: Cross-Validated (10 fold, repeated 3 times) 
  # 
  # Resampling performance over subset size:
  #
  #  Variables Accuracy  Kappa AccuracySD KappaSD Selected
  #          1   0.9583 0.9375    0.06092 0.09139         
  #          2   0.9611 0.9417    0.06086 0.09129        *
  #          3   0.9583 0.9375    0.06092 0.09139         
  #          4   0.9583 0.9375    0.06092 0.09139         
  #
  # The top 2 variables (out of 2):
  #    Petal.Width, Petal.Length
  

• The rfe function in the caret package in R can be used with many different classifier methods (Full list), including:
  • svmLinear: Linear Support Vector Machine (SVM)
  • svmRadial: Radial Support Vector Machine (SVM)
  • knn: k-Nearest Neighbors (k-NN)
  • rpart: Recursive Partitioning (RPART)
  • glm: Generalized Linear Model (GLM)
  • glmnet: Lasso and Elastic-Net Regularized Generalized Linear Models (GLMNET)
  • xgbLinear: Extreme Gradient Boosting (XGBoost) with a linear objective
  • xgbTree: Extreme Gradient Boosting (XGBoost) with a tree-based objective
  • randomForest: Random Forest
  • gbm: Gradient Boosting Machines (GBM)
  • ctree: Conditional Inference Trees (CTREE)

Boruta

• Boruta Explained Exactly How You Wished Someone Explained to You
• Random Forest Feature Selection

knn

• How to code kNN algorithm in R from scratch
• Find K nearest neighbors, starting from a distance matrix
• K Nearest Neighbor : Step by Step Tutorial
  • Assumptions of KNN
  • KNN vs. K-mean
  • knn con
  • A small value of k will lead to a large variance in predictions. Alternatively, setting k to a large value may lead to a large model bias.
  • How to find best K value?

Random forest

• https://en.wikipedia.org/wiki/Random_forest
• randomForest package
  • Error: protect(): protection stack overflow. The trick works on my data with 26 obs and 28110 variables.
  • A complete guide to Random Forest in R
  • 8.5 Permutation Feature Importance from Interpretable Machine Learning by Christoph Molnar
  • 15 Variable Importance from caret
  • How to find the most important variables in R
  • ?importance
  • Random Forest Variable Importance and a simple example with 2 correlated predictors.
  • If you just print the importance object from the model they are the raw importance values. However, when you use the importance function, the default for the scale argument is TRUE which returns the importance values divided by the standard error.
  • Measures of variable importance in random forests
  • Difference between varImp (caret) and importance (randomForest) for Random Forest
  • Utilizing Machine Learning algorithms (GLMnet and Random Forest models) for Genomic Prediction of a Quantitative trait
• Gene selection and classification of microarray data using random forest 2006, and the R package varSelRF
  • The most reliable measure is based on the decrease of classification accuracy when values of a variable in a node of a tree are permuted randomly
  • this measure of variable importance is not the same as a non-parametric statistic of difference between groups, such as could be obtained with a Kruskal-Wallis test)
• Feature Importance in Random Forest

Gradient boost

• https://en.wikipedia.org/wiki/Gradient_boosting
• Gradient boosting in R

GBDT: Gradient Boosting Decision Trees

• An Introduction to Gradient Boosting Decision Trees
• GBDTs & Random Forests As Feature Transformers

caret

• https://cran.r-project.org/web/packages/caret/. In the vignette, the pls method was used (so the pls package needs to be installed in order to run the example).
  • The tune parameter tuneLength controls how many parameter values are evaluated,
  • The tuneGrid argument is used when specific values are desired.
• caret ebook
  • method = 'glmnet' only works for regression and classification problems; see train models by tag
  • The caret Package -> Model Training and Tuning
• Predictive Modeling with R and the caret Package (useR! 2013)
• Caret R Package for Applied Predictive Modeling

• coefficients from glmnet and caret are different for the same lambda.
  • the exact lambda you specified was not used by caret. the coefficients are interpolated from the coefficients actually calculated.
  • when you provide lambda to the glmnet call the model returns exact coefficients for that lambda

  library(caret)
  set.seed(0)
  train_control = trainControl(method = 'cv', number = 10)
  grid = 10 ^ seq(5, -2, length = 100)
  tune.grid = expand.grid(lambda = grid, alpha = 0)
  ridge.caret = train(x[train, ], y[train],
                      method = 'glmnet',
                      trControl = train_control,
                      tuneGrid = tune.grid)
  ridge.caret$bestTune
  
  ridge_full <- train(x, y,
                      method = 'glmnet',
                      trControl = trainControl(method = 'none'), 
                      tuneGrid = expand.grid(
                        lambda = ridge.caret$bestTune$lambda, alpha = 0)
                      )
  coef(ridge_full$finalModel, s = ridge.caret$bestTune$lambda)
  

• R caret train glmnet final model lambda values not as specified

• trainControl()
  • method - "boot", "boot632", "optimism_boot", "boot_all", "cv", "repeatedcv", "LOOCV", "LGOCV", ...
  • number - Either the number of folds or number of resampling iterations
  • repeats - For repeated k-fold cross-validation only: the number of complete sets of folds to compute
  • search - Either "grid" or "random", describing how the tuning parameter grid is determined.
  • seeds - an optional set of integers that will be used to set the seed at each resampling iteration. This is useful when the models are run in parallel.
• train(). A long printout comes from the foreach loop when the nominalTrainWorkflow() function is called.
  • method - A string specifying which classification or regression model to use. Some examples are knn, lm, nnet, rpart, glmboost, ... Possible values are found using names(getModelInfo()).
  • metric - ifelse(is.factor(y_dat), "Accuracy", "RMSE"). By default, possible values are "RMSE" and "Rsquared" for regression and "Accuracy" and "Kappa" for classification. If custom performance metrics are used (via the summaryFunction argument in trainControl, the value of metric should match one of the arguments.
  • maximize - ifelse(metric %in% c("RMSE", "logLoss", "MAE"), FALSE, TRUE)
• Return object from train()
  • results: matrix. Each row = one parameter (alpha, lambda)
  • finalModel$lambda: vector of very long length instead of what we specify. What is this?
  • finalModel$lambdaOpt
  • finalModel$tuneValue$alpha, finalModel$tuneValue$lambda
  • finalModel$a0
  • finalModel$beta

• Regularization: Ridge, Lasso and Elastic Net. Notice the repeatedcv method. If repeatedcv is selected, the performance RMSE is computed by the average of (# CV * # reps) RMSEs for each (alpha, lambda). Then the best tuning parameter (alpha, lambda) is selected with the minimum RMSE in performance. See the related code at Line 679, Line 744,Line 759 where trControl$selectionFunction()=best().

  library(glmnet)  # for ridge regression
  library(dplyr)   # for data cleaning
  
  data("mtcars")
  # Center y, X will be standardized in the modelling function
  y <- mtcars %>% select(mpg) %>% scale(center = TRUE, scale = FALSE) %>% as.matrix()
  X <- mtcars %>% select(-mpg) %>% as.matrix()
  dim(X)  # [1] 32 10
  
  library(caret)
  # Set training control
  train_control <- trainControl(method = "repeatedcv",
                                number = 3,    # 3-fold CV, 2 times
                                repeats = 2,
                                search = "grid",  # "random"
                                verboseIter = TRUE)
  
  # Train the model
  set.seed(123)    # seed for reproducibility
  elastic_net_model <- train(mpg ~ .,
                             data = cbind(y, X),
                             method = "glmnet",
                             preProcess = c("center", "scale"),
                             tuneLength = 4,  # 4 alphas x 4 lambdas
                             trControl = train_control)
  
  # Check multiple R-squared
  y_hat_enet <- predict(elastic_net_model, X)
  cor(y, y_hat_enet)^2
  
  names(elastic_net_model)
  #  [1] "method"       "modelInfo"    "modelType"    "results"      "pred"        
  #  [6] "bestTune"     "call"         "dots"         "metric"       "control"     
  # [11] "finalModel"   "preProcess"   "trainingData" "resample"     "resampledCM" 
  # [16] "perfNames"    "maximize"     "yLimits"      "times"        "levels"      
  # [21] "terms"        "coefnames"    "xlevels"
  elastic_net_model$bestTune   #  alpha and lambda
  
  names(elastic_net_model$finalModel)
  #  [1] "a0"          "beta"        "df"          "dim"         "lambda"     
  #  [6] "dev.ratio"   "nulldev"     "npasses"     "jerr"        "offset"     
  # [11] "call"        "nobs"        "lambdaOpt"   "xNames"      "problemType"
  # [16] "tuneValue"   "obsLevels"   "param"
  length(elastic_net_model$finalModel$lambda)
  # [1] 95
  length(elastic_net_model$finalModel$lambdaOpt)
  # [1] 1
  dim(elastic_net_model$finalModel$beta)
  # [1] 10 95
  which(elastic_net_model$finalModel$lambda == elastic_net_model$finalModel$lambdaOpt)
  # integer(0)      <=== Weird, why
  

• For the Repeated k-fold Cross Validation. The final model accuracy is taken as the mean from the number of repeats.
• Penalized Regression Essentials: Ridge, Lasso & Elastic Net -> Using caret package. tuneLength was used. Note that the results element gives the details for each candidate of parameter combination.

  # Load the data
  data("Boston", package = "MASS")
  # Split the data into training and test set
  set.seed(123)
  training.samples <- Boston$medv %>%
    createDataPartition(p = 0.8, list = FALSE)
  train.data  <- Boston[training.samples, ]
  
  set.seed(123)  # affect CV samples, not affect tuning parameter set
  cv_5 <- trainControl(method = "cv", number = 5) # number of folds
  model <- train(
    medv ~., data = train.data, method = "glmnet",
    trControl = cv_5,
    tuneLength = 10  # 10 alphas x 10 lambdas
  )
  model
  # one row per parameter combination
  model$results
  # Best tuning parameter
  model$bestTune
  

• TuneGrid and TuneLength in Caret
• Exploring the caret package

Compare models

• Bland–Altman plot. It appears in the caret vignette. If we have paired data (e.g. one represents the AUC from the model A, and the other represents the AUC from the model B), we draw a scatterplot that X-axis = average, Y-axis = difference. M vs A plot.