Logistic Regression (glm) — Fit and interpret logistic regression with glm(..., family=binomial) (intro) - R Programming

SECTION 1

Historical Context & Motivation

Much of classical statistics was designed around continuous outcomes — predicting a person's height, a stock's price, or the yield of a chemical reaction. Yet many of the most consequential questions in computer science and data science involve binary outcomes: Will the user click? Will the email be spam? Will the patient be readmitted? Ordinary least squares regression is a poor fit for such problems because it can produce predicted values outside the [0, 1] interval and violates the assumptions of normally distributed errors. Logistic regression was developed precisely to address this gap, providing a principled probabilistic framework for modeling dichotomous response variables.

1838 • The Logistic Function

Pierre-François Verhulst introduces the logistic curve to model population growth with an upper bound, establishing the S-shaped function that would later become the core of logistic regression.

1944 • Logit Model for Bioassay

Joseph Berkson coins the term logit and proposes using the log-odds transformation as a link function, arguing it was more natural than the probit model used in dose–response analysis.

1972 • Generalized Linear Models

John Nelder and Robert Wedderburn publish their landmark paper on Generalized Linear Models (GLMs), unifying logistic, Poisson, and normal regression under a single iteratively reweighted least squares framework.

1993 • R and glm()

Ross Ihaka and Robert Gentleman begin developing R at the University of Auckland. The glm() function becomes a first-class citizen in R's stats package, making logistic regression accessible with a single function call.

2010s • ML Baseline and Interpretability

As deep learning rises, logistic regression persists as the go-to interpretable baseline in machine learning pipelines, regulatory contexts (e.g., credit scoring), and A/B testing frameworks.

The central question logistic regression addresses is deceptively simple: given a set of predictor variables, how do we estimate the probability that a binary outcome equals 1 while respecting the constraint that probabilities must lie between 0 and 1? R's glm(..., family=binomial) syntax provides an elegant, one-line answer rooted in over 180 years of mathematical development.

SECTION 2

Core Principles & Definitions

Logistic regression belongs to the broader family of Generalized Linear Models (GLMs). A GLM extends ordinary linear regression by allowing the response variable to follow any distribution in the exponential family (binomial, Poisson, gamma, etc.) and by connecting the expected value of the response to a linear predictor through a link function. For logistic regression the response is binomial and the link function is the logit. Understanding four foundational ideas will anchor everything that follows.

Binary Response Variable

The outcome y takes only two values, coded 0 or 1. The model estimates P(y = 1 | X), the conditional probability of the 'success' class given the predictors.

The Logit Link Function

The logit transforms a probability p ∈ (0, 1) to the real line via log(p / (1 − p)). This maps the linear predictor Xβ to a valid probability through the inverse logistic function.

Maximum Likelihood Estimation

Unlike OLS, logistic regression uses iteratively reweighted least squares (IRLS) to maximize the log-likelihood. There is no closed-form solution for the coefficients β.

Odds and Odds Ratios

Exponentiated coefficients exp(β) represent multiplicative changes in the odds of the outcome for a one-unit increase in the predictor — the primary tool for interpretation.

✦

KEY TAKEAWAY

Think of logistic regression like a dimmer switch rather than an on/off toggle. The linear predictor Xβ can take any value on the real line (−∞ to +∞), and the logistic function smoothly 'squashes' it into the (0, 1) range — much like a dimmer maps a continuous dial position to a bounded light intensity. The further Xβ moves from zero, the closer the predicted probability gets to 0 or 1, but it never actually reaches either extreme.

SECTION 3

The Logistic (Sigmoid) Curve

The defining visual of logistic regression is the sigmoid curve, which maps any real-valued linear predictor to a probability between 0 and 1. The diagram below contrasts a naive linear fit (which can exceed [0, 1]) with the logistic curve, and annotates the critical inflection point at p = 0.5 where the model's decision boundary typically lies.

The logistic curve (gradient line) maps the linear predictor Xβ to valid probabilities in (0, 1), while the linear fit (dashed red) can produce impossible probability estimates above 1 or below 0. The yellow circle marks the inflection point at p = 0.5, where Xβ = 0.

