🗓️ Week 06
Logistic Regression

STAT 204 – Introduction to Statistical Data Analysis

Prof. Marcela Alfaro Cordoba

Statistics - UCSC

04 Nov 2025

Class Overview

Today’s Focus:

  • When and why we use logistic regression
  • The logistic function and odds ratios
  • Maximum likelihood estimation (theory)
  • Model interpretation
  • Cross-validation and model selection
  • Evaluation metrics for classification

Learning Goals: Understand the theory and practice of modeling binary outcomes

Motivation: Beyond Linear Regression

Linear regression works great for continuous outcomes:

  • Predicting house prices
  • Estimating test scores
  • Forecasting sales

But what about binary outcomes?

  • Disease: Yes/No
  • Pass/Fail
  • Click/Don’t Click
  • Survive/Die

Problem: Linear regression predicts values outside [0,1]!

Why Linear Regression Fails

Example: Predict exam pass (1) or fail (0) from hours studied

Problems:

  • Predictions not bounded to [0, 1]
  • Can’t interpret as probabilities
  • Violates assumptions (non-constant variance, non-normal errors)

The Logistic Function

We need a transformation that maps (-∞, ∞) to (0, 1)

The sigmoid function (logistic function) does exactly this:

\[p(x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}}\]

From Probabilities to Odds

Probability: \(p = P(Y=1)\)

Odds: \(\frac{p}{1-p} = \frac{P(Y=1)}{P(Y=0)}\)

Examples:

  • \(p = 0.5 \rightarrow\) odds = 1 (even odds, “50-50”)
  • \(p = 0.75 \rightarrow\) odds = 3 (3-to-1 odds in favor)
  • \(p = 0.9 \rightarrow\) odds = 9 (9-to-1 odds in favor)

Log-odds (logit): \(\log\left(\frac{p}{1-p}\right)\)

This is our linear predictor!

The Logistic Regression Model

The model links linear predictors to probability:

\[\log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\]

Or equivalently:

\[p = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}}\]

Key differences from linear regression:

  • Response follows Bernoulli distribution (not Normal)
  • We model log-odds, not the outcome directly
  • No closed-form solution—need iterative methods

Maximum Likelihood Estimation

For each observation: \(Y_i \sim \text{Bernoulli}(p_i)\)

Likelihood for observation i: \[L_i = p_i^{y_i}(1-p_i)^{1-y_i}\]

Why this works:

  • If \(y_i = 1\): \(L_i = p_i^1(1-p_i)^0 = p_i\)
  • If \(y_i = 0\): \(L_i = p_i^0(1-p_i)^1 = 1-p_i\)

Total likelihood (assuming independence): \[L(\beta) = \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}\]

Log-Likelihood Function

Log-likelihood (easier to work with): \[\ell(\beta) = \sum_{i=1}^n [y_i \log(p_i) + (1-y_i)\log(1-p_i)]\]

Substitute \(p_i = \frac{1}{1 + e^{-X_i\beta}}\):

\[\ell(\beta) = \sum_{i=1}^n \left[y_i(X_i\beta) - \log(1 + e^{X_i\beta})\right]\]

Goal: Find \(\hat{\beta}\) that maximizes \(\ell(\beta)\)

Problem: No closed-form solution! (Unlike OLS)

Why No Closed-Form Solution?

