R - Logistic Regression

R - A Modern Approach to Logistic Regression with Tidymodels

Logistic Regression is a fundamental classification algorithm used when the response variable is categorical (e.g., True/False, Yes/No, 0/1). Instead of modeling the value of the response directly, it models the probability that the response belongs to a particular category. For a binary response, this probability is calculated using the logistic function (or sigmoid function).

The general mathematical equation for the probability (p) in logistic regression is:

p = 1 / (1 + e^-(β₀ + β₁x₁ + β₂x₂ + ...))

Where:

• p is the probability of the outcome of interest (e.g., P(y=1)).
• x₁, x₂, … are the predictor variables.
• β₀, β₁, … are the model coefficients, which are learned from the data.

While base R’s glm() function is the traditional tool for this, the tidymodels framework offers a modern, tidy, and comprehensive approach to modeling. It encourages best practices like splitting data for validation and provides a consistent interface for many different models.

The Tidymodels Workflow

We will use the built-in mtcars dataset to model whether a car has an automatic (am = 0) or manual (am = 1) transmission based on its engine specifications. Our modern workflow consists of these steps:

1. Setup: Load packages and prepare data.
2. Data Splitting: Divide data into training and testing sets.
3. Model Specification: Define the logistic regression model.
4. Model Fitting: Train the model on the training data.
5. Model Evaluation: Assess model performance on the test data.

Step 1: Setup and Data Preparation

First, we install and load the necessary packages. tidyverse is for data manipulation and tidymodels is for modeling. We will also convert our outcome variable am to a factor, which is the standard way to represent categorical variables in R for modeling.

# Install packages if you haven't already
# install.packages(c("tidyverse", "tidymodels"))

library(tidyverse)
library(tidymodels)

# Prepare the data
car_data <- mtcars %>%
  select(am, cyl, hp, wt) %>%
  mutate(am = factor(am, levels = c(0, 1), labels = c("automatic", "manual")))

# Glimpse the prepared data
glimpse(car_data)

This produces the following result, showing our variables and their types:

Rows: 32
Columns: 4
$ am  <fct> manual, manual, manual, automatic, automatic, automatic, automa…
$ cyl <dbl> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, …
$ hp  <dbl> 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 1…
$ wt  <dbl> 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3.150, …

Step 2: Split Data into Training and Testing Sets

A crucial step in modern machine learning is to split your data. We train the model on the training set and evaluate its performance on the unseen testing set. This gives us an honest estimate of how the model will perform on new data. We use initial_split from the rsample package for this.

# Set a seed for reproducibility
set.seed(123)

# Create the split object
data_split <- initial_split(car_data, prop = 0.75, strata = am)

# Create the training and testing data frames
train_data <- training(data_split)
test_data  <- testing(data_split)

Step 3: Specify and Fit the Model

With parsnip, we define our model independently of its implementation. We specify a logistic_reg() model and set its computational engine to "glm" (the same function used in the original tutorial).

# Specify the model
log_reg_spec <- logistic_reg() %>%
  set_engine("glm")

# Fit the model to the training data
log_reg_fit <- log_reg_spec %>%
  fit(am ~ cyl + hp + wt, data = train_data)

# View the model coefficients in a tidy format
tidy(log_reg_fit)

The tidy() function from the broom package provides a clean summary of our model:

# A tibble: 4 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)   13.9      6.99       1.99  0.0465
2 cyl            0.871     1.13       0.770 0.441
3 hp             0.0218    0.0210     1.04  0.298
4 wt            -7.39      3.40      -2.17  0.0298

Step 4: Evaluate Model Performance

Now, we use our fitted model (log_reg_fit) to make predictions on the test_data and see how well it performs. The yardstick package provides functions for this.

# Make predictions on the test set
predictions <- predict(log_reg_fit, new_data = test_data, type = "class") %>%
  bind_cols(test_data)

# Create a confusion matrix
conf_mat(predictions, truth = am, estimate = .pred_class)

# Calculate accuracy
accuracy(predictions, truth = am, estimate = .pred_class)

The evaluation output will look like this:

# Confusion Matrix
# Truth
# Prediction  automatic manual
#   automatic         4      1
#   manual            1      2

# Accuracy
# A tibble: 1 × 3
#   .metric  .estimator .estimate
#   <chr>    <chr>          <dbl>
# 1 accuracy binary         0.75

Conclusion

Our model correctly classified 6 out of 8 cars in the test set, for an accuracy of 75%. The confusion matrix shows it correctly identified 4 ‘automatic’ and 2 ‘manual’ cars, while misclassifying one of each. Rather than relying only on p-values from the training summary (which can be misleading), this evaluation on unseen data provides a much more reliable measure of the model’s predictive power. The tidymodels workflow enforces these modern best practices, leading to more robust and reliable models.

Further Learning

To improve your model, you could:

• Explore feature engineering with the recipes package to preprocess your variables.
• Use cross-validation with vfold_cv for an even more robust performance estimate.
• Learn more about interpreting model results with ROC curves (roc_curve from yardstick).

The official Tidymodels website is an excellent resource for diving deeper.