Observe that the sigmoid curve is nearly flat in the tails and steepest at the center. This behavior means that changes in the predictor Xβ near zero have the largest impact on the predicted probability, while extreme predictor values shift probability very little. This is a fundamental property of the nonlinear relationship between predictors and probability in logistic regression, and it has direct implications for how we interpret coefficients: a one-unit change in a predictor does not correspond to a constant change in probability.

SECTION 4

Mathematical Framework

The mathematical core of logistic regression involves three related representations: the logit form (log-odds), the probability form (inverse logit), and the log-likelihood that is maximized to estimate coefficients. Understanding all three is essential for interpreting glm() output in R.

LOGIT (LOG-ODDS) FORM

logit(p) = log(p / (1 − p)) = β₀ + β₁x₁ + β₂x₂ + ⋯ + βₖxₖ

Where p = P(y = 1 | X), β₀ = intercept (log-odds when all predictors are 0), and βⱼ = change in log-odds for a one-unit increase in xⱼ, holding all other predictors constant.

PROBABILITY (INVERSE LOGIT) FORM

p = 1 / (1 + e^(−(β₀ + β₁x₁ + ⋯ + βₖxₖ)))

This is the sigmoid function applied to the linear predictor η = Xβ. It guarantees p ∈ (0, 1) for any real η.

ODDS RATIO INTERPRETATION

odds ratio for xⱼ = exp(βⱼ)

If exp(βⱼ) = 1.5, then a one-unit increase in xⱼ multiplies the odds of y = 1 by 1.5 (a 50% increase in odds). If exp(βⱼ) < 1, the odds decrease.

LOG-LIKELIHOOD

ℓ(β) = Σᵢ [ yᵢ log(pᵢ) + (1 − yᵢ) log(1 − pᵢ) ]

The log-likelihood sums over all n observations. R's glm() maximizes this function via IRLS (Iteratively Reweighted Least Squares). The deviance reported in the output equals −2ℓ(β) and serves as the GLM analogue of the residual sum of squares.

⚠

Why Not OLS?

You might wonder why we don't just regress y (0/1) on X using lm(). Three problems arise: (1) predicted values can fall outside [0, 1], (2) the error variance is heteroscedastic (it depends on p), and (3) the relationship between predictors and probability is fundamentally S-shaped, not linear. The logit link and binomial likelihood jointly resolve all three.

SECTION 5

The glm() Function — Syntax & Output Anatomy

In R, logistic regression is fit with a single call to glm() by specifying family = binomial (which defaults to the logit link). The function signature mirrors lm(), accepting a formula and a data frame, but internally it uses IRLS rather than QR decomposition. Below is a detailed anatomy of the function call and its output components.

The diagram breaks the glm() call into its four arguments (formula, data, family, output object) and shows the key components of the summary() output: the coefficients table, deviance information, AIC, and commonly used extractor and interpretation functions.

Output Component	Meaning	R Accessor
`Estimate`	Coefficient β̂ on the log-odds scale	`coef(model)`
`Std. Error`	Standard error of β̂, derived from the Fisher information matrix	`summary(model)$coefficients[,2]`
`z value`	Wald statistic: β̂ / SE(β̂), analogous to the t-statistic in OLS	`summary(model)$coefficients[,3]`
`Pr(>\|z\|)`	Two-sided p-value for the Wald test H₀: βⱼ = 0	`summary(model)$coefficients[,4]`
Null Deviance	−2 × log-likelihood of the intercept-only model	`model$null.deviance`
Residual Deviance	−2 × log-likelihood of the fitted model; drop = improvement	`model$deviance`
AIC	Akaike Information Criterion: −2ℓ + 2k; lower is better	`AIC(model)`

Key components of summary(model) for a logistic regression fit

SECTION 6

Worked Example — Predicting Spam Emails

Suppose we have a data frame email_df with 1000 emails, a binary column spam (1 = spam, 0 = not spam), and two numeric predictors: num_links (count of hyperlinks in the email body) and caps_pct (percentage of uppercase characters). We will fit a logistic regression model, examine the output, and make a prediction.