In OLS: Taking derivative and setting to zero gives: \[\hat{\beta}_{OLS} = (X'X)^{-1}X'y\]

In Logistic Regression: The score equation is: \[\frac{\partial \ell}{\partial \beta} = \sum_{i=1}^n X_i(y_i - p_i) = 0\]

But \(p_i\) depends on \(\beta\) nonlinearly: \[p_i = \frac{1}{1 + e^{-X_i\beta}}\]

This creates a nonlinear system with no closed-form solution.

Solution: Use iterative methods!

Iterative Methods for MLE

Newton-Raphson Algorithm:

  1. Start with initial guess \(\beta^{(0)}\)
  2. Update: \(\beta^{(t+1)} = \beta^{(t)} + [H(\beta^{(t)})]^{-1} \nabla \ell(\beta^{(t)})\)
    • \(\nabla \ell\): gradient (score function)
    • \(H\): Hessian matrix (second derivatives)
  3. Repeat until convergence

Alternatively: IRLS (Iteratively Reweighted Least Squares)

  • Each iteration solves a weighted least squares problem
  • Weights depend on current estimates of \(p_i\)

In R: glm() uses IRLS by default

Properties of MLEs

Under regularity conditions, MLEs are:

  1. Consistent: \(\hat{\beta} \xrightarrow{p} \beta\) as \(n \to \infty\)

  2. Asymptotically normal: \[\sqrt{n}(\hat{\beta} - \beta) \xrightarrow{d} N(0, I(\beta)^{-1})\]

where \(I(\beta)\) is the Fisher Information matrix

  1. Asymptotically efficient: Lowest variance among consistent estimators

  2. Invariant: If \(\hat{\beta}\) maximizes \(\ell\), then \(g(\hat{\beta})\) maximizes likelihood of \(g(\beta)\)

Practical implications:

  • Can construct Wald tests and confidence intervals
  • Standard errors from inverse Hessian
  • Large sample sizes give reliable inference

Example: Computing Probabilities

Given: \(\log\left(\frac{p}{1-p}\right) = -2 + 0.5 \times \text{StudyHours}\)

For a student who studied 6 hours:

# Step 1: Calculate log-odds
hours <- 6
log_odds <- -2 + 0.5 * hours
print(paste("Log-odds:", log_odds))
[1] "Log-odds: 1"
# Step 2: Convert to odds
odds <- exp(log_odds)
print(paste("Odds:", round(odds, 3)))
[1] "Odds: 2.718"
# Step 3: Convert to probability
prob <- 1 / (1 + exp(-log_odds))
# Or equivalently: prob <- odds / (1 + odds)
print(paste("Probability:", round(prob, 3)))
[1] "Probability: 0.731"

Interpretation: Student has 73.1% chance of passing.

Interpreting Coefficients

Remember: Coefficients represent change in log-odds, not probability!

For a one-unit increase in \(x_j\):

  • Log-odds change by \(\beta_j\)
  • Odds multiply by \(e^{\beta_j}\) (odds ratio)

Example: \(\beta_1 = 0.3\) for years of education

  • One additional year of education increases log-odds by 0.3
  • Odds of employment multiply by \(e^{0.3} = 1.35\)
  • Interpretation: “Each additional year of education increases odds of employment by 35%”

Unlike linear regression: Effect on probability depends on starting point!

Example: Odds Ratios

# Coefficient β = 0.3
beta <- 0.3

# Odds ratio
OR <- exp(beta)
print(paste("Odds Ratio:", round(OR, 3)))
[1] "Odds Ratio: 1.35"
# Interpretation
print(paste("Odds multiply by", round(OR, 2)))
[1] "Odds multiply by 1.35"
print(paste("Percent increase:", round((OR-1)*100, 1), "%"))
[1] "Percent increase: 35 %"

Confidence intervals for odds ratios: \[\exp(\beta \pm 1.96 \times SE(\beta))\]

🎯 Activity: Interpret Real Output

Pair Work (6 minutes):


Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -4.2100     0.8200  -5.134  < 0.001 ***
age              0.0450     0.0150   3.000  0.00270 ** 
income           0.0002     0.0001   2.000  0.04550 *  
education        0.3200     0.1200   2.667  0.00766 **

Tasks:

  1. Calculate the odds ratio for education (response variable is employed yes/no)
  2. Interpret the odds ratio in plain language
  3. Which predictor has the strongest effect on odds?
  4. For a 40-year-old with income=50000 and education=16, calculate predicted log-odds (then convert to probability)

Discuss your answers with your partner, then we’ll share them on Ed Discussion.

Model Evaluation: Why Not R²?

In linear regression: R² measures proportion of variance explained

Problem for binary outcomes:

  • Bernoulli variance depends on probability: \(p(1-p)\)
  • R² doesn’t make sense conceptually
  • Need different metrics!

Alternative approaches:

  1. Classification metrics (next)
  2. Likelihood-based measures (AIC, BIC, other)

The Confusion Matrix

Classify predictions using a threshold (usually 0.5)

Predicted: No Predicted: Yes
Actual: No True Negative (TN) False Positive (FP)
Actual: Yes False Negative (FN) True Positive (TP)

Key metrics:

  • Accuracy = \(\frac{TP + TN}{TP + TN + FP + FN}\)
  • Sensitivity (Recall, TPR) = \(\frac{TP}{TP + FN}\) — “How many actual positives did we catch?”
  • Specificity (TNR) = \(\frac{TN}{TN + FP}\) — “How many actual negatives did we correctly identify?”
  • Precision (PPV) = \(\frac{TP}{TP + FP}\) — “Of predicted positives, how many were correct?”

Model Selection: Cross-Validation

Question: Which predictors should we include in our logistic model?

Same challenge as linear regression:

  • Too few predictors: underfitting
  • Too many predictors: overfitting
  • Need to estimate out-of-sample performance

Solution: Cross-validation!

K-Fold Cross-Validation for Classification

Algorithm:

  1. Split data into K folds (typically K=5 or K=10)
  2. For each fold k:
    • Train model on K-1 folds
    • Predict on held-out fold k
    • Calculate performance metric (accuracy, AUC, etc.)
  3. Average metrics across all K folds

Metrics to track:

  • Accuracy: Overall correct classification rate
  • AUC: Area under ROC curve
  • Log-loss: \(-\frac{1}{n}\sum[y_i\log(p_i) + (1-y_i)\log(1-p_i)]\)
  • Brier score: \(\frac{1}{n}\sum(y_i - p_i)^2\)

Cross-Validation in R

library(caret)
library(titanic)
data(titanic_train)

# Remove missing values
titanic <- na.omit(titanic_train)

# Convert Survived to factor with proper names
# 0 = "Died", 1 = "Survived" (or any valid R names)
titanic$Survived_factor <- factor(titanic$Survived, 
                                  levels = c(0, 1),
                                  labels = c("Died", "Survived"))

# Set up 10-fold cross-validation
train_control <- trainControl(
  method = "cv",
  number = 10,
  summaryFunction = twoClassSummary,  # For ROC, Sens, Spec
  classProbs = TRUE,                   # Needed for ROC
  savePredictions = TRUE
)

# Fit model with CV (use Survived_factor, not Survived)
model_cv <- train(
  Survived_factor ~ Pclass + Sex + Age + Fare,
  data = titanic,
  method = "glm",
  family = "binomial",
  trControl = train_control,
  metric = "ROC"  # Optimize for AUC
)

# View results
print(model_cv)
model_cv$results

Comparing Models with Cross-Validation

# Model 1: Simple
model1_cv <- train(
  Survived_factor ~ Pclass + Sex + Age,
  data = titanic,
  method = "glm",
  family = "binomial",
  trControl = train_control,
  metric = "ROC"
)

# Model 2: More complex
model2_cv <- train(
  Survived_factor ~ Pclass + Sex + Age + Fare + SibSp + Parch,
  data = titanic,
  method = "glm",
  family = "binomial",
  trControl = train_control,
  metric = "ROC"
)

# Compare results
results <- resamples(list(Simple = model1_cv, Complex = model2_cv))
summary(results)
bwplot(results)  # Boxplot of performance across folds

Information Criteria: AIC and BIC

Alternative to CV: Use information criteria for quick model comparison

AIC (Akaike Information Criterion): \[\text{AIC} = -2\ell(\hat{\beta}) + 2p\]

BIC (Bayesian Information Criterion): \[\text{BIC} = -2\ell(\hat{\beta}) + p\log(n)\]

where:

  • \(\ell(\hat{\beta})\) is the log-likelihood at MLE
  • \(p\) is the number of parameters
  • \(n\) is the sample size

Lower is better! BIC penalizes complexity more than AIC.

Using AIC/BIC in R

# Fit several models
model1 <- glm(Survived_factor ~ Pclass + Sex, 
              data = titanic, family = binomial)
model2 <- glm(Survived_factor ~ Pclass + Sex + Age, 
              data = titanic, family = binomial)
model3 <- glm(Survived_factor ~ Pclass + Sex + Age + Fare + SibSp, 
              data = titanic, family = binomial)

# Compare AIC
AIC(model1, model2, model3)

# Compare BIC
BIC(model1, model2, model3)

# Stepwise selection using AIC
library(MASS)
full_model <- glm(Survived_factor ~ Pclass + Sex + Age + Fare + SibSp + Parch,
                  data = titanic, family = binomial)
step_model <- stepAIC(full_model, direction = "both", trace = FALSE)
summary(step_model)

🎯 Activity: Model Selection Practice

Pair Work (8 minutes):

library(caret)
library(titanic)
data(titanic_train)
titanic <- na.omit(titanic_train)

# IMPORTANT: Convert to factor with valid names
titanic$Survived_factor <- factor(titanic$Survived,
                                  levels = c(0, 1),
                                  labels = c("Died", "Survived"))

# Your task: Compare at least 3 different models

# Set up CV
train_control <- trainControl(
  method = "cv", number = 10,
  summaryFunction = twoClassSummary,
  classProbs = TRUE
)

# Model A: Your choice (use Survived_factor)
modelA <- train(Survived_factor ~ ..., data = titanic,
                method = "glm", family = "binomial",
                trControl = train_control, metric = "ROC")

# Model B: Your choice
modelB <- train(Survived_factor ~ ..., data = titanic,
                method = "glm", family = "binomial",
                trControl = train_control, metric = "ROC")

# Model C: Your choice
modelC <- train(Survived_factor ~ ..., data = titanic,
                method = "glm", family = "binomial",
                trControl = train_control, metric = "ROC")

# Compare
results <- resamples(list(A = modelA, B = modelB, C = modelC))
summary(results)

Questions: Which model has highest AUC? Does it match AIC/BIC?

ROC Curve and AUC

ROC (Receiver Operating Characteristic) Curve:

  • Plot Sensitivity (TPR) vs. 1-Specificity (FPR) at different thresholds
  • Shows trade-off between true positives and false positives

AUC (Area Under the Curve):

  • Single number summarizing model performance
  • AUC = 0.5: Random guessing
  • AUC = 1.0: Perfect classification
  • AUC > 0.8: Generally considered good

Use ROC/AUC when:

  • Don’t want to commit to a specific threshold
  • Want to compare models overall

Choosing a Threshold

Default: 0.5 — but not always optimal!

Considerations:

  • Medical screening: Lower threshold (higher sensitivity) — don’t want to miss diseases
  • Spam detection: Higher threshold (higher precision) — don’t want false positives
  • Legal applications: Balance based on costs of errors

Methods:

  • Maximize F1 score: \(\frac{2 \times \text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}\)
  • Minimize cost function based on domain knowledge
  • Youden’s index: max(Sensitivity + Specificity - 1)

Logistic Regression in R

Live Coding Demo: Titanic Survival

# Load data
library(titanic)
data(titanic_train)
titanic <- titanic_train

# Explore
head(titanic)
table(titanic$Survived)

# Fit logistic regression
model <- glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare,
             data = titanic,
             family = binomial(link = 'logit'))

