Comparing modelr and CARET

Introduction

For prediction problems, moving from individual packages like modelr and tidymodels to a unified framework like caret (Classification And REgression Training) is a massive quality-of-life upgrade. While modelr is great for tidyverse-aligned EDA, it forces you to manually orchestrate your own cross-validation loops, data transformations, and model syntax across different underlying packages.

caret standardizes this into a single interface. This notebook demonstrates the differences using the Boston housing problem (ISLR, Chapter #4, Exercise #16) using both approaches so you can see the exact code delta.

Preparing the Data

First, we need to load the data and create our binary response variable (crime_above_median) based on the median of crim.

library(tidyverse)
library(ISLR2) # Contains the Boston dataset
library(modelr)
library(caret)
library(MASS)  # For LDA
library(e1071) # For Naive Bayes
library(class) # For KNN
library(DT)    # For interactive tables
library(plotly)# For interactive plots

# Prepare the data
boston_clean <- Boston %>%
  mutate(
    crime_above_median = factor(ifelse(crim > median(crim), "Yes", "No"),
                                levels = c("Yes","No"))
  ) %>%
  dplyr::select(-crim) # Remove original continuous variable to prevent data leakage

# We'll use a subset of predictors for illustration: indus, nox, rm, age, dis
predictors_subset <- c("indus", "nox", "rm", "age", "dis", "crime_above_median")
df_subset <- boston_clean %>% dplyr::select(all_of(predictors_subset))

datatable(df_subset)

The modelr Way (Manual Orchestration)

To get honest out-of-sample error using modelr, we have to manually create cross-validation folds, write custom wrapper functions to handle the vastly different syntax of each algorithm, map the models over the folds, extract predictions, and compute the error.

set.seed(42)

# 1. Create 10-fold CV structure
cv_data <- crossv_kfold(df_subset, k = 10)

# 2. Write custom fitting/predicting functions because syntax varies wildly
fit_predict_cv <- function(train_df, test_df) {
  # Logistic Regression
  glm_mod <- glm(crime_above_median ~ ., data = train_df, family = binomial)
  glm_pred <- ifelse(predict(glm_mod, test_df, type = "response") > 0.5, "Yes", "No")
  
  # LDA
  lda_mod <- lda(crime_above_median ~ ., data = train_df)
  lda_pred <- predict(lda_mod, test_df)$class
  
  # Naive Bayes
  nb_mod <- naiveBayes(crime_above_median ~ ., data = train_df)
  nb_pred <- predict(nb_mod, test_df)
  
  # KNN (Requires matrix extraction and scale normalization inside the loop!)
  train_x <- scale(train_df[, 1:5])
  test_x <- scale(test_df[, 1:5], center = attr(train_x, "scaled:center"), scale = attr(train_x, "scaled:scale"))
  knn_pred <- knn(train = train_x, test = test_x, cl = train_df$crime_above_median, k = 5)
  
  # True values
  actual <- test_df$crime_above_median
  
  # Return classification errors
  tibble(
    logistic = mean(glm_pred != actual),
    lda      = mean(lda_pred != actual),
    nb       = mean(nb_pred != actual),
    knn      = mean(knn_pred != actual)
  )
}

# 3. Map the function over all 10 folds and calculate mean error
modelr_results <- cv_data %>%
  mutate(
    train_df = map(train, as_tibble),
    test_df  = map(test, as_tibble),
    errors   = map2(train_df, test_df, fit_predict_cv)
  ) %>%
  unnest(errors) %>%
  summarise(across(logistic:knn, mean)) %>%
   pivot_longer(
    cols = everything(), 
    names_to = "Model", 
    values_to = "Error"
  ) %>%
  arrange(Error) # Sorts from lowest error to highest error

datatable(modelr_results)

Pain Points with modelr:

• Syntax Fragmentation: glm() uses type = "response", lda() returns a list where you need $class, and knn() doesn’t even use standard formula syntax (it requires separate data matrices).
  
• Pre-processing Leakage: Pre-processing like scaling for KNN must be done inside each CV loop to avoid data leakage, adding significant boilerplate.
  
• Bookkeeping: You are entirely responsible for collecting vectors, calculating percentages, and tracking results.

The caret Way (Unified Streamlining)

caret solves this by introducing a single unified interface (train()) and a control structure (trainControl()). It standardizes data splitting, pre-processing, hyperparameter tuning, and variable evaluation.

set.seed(42)

# 1. Define the CV strategy once
cv_control <- trainControl(
  method = "cv", 
  number = 10,
  classProbs = TRUE # Required for probabilistic models
)

# 2. Fit all models using identical syntax
model_glm <- train(crime_above_median ~ ., data = df_subset, method = "glm", family = "binomial", trControl = cv_control)
model_lda <- train(crime_above_median ~ ., data = df_subset, method = "lda", trControl = cv_control)
model_nb  <- train(crime_above_median ~ ., data = df_subset, method = "nb", trControl = cv_control)

# KNN gets automatic scaling/centering and hyperparameter tuning natively!
model_knn <- train(crime_above_median ~ ., data = df_subset, method = "knn", 
                   preProcess = c("center", "scale"), 
                   tuneGrid = expand.grid(k = c(1, 3, 5, 7, 9)),
                   trControl = cv_control)

# 3. Compare out-of-sample metrics instantly
results <- resamples(list(
  Logistic = model_glm, 
  LDA = model_lda, 
  NaiveBayes = model_nb, 
  KNN = model_knn
))

# Visual comparison of model metrics across folds
bwplot(results, metric = c("Accuracy", "Kappa"), main = "Model Performance Comparison")

The bwplot() function creates a side-by-side boxplot. This lets you see the distribution of Accuracy and Kappa values across all 10 folds for each model, highlighting stability and spread.

Do not look just at the average. Look at the boxplot width:

• Narrow boxplots mean the model is stable and predictable across different slices of data.
  
• Wide boxplots mean the model is volatile, highly sensitive to training data, and risks overfitting.

You can also print the confusion matrix for each model by calling confusionMatrix() on the output of your train() function:

# Logistic Regression
confusionMatrix(model_glm)

Cross-Validated (10 fold) Confusion Matrix 

(entries are percentual average cell counts across resamples)
 
          Reference
Prediction  Yes   No
       Yes 42.3  7.3
       No   7.7 42.7
                            
 Accuracy (average) : 0.8498

# LDA
confusionMatrix(model_lda)

Cross-Validated (10 fold) Confusion Matrix 

(entries are percentual average cell counts across resamples)
 
          Reference
Prediction  Yes   No
       Yes 39.5  5.9
       No  10.5 44.1
                           
 Accuracy (average) : 0.836

# Naive Bayes
confusionMatrix(model_nb)

Cross-Validated (10 fold) Confusion Matrix 

(entries are percentual average cell counts across resamples)
 
          Reference
Prediction  Yes   No
       Yes 42.7  7.5
       No   7.3 42.5
                            
 Accuracy (average) : 0.8518

# KNN (This automatically grabs the confusion matrix for the best tuned 'k')
confusionMatrix(model_knn)

Cross-Validated (10 fold) Confusion Matrix 

(entries are percentual average cell counts across resamples)
 
          Reference
Prediction  Yes   No
       Yes 44.7  3.0
       No   5.3 47.0
                           
 Accuracy (average) : 0.917

When you run confusionMatrix(model_glm) on a cross-validated object, caret will print out a cross-validated confusion matrix showing percentages or overall counts across all folds, rather than a single evaluation on a held-out test set.

You can also see a breakdown of the Accuracy and Kappa metrics (including the minimum, median, mean, and maximum scores achieved across the folds):

# Print a statistical summary of all models across the 10 folds
summary(results)

Call:
summary.resamples(object = results)

Models: Logistic, LDA, NaiveBayes, KNN 
Number of resamples: 10 

Accuracy 
                Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
Logistic   0.8000000 0.8250000 0.8515686 0.8496471 0.8627451 0.9215686    0
LDA        0.7647059 0.8116516 0.8415686 0.8361418 0.8600000 0.9000000    0
NaiveBayes 0.8000000 0.8407843 0.8613725 0.8516365 0.8647247 0.9000000    0
KNN        0.8431373 0.9000000 0.9207843 0.9170588 0.9362745 0.9803922    0

Kappa 
                Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
Logistic   0.6000000 0.6500000 0.7037979 0.6991451 0.7250673 0.8425926    0
LDA        0.5306748 0.6233371 0.6828351 0.6722813 0.7200000 0.8000000    0
NaiveBayes 0.6000000 0.6811747 0.7225866 0.7030067 0.7293702 0.8000000    0
KNN        0.6866359 0.8000000 0.8415385 0.8340798 0.8722958 0.9607994    0

When you run summary(results) in caret, it will output both Accuracy and Kappa statistics automatically based on the out-of-sample folds, sparing you from writing a single line of matrix evaluation algebra.

Key Contrasts between modelr and caret

This table summarize the architectural differences:

Feature	modelr + Base R	caret
Formula Syntax	Divergent (glm vs lda vs knn matrix inputs)	Unified (Always train(formula, method = "..."))
Cross-Validation	Manual looping via map2() or for loops	Declarative via trainControl()
Data Pre-processing	Manual scaling inside the loop to avoid leakage	Automated via preProcess = c("center", "scale")
Hyperparameter Tuning	Requires nested loops	Built-in via tuneGrid or random search
Model Comparison	Manual calculation and tibble aggregation	Instant via resamples() wrapper

Quick Code Visualizer

This pipeline emphasizes how much smoother caret handles the process

Workflow Comparison: modelr vs. caret

Interactive Plot Code (ggplotly)

caret stores the cross-validation fold metrics inside results$values. We can use tidyr::pivot_longer() to clean this up, group by the model and metric to find the average performance across all 10 folds, and then feed it into ggplot2 and plotly.

# 1. Extract and clean the cross-validation metrics
results_tidy <- results$values %>%
  # Resamples outputs columns like "Logistic~Accuracy" and "Logistic~Kappa"
  pivot_longer(
    cols = -Resample,
    names_to = c("Model", "Metric"),
    names_sep = "~"
  ) %>%
  # Calculate the mean across all 10 cross-validation folds
  group_by(Model, Metric) %>%
  summarise(Mean_Value = mean(value), .groups = "drop")

# 2. Build the ggplot (sorted dynamically using reorder)
gg_metrics <- ggplot(results_tidy, aes(
  x = reorder(Model, Mean_Value), 
  y = Mean_Value, 
  fill = Model,
  text = paste("Model:", Model, "<br>Mean", Metric, ":", round(Mean_Value, 4))
)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~Metric, scales = "free_x") + 
  coord_flip() +  # Keeps model names highly readable on the y-axis
  labs(
    title = "Out-of-Sample Model Comparison (10-Fold CV)",
    x = "Model (Sorted from lowest to highest performance)",
    y = "Cross-Validated Mean Value"
  ) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set2")

# 3. Convert to an interactive plotly object
# We specify tooltip = "text" to give students custom, clean hover labels
ggplotly(gg_metrics, tooltip = "text")

Metric Deep Dive: Accuracy vs. Kappa

When caret evaluates a classification task, it defaults to two primary metrics: Accuracy and Kappa. Look at both to avoid pitfalls. Basically, accuracy asks “What did I get right?” while Kappa asks “Did I actually learn something, or did I get lucky?”

Accuracy: The Intuitive Baseline

\[Accuracy = \frac{\text{True Positives} + \text{True Negatives}}{\text{Total Observations}}\]

• What it is: The simple percentage of correct predictions out of all total predictions made.
  
• When to use it: When your classes are perfectly balanced (50/50 split). Because we split the Boston dataset explicitly at the median of the crime rate, our target classes are exactly balanced (50% above median, 50% below). In this specific scenario, Accuracy is a safe and highly reliable metric.
  
• The Flaw: Accuracy completely breaks down on imbalanced datasets. If 95% of a dataset belongs to Class A, a broken model that blindly guesses “Class A” every single time will achieve a 95% accuracy. It looks brilliant on paper but is practically useless.

Kappa (Cohen’s Kappa): The Skeptic’s Metric

\[\kappa = \frac{p_o - p_e}{1 - p_e}\]

Where \(p_o\) is the observed accuracy, and \(p_e\) is the expected accuracy if the model were just randomly guessing based on class distributions.

• What it is: A measure of how much better your classification model is performing relative to a completely random classifier guessing based on chance alone. Kappa values range from -1 to 1:
  • 1.0: Perfect agreement between predictions and reality.
    
  • 0.0: The model is performing exactly as well as random guessing.
    
  • Negative: The model is actually doing worse than random guessing.
    
• When to use it: Essential when dealing with imbalanced classes (e.g., credit card fraud detection where 0.1% of transactions are fraudulent, or rare disease diagnosis).
  
• Why it matters here: Even though our Boston subset is split 50/50, displaying Kappa alongside Accuracy provides a “sanity check.” If a model has 85% accuracy and a high Kappa (e.g., 0.70), you know it has genuinely learned patterns rather than exploiting a distribution quirk.