Review
- Gamma regression from last week: positive, right-skewed response; log link; constant dispersion; exp(b) - 1 as a percentage effect
- The three unifying features of generalized linear models (slides): a linear predictor, a random component from the exponential family, and a link function
- Deviance, residual deviance, and deviance residuals carry over unchanged; only the family and the link change this week
- Discuss the gamma regression multiple-choice questions
- Box-Cox Transformation and smearing factors
Required Packages and Data: ROCR, default.csv (10,000 credit-card customers: default, student, balance, income)
From Numbers to Yes/No Outcomes (15 min)
- Which model works for predicting which business variables (slides): default, churn, fraud, response to a campaign, disease
- Why not OLS? The linear probability model predicts impossible probabilities and violates normality and constant variance
credit = read.csv("default.csv", stringsAsFactors = TRUE)
credit$y = ifelse(credit$default == "Yes", 1, 0)
lpm = lm(y ~ balance, data = credit)
sum(fitted(lpm) < 0) # 3123 negative "probabilities"
max(fitted(lpm)) # only 0.27, even for the highest balances
- Bernoulli
- outcome: y = 1 with probability p, y = 0 with probability 1 - p
- Expected value E(y ) = 1 * p + 0 * (1 - p) = p.
- which link function or its inverse function g() can turn an arbitrary number a + b1 x1 + ... + bk xk into a probability such that p = g(a + b1 x1 + ... + bk xk)? A natural choice for g() is logistic (same as Sigmoid or Softmax for categorical probability distributions).
- The link function is log-odds (logit)
- Odds = p / (1 - p); log-odds (logit) = log(p / (1 - p)), which ranges over the whole real line
p = c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)
data.frame(p, odds = p/(1-p), logodds = log(p/(1-p)))
curve(1/(1 + exp(-x)), -6, 6) # logistic (S) curve
curve(log(x/(1 - x)), 0.01, 0.99) # logit curve
Key concept: The logistic regression model is log(p/(1-p)) = a + b1 x1 + ... + bk xk, or equivalently p = 1/(1 + exp(-(a + b1 x1 + ...))). The model is linear in the log-odds, not in the probability.
- Assumptions of logistic regression: binary (Bernoulli) response, independent cases, log-odds linear in the predictors, and no severe collinearity. No normality and no constant-variance assumption
Concept of MLE for Binary Data (20 min)
- Likelihood of one observation: p if y = 1 and 1 - p if y = 0; log-likelihood: y log(p) + (1 - y) log(1 - p)
- Computing the log-likelihood of observations y = c(1, 0, 1, 1, 0) under different values of p; which p gives the largest value?
LogLik.ob = function(y, p) y * log(p) + (1 - y) * log(1 - p)
LogLik.obs = function(y, p) sum(LogLik.ob(y, p))
y5 = c(1, 0, 1, 1, 0)
sapply(c(0.4, 0.5, 0.6, 0.7), function(p) LogLik.obs(y5, p))
# -3.77 -3.47 -3.37 -3.48 -> maximum at p = 0.6 = 3/5
- In logistic regression, p depends on x through a and b; maximum likelihood estimation (MLE) chooses a and b to maximize the total log-likelihood
- Comparisons of OLS and MLE: OLS minimizes squared errors and has a formula; MLE maximizes the likelihood and needs an iterative algorithm (Fisher scoring)
- Doing MLE ourselves with optim() and checking it against glm() later (balance measured in $100s):
negLL = function(par, x, y) {
p = 1 / (1 + exp(-(par[1] + par[2] * x)))
-LogLik.obs(y, p)
}
optim(c(0, 0), negLL, x = credit$balance/100, y = credit$y, method = "BFGS")$par
# -10.651 0.550
Explore the Data (10 min)
t = table(credit$default) # No 9667, Yes 333
prop.table(t) # default rate 0.0333, only 3.3% default: an imbalanced response, which matters when we evaluate the classifier
tapply(credit$balance, credit$default, mean)
tapply(credit$income, credit$default, mean)
aggregate(cbind(balance, income) ~ default, credit, mean) #instead of tapply()
boxplot(balance ~ default, data = credit) # defaulters: higher balances
Build and Interpret a One-Predictor Model (20 min)
m1 = glm(y ~ balance, family = binomial(link = "logit"), data = credit)
summary(m1)
# Estimate Std. Error z value
# (Intercept) -10.65 0.3612 -29.49
# balance 0.005499 0.0002204 24.95
- Regression line: log(p/(1-p)) = -10.65 + 0.005499 balance; z-tests (not t-tests) because the binomial dispersion is fixed at 1
- fitted(m1) -- predictions on the training data set
- confint.default(m1) -- 95% confidence intervals of estimated parameters
Key concept: exp(b) is the odds ratio: each one-unit increase of x multiplies the odds of y = 1 by exp(b). It does not multiply the probability.
# +$100 of balance multiplies the odds by 1.733 (+73%)
exp(100 * coef(m1)["balance"])
exp(100 * confint.default(m1)["balance", ]) # 1.660 1.810
predict(m1, data.frame(balance = c(1000, 1500, 2000)), type = "response")
# 0.006 0.083 0.586
-coef(m1)[1] / coef(m1)[2] # which balance makes logodds = 0 or p = 0.5: balance $1,937
- Remember type = "response" for probabilities; the default returns the log-odds
- The same $100 matters little at a $1,000 balance and a lot at $2,000: the effect on the probability depends on where you are on the S-curve
plot(credit$balance, credit$y)
curve(predict(m1, data.frame(balance = x), type = "response"), col = "red", add = TRUE)
Measuring Fit: Deviance (20 min)
- Deviance = -2 x log-likelihood of the fitted model (for 0/1 data the saturated model has log-likelihood 0)
- Null deviance: deviance of the model with no predictors, where every p equals the overall default rate
deviance.obs = function(y, p) -2 * LogLik.obs(y, p)
D = deviance.obs(credit$y, fitted(m1)) # 1596.5 = m1$deviance
D0 = deviance.obs(credit$y, rep(mean(credit$y), nrow(credit)))
# 2920.6 = m1$null.deviance
1 - D / D0 # pseudo R-squared: 0.453
pchisq(D0 - D, df = 1, lower.tail = FALSE) # chi-squared test: p < 2e-16
AIC(m1) # 1600.5 = D + 2k with k = 2
- Deviance residuals: each case's signed contribution; their squares sum to the residual deviance
rd = sign(credit$y - fitted(m1)) * sqrt(-2 * LogLik.ob(credit$y, fitted(fit)))
sum(rd^2) # 1596.5
all.equal(rd, residuals(fit), check.attributes = FALSE) # TRUE
- Goodness of model: pseudo R-squared (larger is better), chi-squared test against the null model, AIC (smaller is better)
Multiple Logistic Regression and Confounding (20 min)
m2 = glm(y ~ balance + income + student, family = binomial, data = credit)
summary(m2) # income: p = 0.71
exp(cbind(OR = coef(m2), confint.default(fit2)))
m3 = glm(y ~ balance + student, family = binomial, data = credit)
anova(m3, m2, test = "Chisq") # p = 0.71: drop income
AIC(m1, m3, m2) # 1600.5 1577.7 1579.5
- A paradox: students default more overall but less at the same balance
# alone, students look riskier: odds ratio 1.50
exp(coef(glm(y ~ student, binomial, credit))[2])
# at the same balance, students are safer: odds ratio 0.49
exp(coef(fit3)["studentYes"])
aggregate(balance ~ student, credit, mean) # students: higher balances
Key concept: Confounding: balance is related to both student status and default, so the student effect reverses once balance is held fixed. Interpret each coefficient as holding the other predictors constant.
Evaluating the Classifier (30 min)
- Split the data, fit on the training set, and predict probabilities for the test set
set.seed(2026)
case_select = sample(1:nrow(credit), size = 0.8 * nrow(credit))
train = credit[case_select, ]; test = credit[-case_select, ]
fitT = glm(y ~ balance + student, family = binomial, data = train)
pred = predict(fitT, test, type = "response")
- Confusion matrix at the cutoff 0.5:
prediction = ifelse(pred >= 0.5, "Yes", "No")
cm = table(Actual = test$default, prediction)
prediction
Actual No Yes
No 1923 3
Yes 52 22
mean((pred >= 0.5) == test$y) # accuracy 0.972
1 - mean(test$y) # always predicting "No": 0.963
- Concepts:
- False Positive (FP) = Type I Error
- False Negative (FN) / Type II Error
- Accuracy measures the overall proportion of correct predictions among all total predictions.
- Use kappa instead of accuracy in this case: predicting every case as No, leading accuracy 97%, but kappa = 0.
- Precision answers: "When the model predicts a positive outcome, how often is it correct?" It focuses on minimizing False Positives.
- Sensitivity (= Recall = true positive rate) answers: "Of all actual positive cases, how many did the model manage to catch?" It focuses on minimizing False Negatives
- Accuracy looks excellent, but sensitivity is only 26/74 = 35%: the model misses most defaulters.
- Conflicting Goals between Precision and Sensitivity: classify all cases as negative = no false positive; classifying all cases as positive = no false negative.
- Balance measures: F1 = 2 precision x recall / (precision + sensitivity)
-
Specificity (= true negative rate) answers the question: "Of all the individuals who are actually negative, how many did the model correctly identify as negative?
-
ROC curve: trade-off curve between negative positive rate (1 - specificity) and true positive rate (recall or sensitivity)
- Accuracy against the cutoff, the ROC curve between false positive rate and true positive rate, and AUC with ROCR:
library(ROCR)
predObj = prediction(pred, test$y)
perf = performance(predObj, "acc"); plot(perf)
max_ind = which.max(slot(perf, "y.values")[[1]])
slot(perf, "x.values")[[1]][max_ind] # best-accuracy cutoff: 0.36
roc = performance(predObj, "tpr", "fpr")
plot(roc, colorize = TRUE); abline(a = 0, b = 1, lty = 3)
performance(predObj, "auc")@y.values[[1]] # AUC 0.949
Key concept: The ROC curve shows the trade-off between sensitivity and the false positive rate over all cutoffs; AUC is the probability that a random defaulter gets a higher score than a random non-defaulter.
- Choosing a cutoff that reflects business costs: if a missed default costs 10 times a false alarm, the theoretical cutoff is 1/(1 + 10) = 0.09
pc = performance(predObj, "cost", cost.fp = 1, cost.fn = 10)
pc@x.values[[1]][which.min(pc@y.values[[1]])] # about 0.097
table(Predicted = ifelse(pred >= 1/11, "Yes", "No"), Actual = test$default)
# catches 58 of 74 defaulters (78%) at the cost of 122 false alarms
- The cutoff is a business decision, not a statistical one
Wrap-Up (5 min)
- Workflow: explore, fit with glm(..., family = binomial), interpret odds ratios, check deviance and AIC, then evaluate on test data with a confusion matrix, ROC, AUC and a cost-based cutoff
- Preview of next week
Homework
- Reading: Lecture Note 14 -- Logistic Regression
- Writing:
- Multiple Choice Questions (on ecourse.org): due before the next week
- Hands-on: Exercise 9 (heart disease prediction) of Lecture Note 14: due before the next week
|