# View results
summary(model)

# Odds ratios
exp(coef(model))
exp(confint(model))

# Predictions (probabilities)
titanic$pred_prob <- predict(model, type = 'response')

# Classify with threshold 0.5
titanic$pred_class <- ifelse(titanic$pred_prob > 0.5, 1, 0)

# Confusion matrix
table(Actual = titanic$Survived, Predicted = titanic$pred_class)

# Calculate metrics
library(caret)
confusionMatrix(factor(titanic$pred_class), 
                factor(titanic$Survived), 
                positive = "1")

# ROC curve
library(pROC)
roc_obj <- roc(titanic$Survived, titanic$pred_prob)
plot(roc_obj, main = paste("AUC =", round(auc(roc_obj), 3)))
auc(roc_obj)

Student Practice (10 minutes)

Your Turn:

library(titanic)
library(caret)
data(titanic_train)
titanic <- na.omit(titanic_train)

# For caret: Convert to factor with valid names
titanic$Survived_factor <- factor(titanic$Survived,
                                  levels = c(0, 1),
                                  labels = c("Died", "Survived"))

# 1. Fit a logistic regression model with your chosen predictors
your_model <- glm(Survived_factor ~ ...,  # glm() is fine with 0/1
                  data = titanic,
                  family = binomial)

