Machine Learning models for Breast Cancer disease prediction
Badawi Amin Muhammed
2025-11-23
Introduction
This article demostrates my work on building machine-learning models to predict the course of different diseases. It moves from building a model, evaluating its performance, and answering or addressing different disease related questions using machine learning.
Setup
All analyses are done in R using RStudio. some packages might not install and load due to CRAN conflicting changing versions and systems comapatibilty.
All figures are produced with ggplot2.
The dataset used for analyses in this project is the Breast Cancer Wisconsin (Diagnostic) Dataset. The data was downloaded from the UC Irvine Machine Learning Repository (in three categories).
note: use python for easier scaping of the repository.
The first dataset looks at the predictor classes:
- malignant or
- benign breast mass.
The features characterize cell nucleus properties and were generated from image analysis of fine needle aspirates (FNA) of breast masses:
-Sample ID (code number)
-Clump thickness
-Uniformity of cell size
-Uniformity of cell shape
-Marginal adhesion
-Single epithelial cell size
-Number of bare nuclei
-Bland chromatin
-Number of normal nuclei
-Mitosis
-Classes, i.e. diagnosis
bc_data <- read.csv("C:/Users/hp/Documents/breast+cancer+wisconsin+diagnostic/winscosin_breast_cancer.csv",
header =TRUE,
sep = ",")
colnames(bc_data) <- c("sample_code_number",
"radius",
"clump_thickness",
"uniformity_of_cell_size",
"uniformity_of_cell_shape",
"marginal_adhesion",
"single_epithelial_cell_size",
"bare_nuclei",
"bland_chromatin",
"normal_nucleoli",
"mitosis",
"classes")
bc_data$classes <- ifelse(bc_data$classes == "B", "benign",
ifelse(bc_data$classes == "M", "malignant", NA))
bc_data[bc_data == "?"] <- NA
# how many NAs are in the data
length(which(is.na(bc_data)))## [1] 0# how many samples would we loose, if we removed them?
nrow(bc_data)## [1] 569in case of Missing values, they can be imputed with the
mice package.
# impute missing data
#library(mice)
#bc_data[,2:11] <- apply(bc_data[, 2:11], 2, function(x) as.numeric(as.character(x)))
#dataset_impute <- mice(bc_data[, 2:11], print = FALSE)
#bc_data <- cbind(bc_data[, 11, drop = FALSE], mice::complete(dataset_impute, 1))
bc_data$classes <- as.factor(bc_data$classes)
# how many benign and malignant cases are there?
summary(bc_data$classes)## benign malignant
## 357 212Data exploration
-Response variable for classification
library(ggplot2)
ggplot(bc_data, aes(x = classes, fill = classes)) +
geom_bar()
target classification
Response variable for regression
ggplot(bc_data, aes(x = clump_thickness)) +
geom_histogram(bins = 9)
distribution plot of responses
Principal Component Analysis
library(ggplot2)
library(ellipse)
# --- 1. Prepare predictors & class labels ---
# Assumption: rows = samples, predictors in cols 2:11, classes column named 'classes'
predictors <- bc_data[, 2:11]
# ensure predictors are numeric
predictors <- as.data.frame(lapply(predictors, function(x) {
if(is.factor(x)) as.numeric(as.character(x)) else as.numeric(x)
}))
groups <- factor(bc_data$classes) # ensure factor
# quick sanity checks
if(nrow(predictors) != length(groups)) stop("Number of samples in predictors and classes differ.")
# --- 2. PCA (no transpose) ---
pca_res <- prcomp(predictors, center = TRUE, scale. = TRUE)
# proportion of variance
var_exp <- (pca_res$sdev^2) / sum(pca_res$sdev^2)
# --- 3. Prepare scores dataframe (PC1, PC2) and attach groups ---
pca_scores <- as.data.frame(pca_res$x[, 1:2])
colnames(pca_scores) <- c("PC1", "PC2")
pca_scores$groups <- groups
# --- 4. Centroids by group ---
centroids <- aggregate(cbind(PC1, PC2) ~ groups, pca_scores, mean)
# --- 5. Confidence ellipses per group (95%) ---
conf.rgn <- do.call(rbind, lapply(levels(pca_scores$groups), function(g) {
sub <- pca_scores[pca_scores$groups == g, c("PC1", "PC2")]
# ellipse needs at least 3 points to compute a sensible covariance; skip otherwise
if (nrow(sub) < 3) return(NULL)
ell <- ellipse::ellipse(cov(sub), centre = as.numeric(centroids[centroids$groups == g, 2:3]), level = 0.95)
data.frame(groups = g, PC1 = ell[, 1], PC2 = ell[, 2], stringsAsFactors = FALSE)
}))
# Keep factor levels consistent
if (!is.null(conf.rgn)) conf.rgn$groups <- factor(conf.rgn$groups, levels = levels(pca_scores$groups))then we plot
Plot of PCA
ggplot(data = pca_scores, aes(x = PC1, y = PC2, group = groups, color = groups)) +
{ if (!is.null(conf.rgn)) geom_polygon(data = conf.rgn, aes(x = PC1, y = PC2, fill = groups), alpha = 0.2) } +
geom_point(size = 2, alpha = 0.7) +
scale_color_brewer(palette = "Set1") +
labs(
color = "",
fill = "",
x = paste0("PC1: ", round(var_exp[1] * 100, 2), "% variance"),
y = paste0("PC2: ", round(var_exp[2] * 100, 2), "% variance")
) +
theme_minimal()
plot of PCA on Classes
Features
library(tidyr)
gather(bc_data, x, y, clump_thickness:mitosis) %>%
ggplot(aes(x = y, color = classes, fill = classes)) +
geom_density(alpha = 0.3) +
facet_wrap( ~ x, scales = "free", ncol = 3)
plot of feautures
Machine Learning packages for R
the doParallel is one the most used packages for
clustering configuration
# configure multicore
library(doParallel)
cl <- makeCluster(detectCores())
registerDoParallel(cl)
library(caret)Training, validation and test data we train, test and validate our model as follows:
set.seed(42)
index <- createDataPartition(bc_data$classes, p = 0.7, list = FALSE)
train_data <- bc_data[index, ]
test_data <- bc_data[-index, ]Load package dplyr
library(dplyr)
rbind(data.frame(group = "train", train_data),
data.frame(group = "test", test_data)) %>%
gather(x, y, clump_thickness:mitosis) %>%
ggplot(aes(x = y, color = group, fill = group)) +
geom_density(alpha = 0.3) +
facet_wrap( ~ x, scales = "free", ncol = 3)
Regression
now we define our model:
set.seed(42)
model_glm <- caret::train(clump_thickness ~ .,
data = train_data,
method = "glm",
preProcess = c("scale", "center"),
trControl = trainControl(method = "repeatedcv",
number = 10,
repeats = 10,
savePredictions = TRUE,
verboseIter = FALSE))Model Results
model_glm## Generalized Linear Model
##
## 399 samples
## 11 predictor
##
## Pre-processing: scaled (11), centered (11)
## Resampling: Cross-Validated (10 fold, repeated 10 times)
## Summary of sample sizes: 359, 359, 359, 359, 359, 359, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 3.863022 0.2788419 3.019432predict on test data:
predictions <- predict(model_glm, test_data)precitors vs resduals
# model_glm$finalModel$linear.predictors == model_glm$finalModel$fitted.values
data.frame(residuals = resid(model_glm),
predictors = model_glm$finalModel$linear.predictors) %>%
ggplot(aes(x = predictors, y = residuals)) +
geom_jitter() +
geom_smooth(method = "lm")
model prediction
plot of residuals:
# y == train_data$clump_thickness
data.frame(residuals = resid(model_glm),
y = model_glm$finalModel$y) %>%
ggplot(aes(x = y, y = residuals)) +
geom_jitter() +
geom_smooth(method = "lm")
plot of resuduals
Actual vs Predicted:
data.frame(actual = test_data$clump_thickness,
predicted = predictions) %>%
ggplot(aes(x = actual, y = predicted)) +
geom_jitter() +
geom_smooth(method = "lm")
plot of actual vs predicted
Classification
Decision trees
using the rpart package:
library(rpart)
library(rpart.plot)
set.seed(42)
fit <- rpart(classes ~ .,
data = train_data,
method = "class",
control = rpart.control(xval = 10,
minbucket = 2,
cp = 0),
parms = list(split = "information"))
rpart.plot(fit, extra = 100)
decision tree
Random Forests
Random Forests predictions are based on the generation of multiple classification trees. They can be used for both, classification and regression tasks. Here, I show a classification task:
set.seed(42)
model_rf <- caret::train(classes ~ .,
data = train_data,
method = "rf",
preProcess = c("scale", "center"),
trControl = trainControl(method = "repeatedcv",
number = 10,
repeats = 10,
savePredictions = TRUE,
verboseIter = FALSE))When we specify savePredictions = TRUE, we can access
the cross-validation results with model_rf$pred.
model_rf$finalModel$confusion## benign malignant class.error
## benign 238 12 0.04800000
## malignant 10 139 0.06711409-Feature Importance
imp <- model_rf$finalModel$importance
imp[order(imp, decreasing = TRUE), ]## bland_chromatin uniformity_of_cell_size
## 39.715577 28.518652
## bare_nuclei uniformity_of_cell_shape
## 27.623267 24.423469
## radius single_epithelial_cell_size
## 21.823507 14.032896
## clump_thickness sample_code_number
## 9.876583 6.157004
## marginal_adhesion mitosis
## 5.658340 4.377050
## normal_nucleoli
## 3.973324estimating variable importance and plot:
# estimate variable importance
importance <- varImp(model_rf, scale = TRUE)
plot(importance)
plot of variable importance
predicting test data:
confusionMatrix(predict(model_rf, test_data), test_data$classes)## Confusion Matrix and Statistics
##
## Reference
## Prediction benign malignant
## benign 100 7
## malignant 7 56
##
## Accuracy : 0.9176
## 95% CI : (0.8657, 0.9542)
## No Information Rate : 0.6294
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8235
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9346
## Specificity : 0.8889
## Pos Pred Value : 0.9346
## Neg Pred Value : 0.8889
## Prevalence : 0.6294
## Detection Rate : 0.5882
## Detection Prevalence : 0.6294
## Balanced Accuracy : 0.9117
##
## 'Positive' Class : benign
## result of predictions:
results <- data.frame(actual = test_data$classes,
predict(model_rf, test_data, type = "prob"))
results$prediction <- ifelse(results$benign > 0.5, "benign",
ifelse(results$malignant > 0.5, "malignant", NA))
results$correct <- ifelse(results$actual == results$prediction, TRUE, FALSE)
ggplot(results, aes(x = prediction, fill = correct)) +
geom_bar(position = "dodge")
plot of prediction against benign class:
#prediction against benign class
ggplot(results, aes(x = prediction, y = benign, color = correct, shape = correct)) +
geom_jitter(size = 3, alpha = 0.6)
Feature Selection
Performing feature selection on the whole dataset would lead to prediction bias, we therefore need to run the whole modeling process on the training data alone!
Correlations between all features are calculated and visualised with
the corrplot package. we are then removing all features
with a correlation higher than 0.7, keeping the feature with the lower
mean.
library(corrplot)
# calculate correlation matrix
corMatMy <- cor(train_data[2:11])
corrplot(corMatMy, order = "hclust")
correlation filter
#Apply correlation filter at 0.70,
highlyCor <- colnames(train_data[, -1])[findCorrelation(corMatMy, cutoff = 0.7, verbose = TRUE)]## Compare row 8 and column 7 with corr 0.939
## Means: 0.655 vs 0.479 so flagging column 8
## Compare row 7 and column 6 with corr 0.896
## Means: 0.605 vs 0.439 so flagging column 7
## Compare row 3 and column 1 with corr 0.998
## Means: 0.519 vs 0.399 so flagging column 3
## Compare row 1 and column 4 with corr 0.987
## Means: 0.427 vs 0.37 so flagging column 1
## All correlations <= 0.7# which variables are flagged for removal?
highlyCor## [1] "bland_chromatin" "bare_nuclei"
## [3] "uniformity_of_cell_size" "radius"remove highly correlated variables:
#then we remove these variables
train_data_cor <- train_data[, which(!colnames(train_data) %in% highlyCor)]
set.seed(7)Genetic Algorithm (GA)
The Genetic Algorithm (GA) has been developed based on evolutionary principles of natural selection: It aims to optimize a population of individuals with a given set of genotypes by modeling selection over time. In each generation (i.e. iteration), each individual’s fitness is calculated based on their genotypes. Then, the fittest individuals are chosen to produce the next generation. This subsequent generation of individuals will have genotypes resulting from (re-) combinations of the parental alleles. These new genotypes will again determine each individual’s fitness. This selection process is iterated for a specified number of generations and (ideally) leads to fixation of the fittest alleles in the gene pool.
This concept of optimization can be applied to non-evolutionary models as well, like feature selection processes in machine learning.
library(tidymodels)
library(GA)
set.seed(27)
# Prepare data
bc_data <- bc_data %>%
mutate(classes = factor(classes, levels = c("malignant", "benign")))
set.seed(27)
folds <- vfold_cv(bc_data, v = 5, strata = classes)
ga_fitness <- function(bitstring, data, folds) {
# Which features were selected?
selected <- which(bitstring == 1)
# If GA selects 0 features → penalise it
if (length(selected) == 0) return(0)
# Subset data
df <- data %>% select(all_of(selected), classes)
# Define model
mod <- rand_forest(trees = 500) %>%
set_engine("randomForest") %>%
set_mode("classification")
# CV performance
rs <- fit_resamples(
mod,
classes ~ .,
resamples = folds,
metrics = metric_set(accuracy),
control = control_resamples(verbose = FALSE)
)
# Return mean accuracy
collect_metrics(rs) %>%
filter(.metric == "accuracy") %>%
pull(mean)
}
set.seed(27)
ga_result <- ga(
type = "binary",
fitness = function(bitstring) ga_fitness(bitstring, bc_data, folds),
nBits = ncol(bc_data) - 1, # number of predictors
popSize = 20,
maxiter = 30,
pmutation = 0.2,
run = 10
)
selected_features <- which(ga_result@solution[1,] == 1)
selected_features_names <- names(bc_data)[selected_features + 1] # +1 offset
selected_features_names## [1] "radius" "uniformity_of_cell_shape"
## [3] "marginal_adhesion" "bare_nuclei"
## [5] "normal_nucleoli" "classes"final_data <- bc_data %>% select(all_of(selected_features_names), classes)
rf_model <- rand_forest(trees = 500) %>%
set_engine("randomForest") %>%
set_mode("classification") %>%
fit(classes ~ ., data = final_data)plot of MEAN fitness by generation:
plot(ga_result)
Grid Search
The R package h2o was supposed to be used beacause of its features, however, due to version constriant on the latest version of R, we had to resort to:
Use tidymodels + ranger + dials + tune + workflows
This gives us:
✔ Hyperparameter grid search (better than H2O Grid)
✔ Feature importance and partial dependence (better plots)
✔ Data reshaping with tidyr (better than reshape2)
✔ Classification with ranger (fast implementation of Random Forest)
reshape2 is also outdated, tidyr replaces it entirely.
library(tidymodels)
library(ranger)
library(ggplot2)
library(tidyr)
#Step 2 — Split data
# dataset is bc_data and class column is classes:
set.seed(27)
data_split <- initial_split(bc_data, prop = 0.8, strata = classes)
train_data <- training(data_split)
test_data <- testing(data_split)
#Step 3 — Preprocess (recipes)
rec <- recipe(classes ~ ., data = train_data) %>%
step_normalize(all_predictors()) %>%
step_zv(all_predictors())
#Step 4 — Specify Random Forest model (ranger)
rf_model <- rand_forest(
mtry = tune(),
min_n = tune(),
trees = tune()
) %>%
set_engine("ranger", importance = "impurity") %>%
set_mode("classification")
#Step 5 — Define workflow
wf <- workflow() %>%
add_model(rf_model) %>%
add_recipe(rec)
set.seed(27)
#Step 6 — Hyperparameter Grid Search (H2O alternative)
rf_grid <- grid_regular(
mtry(range = c(2, 10)),
min_n(range = c(3, 20)),
trees(range = c(200, 400)),
levels = 3
)Cross-validation:
cv <- vfold_cv(train_data, v = 5, strata = classes)
tuned_rf <- tune_grid(
wf,
resamples = cv,
grid = rf_grid,
metrics = metric_set(roc_auc, accuracy)
)
cv## # 5-fold cross-validation using stratification
## # A tibble: 5 × 2
## splits id
## <list> <chr>
## 1 <split [363/91]> Fold1
## 2 <split [363/91]> Fold2
## 3 <split [363/91]> Fold3
## 4 <split [363/91]> Fold4
## 5 <split [364/90]> Fold5Selection of best model
Select the highest performing model according to ROC-AUC.
Why it matters:
• Guarantees maximum diagnostic power
• Eliminates underperforming configurations
best_params <- select_best(tuned_rf, metric="roc_auc")
best_params## # A tibble: 1 × 4
## mtry trees min_n .config
## <int> <int> <int> <chr>
## 1 10 400 3 pre0_mod25_post0fiitng of final model
final_rf <- finalize_workflow(wf, best_params)
#Fit final model
rf_fit <- fit(final_rf, train_data)
rf_fit## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 2 Recipe Steps
##
## • step_normalize()
## • step_zv()
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Ranger result
##
## Call:
## ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~10L, x), num.trees = ~400L, min.node.size = min_rows(~3L, x), importance = ~"impurity", num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE)
##
## Type: Probability estimation
## Number of trees: 400
## Sample size: 454
## Number of independent variables: 11
## Mtry: 10
## Target node size: 3
## Variable importance mode: impurity
## Splitrule: gini
## OOB prediction error (Brier s.): 0.04072395plotting feauture importance
#Plot Feature Importance (H2O alternative)
imp <- rf_fit %>%
extract_fit_engine() %>%
ranger::importance()
imp_df <- data.frame(
feature = names(imp),
importance = imp
)
ggplot(imp_df, aes(x = reorder(feature, importance), y = importance)) +
geom_col() +
coord_flip() +
labs(x = "Feature", y = "Importance")
Melt and Plot (reshape2 alternative)
melted <- bc_data %>%
pivot_longer(
cols = -classes,
names_to = "feature",
values_to = "value"
)plotting of melted feautures
ggplot(melted, aes(x = feature, y = value, color = classes)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90))
plotting of melted
Evaluation of Test Set
pred <- predict(rf_fit, test_data, type = "prob")
roc_auc_vec(test_data$classes, pred$.pred_malignant)## [1] 0.9896641with a ROC_AUC of 0.98, we can scientifically infer that the
modelsignals very good feature quality, good model architecture and
fit.
hence, our classifier can achieve high true positive while
maintainig a low false positive rate.
Conclusion
— Significance of the Modelling Approach
The workflow presented combines statistical dimensionality reduction, evolutionary optimisation, and modern machine learning engineering to achieve a high-performing and interpretable predictive system for disease classification. It demonstrates an intelligent and transparent approach to model development, ensuring that:
✔ The data are cleaned, normalised, and scientifically treated
✔ Signal-to-noise ratios are improved through PCA and GA
✔ Model hyperparameters are optimised through systematic search
✔ Predictive reliability is validated using stratified resampling
✔ Results are interpretable to clinicians, researchers, and stakeholders
The integration of GA-optimised features and tuned Random Forest models provides an innovative analytical framework that exceeds traditional modelling pipelines. Importantly, by leveraging tidymodels—as a replacement for deprecated or unstable packages like H2O and pcaGoPromoter—the final workflow becomes far more reproducible, robust, and suitable for large-scale clinical or epidemiological deployments.
This pipeline not only predicts disease risk but also reveals the biological patterns that distinguish patient groups, offering practical relevance for diagnostics, triage, and public-health surveillance.
Badawi Aminu Muhammed
08065440075
official.badawy@gmail.com