Fitting and Interpreting a Spam Classifier with glm()

1
Step 1 — Fit the Model

We call glm() with the formula spam ~ num_links + caps_pct and family = binomial. The R code is: spam_model <- glm(spam ~ num_links + caps_pct, data = email_df, family = binomial)

R returns a glm object containing coefficient estimates, deviances, and convergence info.
2
Step 2 — Examine the Summary

Calling summary(spam_model) yields (hypothetical output): (Intercept) -3.20 0.45 -7.11 <0.001 *** num_links 0.52 0.08 6.50 <0.001 *** caps_pct 0.10 0.03 3.33 0.0009 *** Null deviance: 1200.5 on 999 df. Residual deviance: 850.3 on 997 df. AIC: 856.3.

Both predictors are statistically significant (p < 0.001). The model reduces deviance by 350.2 from the null model.
3
Step 3 — Interpret Coefficients as Odds Ratios

Exponentiate with exp(coef(spam_model)): • exp(0.52) ≈ 1.68 — each additional link multiplies the odds of spam by 1.68 (68% increase). • exp(0.10) ≈ 1.11 — each 1 percentage-point increase in uppercase characters multiplies the odds by 1.11 (11% increase).

exp(β_num_links) ≈ 1.68, exp(β_caps_pct) ≈ 1.11
4
Step 4 — Predict for a New Email

Consider an email with 5 links and 15% uppercase characters. The linear predictor is: η = −3.20 + 0.52 × 5 + 0.10 × 15 = −3.20 + 2.60 + 1.50 = 0.90 Apply the inverse logit: p̂ = 1 / (1 + e^−0.90) ≈ 1 / (1 + 0.4066) ≈ 0.711. In R: predict(spam_model, newdata = data.frame(num_links=5, caps_pct=15), type='response')

Predicted P(spam) ≈ 0.711. At a 0.5 decision threshold, this email would be classified as spam.
5
Step 5 — Assess Model Fit

We compare null and residual deviance. The difference (1200.5 − 850.3 = 350.2) on 2 df is highly significant by a likelihood-ratio test: 1 - pchisq(350.2, df=2) ≈ 0, confirming that the predictors collectively improve the model. We can also compute McFadden's pseudo-R²: 1 − (850.3 / 1200.5) ≈ 0.292, indicating a moderate fit.

Deviance reduction of 350.2 on 2 df (p ≈ 0). McFadden's pseudo-R² ≈ 0.29.

SECTION 7

Strengths, Limitations & When to Use Logistic Regression

Logistic regression occupies a sweet spot in the modeling toolkit: it is powerful enough to reveal statistically significant relationships and produce calibrated probability estimates, yet simple enough to be fully interpretable. However, like any model, it comes with trade-offs that every practitioner must understand.

Strengths	Limitations
Produces well-calibrated probabilities, not just class labels — critical for risk scoring and decision-making under uncertainty.	Assumes a linear relationship between predictors and the log-odds; cannot automatically capture nonlinear patterns or interactions.
Coefficients have a direct odds-ratio interpretation, making results transparent and explainable to stakeholders.	Sensitive to multicollinearity — highly correlated predictors inflate standard errors and destabilize coefficient estimates.
Computationally efficient; IRLS typically converges in under 25 iterations even on moderately large data sets.	Can suffer from complete or quasi-complete separation when a predictor perfectly predicts the outcome, causing coefficients to diverge.
Well-studied statistical properties — asymptotic normality of MLEs, Wald tests, likelihood-ratio tests, and confidence intervals are all readily available.	Basic form handles only binary outcomes; ordinal or multinomial responses require extensions (ordinal logistic, multinomial logistic).
Acts as a strong baseline model in ML pipelines; outperforms complex models when the true decision boundary is approximately linear.	No built-in regularization in base R's glm(); use glmnet for L1/L2 penalties when p is large relative to n.

Strengths and limitations of logistic regression via glm()

🎯

WHEN TO REACH FOR LOGISTIC REGRESSION