# 2. Examine coefficients and calculate odds ratios
summary(your_model)
exp(coef(your_model))
exp(confint(your_model))

# 3. Compare with stepwise selection
library(MASS)
full_model <- glm(Survived_factor ~ Pclass + Sex + Age + Fare + SibSp + Parch, data = titanic, family = binomial)
step_model <- stepAIC(full_model, direction = "both")

# 4. Generate predictions
pred_probs <- predict(your_model, type = 'response')

# 5. Create confusion matrix (threshold = 0.5)
pred_class <- ifelse(pred_probs > 0.5, 1, 0)
table(Actual = titanic$Survived_factor, Predicted = pred_class)

# 6. Calculate metrics and ROC/AUC
library(caret)
confusionMatrix(factor(pred_class), factor(titanic$Survived), 
                positive = "1")

library(pROC)
roc_result <- roc(titanic$Survived_factor, pred_probs)
plot(roc_result)
auc(roc_result)

# 7. Use cross-validation to compare models
# IMPORTANT: Use Survived_factor for caret!
train_control <- trainControl(method = "cv", number = 10,
                              summaryFunction = twoClassSummary,
                              classProbs = TRUE)
cv_model <- train(Survived_factor ~ ..., data = titanic,
                  method = "glm", family = "binomial",
                  trControl = train_control, metric = "ROC")
