This tutorial guides you through implementing the two-stage sampling aggregation (\(\mbox{TSMA}\)), designed to handle unbalanced datasets. The \(\mbox{TSMA}\) framework incorporates bagging concepts, resampling techniques, and machine learning techniques aimed at significantly improving the accuracy of predictive models.
Prerequisites: R environment with necessary libraries:
randomForest, caret, ROCR,
tidyverse, mltools.
# Remove all objects from the current workspace to start with a clean environment
rm(list = ls())
# Load necessary libraries for the tutorial
library(caret) # Load the caret package for data splitting, training, and performance evaluation
library(ROSE) # Load the ROSE package for undersampling
library(ROCR) # Load the ROCR package for ROC analysis and performance metrics
library(tidyverse) # Load the tidyverse package for data manipulation and visualization
library(mltools) # Load the mltools package for machine learning tools and utilities
# Define a function 'metrics' to calculate various evaluation metrics
metrics <- function(pred.class, pred.prob, ref) {
# pred.class: Predicted class labels
# pred.prob: Predicted probabilities for the positive class
# ref: Actual class labels (reference)
# Calculate the confusion matrix using the caret package
cm <- caret::confusionMatrix(data = pred.class, reference = ref, positive = '1')
# Extract accuracy from the confusion matrix
accuracy = cm$overall[['Accuracy']]
# Extract precision from the confusion matrix
precision = cm$byClass[['Precision']]
# Extract recall from the confusion matrix
recall = cm$byClass[['Recall']]
# Extract specificity from the confusion matrix
specificity = cm$byClass[['Specificity']]
# Calculate the balanced accuracy
balanced_acc = cm$byClass[['Balanced Accuracy']]
# Calculate the F1 measure
f1_measure = cm$byClass[['F1']]
# Create a prediction object for ROC analysis
pred_model <- prediction(pred.prob, ref)
# Calculatevthe Area Under the Precision-Recall Curve (AUCPR)
aucpr <- performance(pred_model, measure = "aucpr")@y.values[[1]]
# Calculate Cohen's Kappa
kappa <- cm$overall[['Kappa']]
# Calculate and round the Area Under the Receiver Operating Characteristic (AUROC)
auroc <- performance(pred_model, measure = "auc")@y.values[[1]]
# Calculate Matthews Correlation Coefficient (MCC)
mcc <- mcc(pred.class, ref)
# Calculate the geometric mean of recall and specificity
g_mean1 <- sqrt(recall * specificity)
# Calculate the geometric mean of recall and precision
g_mean2 <- sqrt(recall * precision)
# Compile all the calculated metrics into a data frame
output <- data.frame(accuracy = accuracy, precision = precision, recall = recall,
specificity = specificity, balanced_acc = balanced_acc,
f1_measure = f1_measure, aucpr = aucpr, kappa = kappa,
auroc = auroc, mcc = mcc, g_mean1 = g_mean1, g_mean2 = g_mean2)
# Return the data frame with all metrics
return(output)
}# This will read the 'unbalanced_disease' data into the environment
unbalanced_disease <- readRDS("~/Documents/camrsa/Data/unbalanced_disease.rds")
# Display the first 5 rows and 6 columns of the loaded dataset for a quick overview
unbalanced_disease[1:12, 1:7]## case_control sex snps_1 snps_2 snps_3 snps_4 snps_5
## 1 0 Female 2 2 2 2 2
## 2 0 Male 2 2 2 2 2
## 3 0 Female 2 2 2 2 2
## 4 0 Male 2 2 2 2 2
## 5 0 Male 2 2 2 2 2
## 6 0 Female 2 2 2 2 2
## 7 0 Male 2 2 2 2 2
## 8 0 Male 2 2 2 2 2
## 9 0 Male 2 2 2 2 2
## 10 0 Male 2 2 2 2 2
## 11 0 Female 2 2 2 2 2
## 12 1 Male 2 2 2 2 2# Create a frequency table for the 'case_control' column to understand its distribution
table(unbalanced_disease$case_control)##
## 0 1
## 1527 55The example data is an unbalanced dataset. Within the dataset, case_control represents the outcome variable, which includes 55 cases and 1,527 controls, resulting in an unbalance rate of 3%. Meanwhile, sex is the phenotypic factor that, along with 1500 genotypes, is used as a predictor.
Divide your dataset: Split your unbalanced data into training and testing sets, maintaining the ratio of major to minor classes as in the original dataset.
# Set seed for reproducibility
set.seed(123)
# Create indices for a training set containing 80% of the data, ensuring the proportion of 'case_control' is maintained
train.rows <- createDataPartition(y = unbalanced_disease$case_control, p = 0.8, list = FALSE)
# Subset the dataset to create a training set using the indices, comprising 80% of the original data
unbalanced_disease_train <- unbalanced_disease[train.rows,]
# Create a testing set with the remaining 20% of the data not included in the training set
unbalanced_disease_test <- unbalanced_disease[-train.rows,]# Train an Elastic Net model
elastic <- train(
case_control ~ ., data = unbalanced_disease_train, method = "glmnet",
trControl = trainControl("cv", number = 5), # Set up cross-validation with 5 folds
tuneLength = 10) # Specify the length of the tuning parameter grid
# Use the trained Elastic Net model to make predictions on the test dataset
predictions <- predict(elastic, unbalanced_disease_test)
# Predict the probabilities of the positive class for the test dataset using the Elastic Net model
predict <- predict(elastic, unbalanced_disease_test, type = "prob")[,2]
# Calculate various performance metrics (like accuracy, precision, recall, etc.)
# for the Elastic Net model using the predicted classes and probabilities
el <- metrics(predictions, predict, unbalanced_disease_test$case_control)# Apply under-sampling to balance the dataset, ensuring equal representation of each class in 'case_control'
under_caret_el <- ovun.sample(case_control ~ ., data = unbalanced_disease_train , method = "both", p=0.5)$data
# Train an Elastic Net model using the balanced dataset
# - Specify the model formula where 'case_control' is the response variable and all other columns are predictors
# - Set cross-validation control using 'cv' with 5 folds for model validation
# - Define the length of the grid for hyperparameter tuning (tuneLength = 10)
elastic_under_el <- train(
case_control ~ ., data = under_caret_el, method = "glmnet",
trControl = trainControl("cv", number = 5),
tuneLength = 10)
# Use the trained Elastic Net model with under-sampling to make predictions on the test set
predictions_el <- elastic_under_el %>% predict(unbalanced_disease_test)
# Predict the probabilities using the Elastic Net model and select the probabilities of the positive class
predict_el <- predict(elastic_under_el, unbalanced_disease_test, type = "prob")[,2]
# Evaluate the model using various metrics
un_el_downcaret <- metrics(predictions_el, predict_el, unbalanced_disease_test$case_control)- Stage 1 Bootstrap Samples
- Stage 2 Resampling
- Model Training, Aggregation and Prediction
#Prepare for the Aggregation step
# Initialize a matrix to store the estimated probabilities for each bootstrap sample
estimated_probabilities_tsma <- matrix(NA, nrow = nrow(unbalanced_disease_test), ncol = 100)
# Stage 1: Iterate over 100 bootstrap samples
for(j in 1:100){
# Stage 1: Bootstrap Samples
# Create a bootstrap sample from the training dataset (90% of the training data)
train_sample <- unbalanced_disease_train[sample(nrow(unbalanced_disease_train),
size = round(0.90 * nrow(unbalanced_disease_train))), ]
# Stage 2: Resampling
# Apply under-sampling to balance the bootstrap sample
under_caret <- ovun.sample(case_control ~ ., data = train_sample , method = "both", p=0.5)$data
# Model Training: Train an Elastic Net model on the downsampled data
elastic_under_tsma <- train(
case_control ~ ., data = under_caret, method = "glmnet",
trControl = trainControl("cv", number = 5),
tuneLength = 10)
# Predict probabilities for the test set using the trained model
predicted_tsma <- predict(elastic_under_tsma, unbalanced_disease_test, type = "prob")[,2]
# Store the predicted probabilities in the matrix
estimated_probabilities_tsma[, j] <- predicted_tsma
}
# Aggregation and Prediction
# Calculate the mean of the estimated probabilities for each observation
pred.prob_tsma <- rowMeans(estimated_probabilities_tsma, na.rm = TRUE)
# Convert the mean probabilities to class predictions based on a threshold of 0.5
pred.class_tsma <- as.factor(ifelse(pred.prob_tsma > 0.5, 1, 0))
# Evaluate the tsma model using various metrics
un_el_tsma <- metrics(pred.class_tsma, pred.prob_tsma, unbalanced_disease_test$case_control)# Combine the results from the two models (Elastic Net with undersampling and TSMA-EN) into a single data frame
compare_result <- rbind(el, un_el_downcaret, un_el_tsma)
# Rename the columns of the resulting data frame for clarity
colnames(compare_result) = c("Accuracy", "Precision", "Recall", "Specificity", "Balanced Accuracy",
"F1 score", "AUC-PR", "Cohen's Kappa", "AUC-ROC", "MCC", "G Mean1", "G Mean2")
# Assign descriptive row names to identify each model's results
rownames(compare_result) = c("EN Alone", "EN with Under", "TSMA-EN")
# Display the comparison of performance metrics between the two models
compare_result## Accuracy Precision Recall Specificity Balanced Accuracy
## EN Alone 0.9651899 NA 0.0000000 1.0000000 0.5000000
## EN with Under 0.9505063 0.76 0.8181818 0.9508197 0.9074516
## TSMA-EN 0.9905063 0.90 1.0000000 0.9967213 0.9983607
## F1 score AUC-PR Cohen's Kappa AUC-ROC MCC G Mean1
## EN Alone NA 0.9648164 0.0000000 0.9976155 0.0000000 0.00000
## EN with Under 0.8571429 0.9103142 0.8522444 0.9949329 0.8532799 0.90305
## TSMA-EN 0.9583607 0.9649679 0.9548829 0.9973174 0.9985359 0.99835
## G Mean2
## EN Alone NA
## EN with Under 0.8581163
## TSMA-EN 0.9574271From the results, we observe a comparison between three models: EN Alone, EN combined with UNDER and our proposed model \(\mbox{TSMA}_{\text{EN}}\), which incorporates EN combined with UNDER into our two-stage sampling approach. Based on the metrics provided, \(\mbox{TSMA}_{\text{EN}}\) demonstrates significant improvements across various performance indicators.
These improvements, as shown in the results section, suggest that our proposed \(\mbox{TSMA}_{\text{EN}}\) model is more accurate and reliable, with significant enhancements in precision, F1 score, Cohen’s Kappa, MCC, G Mean2 and other critical metrics that are important for model evaluation.
By adopting the TSMA system, you can implement the TSMA method in R to handle unbalanced datasets effectively. This approach is particularly beneficial in scenarios where class imbalance significantly affects the performance of predictive models.
The example data utilized for the tutorial is located in the Data folder.