Think of logistic regression as the printf debugging of classification: it is never glamorous, but it is the first tool you should try because it establishes a clear, interpretable baseline. If a neural network achieves 93% accuracy on your spam dataset but logistic regression gets 90%, you need a strong justification for the added complexity. In regulated domains like finance or healthcare, the interpretability of odds ratios is not optional — it is a requirement.

SECTION 8

Connection to Advanced Methods

Logistic regression via glm() is the starting point for a family of increasingly flexible models. Understanding how it connects to these extensions helps you know when to graduate to more powerful tools and what you trade away when you do.

Feature	glm() Logistic Regression	Advanced Extension
Regularization	No penalty; MLE only	glmnet: L1 (Lasso), L2 (Ridge), Elastic Net for variable selection and shrinkage
Nonlinear effects	Linear in log-odds; manually add polynomials or interactions	GAMs (mgcv::gam) fit smooth splines for each predictor automatically
Response categories	Binary (0/1) only	nnet::multinom for multinomial; MASS::polr for ordinal outcomes
Random effects	Fixed effects only; no grouping structure	lme4::glmer for mixed-effects logistic regression with hierarchical/clustered data
Prediction power	Strong baseline; linear decision boundary	Random forests, gradient boosting (xgboost), neural nets — nonlinear boundaries, higher capacity, less interpretability

From glm() logistic regression to advanced classifiers

A key insight is that logistic regression's loss function — the negative log-likelihood — is identical to the binary cross-entropy loss used in neural networks. When you train a neural network for binary classification with a single sigmoid output unit and no hidden layers, you are literally fitting logistic regression. This equivalence means that every concept you learn here — odds ratios, log-likelihood, deviance — transfers directly into deep learning diagnostics. Mastering glm() is therefore not just a statistics exercise; it is building foundational intuition for modern machine learning.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Why does logistic regression use the logit link function instead of modeling P(y = 1) directly as a linear combination of the predictors? Explain both the mathematical constraint being addressed and the interpretive benefit that the logit provides.

Show answer

A linear model P(y = 1) = β₀ + β₁x₁ + ⋯ can produce predicted probabilities outside the valid [0, 1] interval. The logit function, log(p / (1 − p)), maps (0, 1) to (−∞, +∞), allowing the right-hand side to be an unrestricted linear predictor. The inverse logit then maps back to valid probabilities. The interpretive benefit is that coefficients on the logit scale correspond directly to changes in log-odds; exponentiating them yields odds ratios, which have a natural multiplicative interpretation: 'each one-unit increase in x multiplies the odds by exp(β).'

Question 2: BASIC CALCULATION

A logistic regression model has intercept β₀ = −1.5 and one predictor with β₁ = 0.8. Calculate: (a) the log-odds when x₁ = 3, (b) the odds, (c) the predicted probability P(y = 1), and (d) the odds ratio associated with a one-unit increase in x₁.

Show answer

(a) log-odds = −1.5 + 0.8 × 3 = −1.5 + 2.4 = 0.9. (b) odds = exp(0.9) ≈ 2.460. (c) p = 1 / (1 + exp(−0.9)) = 1 / (1 + 0.4066) ≈ 0.711, or equivalently 2.460 / (1 + 2.460) ≈ 0.711. (d) The odds ratio is exp(0.8) ≈ 2.226, meaning that each one-unit increase in x₁ multiplies the odds of y = 1 by approximately 2.23 — more than doubling them.

Question 3: INTERMEDIATE

You fit two nested logistic regression models in R: m1 <- glm(admit ~ gre, data=grad, family=binomial) m2 <- glm(admit ~ gre + gpa + rank, data=grad, family=binomial) Model m1 has residual deviance 486.1 (df = 398) and m2 has residual deviance 458.5 (df = 394). Perform a likelihood-ratio test to determine whether the additional predictors (gpa, rank) significantly improve the model. Show the R code and interpret the result.

Show answer