print(cv_model)

Common Pitfalls to Avoid

❌ Don’t:

  • Interpret coefficients as effects on probability
  • Use R² from linear regression
  • Ignore class imbalance
  • Use only accuracy (misleading with imbalance!)
  • Forget to check for multicollinearity

✓ Do:

  • Interpret as odds ratios
  • Use classification metrics (sensitivity, specificity, AUC)
  • Consider stratified sampling or weights
  • Report multiple metrics
  • Check VIF for predictors

Dealing with Class Imbalance

Problem: 95% negative, 5% positive → model predicts all negative, gets 95% accuracy!

Solutions:

  1. Resampling:
    • Oversample minority class
    • Undersample majority class
    • SMOTE (Synthetic Minority Oversampling)
  2. Cost-sensitive learning:
    • Weight observations differently
    • Penalize errors on minority class more
  3. Adjust threshold:
    • Lower threshold to catch more positives
  4. Use better metrics:
    • F1 score, precision-recall curve, balanced accuracy

Comparing to Linear Regression

Feature Linear Regression Logistic Regression
Outcome Continuous Binary (0/1)
Distribution Normal Bernoulli
Link function Identity Logit
Interpretation Direct effect on Y Effect on log-odds
Estimation OLS (closed-form) MLE (iterative)
Predictions Any real number Probability [0,1]
Evaluation R², RMSE Accuracy, AUC, sensitivity

Summary

Key Takeaways:

  • Logistic regression is for binary outcomes
  • Models log-odds as linear function of predictors
  • Coefficients represent effects on log-odds, not probability
  • Odds ratios (\(e^\beta\)) are easier to interpret
  • Use classification metrics: accuracy, sensitivity, specificity, AUC
  • Confusion matrix shows how well we classify
  • ROC curve shows trade-offs across thresholds
  • Threshold choice depends on costs of different errors

Connection to what you know:

  • Similar to linear regression in structure
  • Different link function and distribution
  • Foundation for more complex models (multinomial, ordinal)

Final Questions & Next Steps

For next class:

  • Project Proposal Presentations

Office Hours: After class today