The test statistic is the difference in deviances: G² = 486.1 − 458.5 = 27.6, on 398 − 394 = 4 degrees of freedom (3 for rank if it has 4 levels encoded as 3 dummies, plus 1 for gpa). In R: anova(m1, m2, test='Chisq'). Under H₀ (the restricted model is adequate), G² ~ χ²(4). Computing pchisq(27.6, df=4, lower.tail=FALSE) ≈ 1.5 × 10⁻⁵, which is far below any conventional significance level. We reject H₀ and conclude that gpa and rank together contribute significant explanatory power beyond gre alone. The deviance dropped by about 5.7%, corresponding to a meaningful improvement in log-likelihood.

Question 4: APPLIED

You are building a customer churn model for a SaaS product. After fitting churn_mod <- glm(churned ~ tenure_months + support_tickets + plan_tier, data=customers, family=binomial), you obtain: β_tenure_months = −0.08, β_support_tickets = 0.35, and for plan_tier (reference = 'free'): β_basic = −0.50, β_premium = −1.20. Interpret each coefficient in business terms. Then write the R code to predict churn probability for a premium-tier customer with 18 months of tenure and 3 support tickets.

Show answer

Interpretation: exp(−0.08) ≈ 0.923 — each additional month of tenure reduces the odds of churn by about 7.7%, consistent with loyalty effects. exp(0.35) ≈ 1.42 — each additional support ticket increases churn odds by 42%, suggesting unresolved issues drive attrition. exp(−0.50) ≈ 0.61 — basic-tier customers have 39% lower churn odds than free-tier users, and exp(−1.20) ≈ 0.30 — premium customers have 70% lower odds, indicating strong retention for paid plans. Prediction code: predict(churn_mod, newdata = data.frame(tenure_months = 18, support_tickets = 3, plan_tier = 'premium'), type = 'response'). This returns the estimated P(churned = 1). Using the coefficients: η = β₀ + (−0.08)(18) + (0.35)(3) + (−1.20) = β₀ − 1.44 + 1.05 − 1.20 = β₀ − 1.59. The exact probability depends on the intercept.

Question 5: CRITICAL THINKING

When fitting a logistic regression model, R warns: 'glm.fit: fitted probabilities numerically 0 or 1 occurred.' (a) What does this warning indicate about the data? (b) What consequences does it have for the estimated coefficients? (c) Propose two strategies for diagnosing and addressing the issue, with R code or packages where appropriate.

Show answer

(a) This warning signals complete or quasi-complete separation: at least one predictor (or linear combination of predictors) perfectly or nearly perfectly separates the two classes. For example, all observations with x > threshold might have y = 1. (b) The MLE does not exist in the finite sense — the algorithm pushes the relevant coefficient toward ±∞ to achieve a likelihood approaching the maximum. In practice, the estimate will be extremely large in magnitude with an inflated standard error, and Wald tests become unreliable. (c) Strategy 1: Penalized likelihood via Firth's bias-reduced logistic regression (package brglm2: glm(y ~ x, family=binomial, method='brglmFit')), which adds a Jeffreys prior penalty that guarantees finite estimates. Strategy 2: If the separating variable is categorical with sparse cells, consider collapsing categories or using exact logistic regression (package elrm or logistf). You can diagnose separation by inspecting a cross-tabulation (table(y, x)) for zero cells, or by using the detectSeparation package: glm(y ~ x, family=binomial, method='detect_separation').

SUMMARY

Lesson Summary

Summary

Logistic regression models binary outcomes by applying the logit link function — log(p / (1 − p)) — to map a linear predictor Xβ to probabilities in (0, 1) via the sigmoid function. Coefficients are estimated through maximum likelihood estimation (IRLS) rather than OLS, and the primary interpretation tool is the odds ratio: exp(βⱼ) gives the multiplicative change in odds for a one-unit increase in xⱼ.

In R, the entire workflow fits into glm(y ~ x, data, family = binomial). The summary() output reports the coefficients table (Wald z-tests), null and residual deviance for assessing overall fit, and AIC for model comparison. Predictions on the probability scale require predict(..., type = 'response'). This model serves as the interpretable, efficient baseline classifier from which all more complex binary classification methods — regularized regression, GAMs, tree ensembles, and neural networks — depart.