Logistic Regression and Classification: Method

Logistic Regression and Classification

The model

A binary response y takes the value 1 (the event) or 0. Logistic regression models the probability of the event as a function of the predictors through the log odds: ln[p / (1 − p)] = η = β₀ + β₁x₁ + … + β_kx_k, so that p = e^η / (1 + e^η) always lies between 0 and 1. Each coefficient is the change of the log odds per unit of its predictor, and e^β is the odds ratio: the factor by which the odds p / (1 − p) are multiplied. A categorical predictor with m levels enters as m − 1 dummy variables (0/1 indicators of every level but the first, the reference), so a dummy's coefficient compares its level with the reference and its exponential is the odds ratio of the two levels. The probit model replaces the logistic function by the standard normal distribution function, p = Φ(η); its coefficients are about 1.6 to 1.8 times smaller than the logit ones and have no odds-ratio reading, but the two fits are almost always equivalent in their predictions.

Maximum likelihood

Nothing is minimized in the least-squares sense; the coefficients maximize the likelihood of the observed 0/1 outcomes, L(β) = Π p_i^y_i (1 − p_i)^{1 − y_i}, or its logarithm ℓ. The tool uses Fisher scoring, the iteratively reweighted least squares of R's glm: starting from p_i = (y_i + ½)/2 it solves the weighted normal equations XᵀWXβ = XᵀWz with weights W = diag[p(1 − p)] (for the probit φ(η)² / (p(1 − p))) and working responses z = η + (y − p) / (dp/dη), recomputes p, and repeats until the deviance changes by less than 10⁻⁸ relative (R's rule; a few extra steps then polish the estimates to working precision). The deviance D = −2ℓ plays the role of the residual sum of squares: the residual deviance is that of the fitted model, the null deviance that of the model with the intercept alone, whose fitted probability is the overall event fraction. The deviance residual of an observation is ±√d_i, the square root of its contribution to D with the sign of y − p; R prints their five-number summary and so does the page.

Standard errors, tests and intervals

The covariance matrix of the estimates is the inverse Fisher information (XᵀWX)⁻¹ at the solution; the standard errors are the square roots of its diagonal. The Wald statistic z = b / SE is compared with the standard normal, and the Wald interval is b ± z_α/2 SE; exponentiating the limits gives the interval of the odds ratio. For small samples or large coefficients the Wald interval is known to be too narrow and symmetric where the likelihood is not: the profile-likelihood interval (R's confint for a glm) holds one coefficient fixed at a trial value, refits the others, and takes as limits the two values at which the deviance rises by the χ²₁ quantile, 3.84 at 95 %, above its minimum. The tool finds those two roots by bisection on the refitted deviance and prints them beside the Wald limits. The likelihood-ratio test of all predictors together is LR = D_null − D_res, χ² with as many degrees of freedom as there are predictor terms (the dummies counted separately); it is the analogue of the overall F test. AIC = D + 2p and BIC = D + p ln n penalize the fit by the number of coefficients p and compare models of the same data, lower being better. Pseudo R² measures restate the likelihood gain: McFadden 1 − ℓ/ℓ₀, Cox–Snell 1 − exp[2(ℓ₀ − ℓ)/n], which cannot reach 1, and Nagelkerke, the Cox–Snell value divided by its maximum. None is the fraction of variance explained; McFadden values between 0.2 and 0.4 already indicate a good fit to individual outcomes.

Separation and Firth's correction

When a predictor, or a combination of predictors, splits the events from the non-events completely (every event above a value, every non-event below), the likelihood keeps rising as the coefficient grows and the maximum-likelihood estimate does not exist; with quasi-complete separation the same happens along one direction with ties on the boundary. The iterations then run to the limit with fitted probabilities numerically 0 or 1 and enormous coefficients and standard errors; the tool, like R, detects the probabilities and warns. Firth's penalised likelihood ℓ + ½ ln |XᵀWX| (the Jeffreys prior) removes the first-order bias of the estimates and always has a finite maximum; the modified score equations Xᵀ(y − p + h(½ − p)) = 0, with h the leverages, are solved by Newton–Raphson, and the standard errors come from the same information matrix. The estimates are shrunk toward zero relative to the maximum likelihood ones; for a saturated one-predictor model they equal the log odds ratio of the 2 × 2 table with ½ added to every cell. Choose the Firth model when the separation warning appears, when the events are rare, or when a level of a categorical predictor has no events.

Goodness of fit: the Hosmer–Lemeshow test

The deviance of a model for individual 0/1 outcomes has no χ² distribution, so lack of fit is judged by grouping. The Hosmer–Lemeshow test sorts the observations by fitted probability into g groups (deciles of risk by default), compares the observed and expected numbers of events and non-events in each, and sums (O − E)²/E over the 2g cells; the statistic is χ² with g − 2 degrees of freedom. The groups are cut at the quantiles of the fitted values as R's hoslem.test does, so tied fitted values (a model with only categorical predictors) can merge groups, and the page says how many were formed. A small p-value says the probabilities are wrong somewhere in their range, which usually calls for a transformed predictor or an interaction, not for discarding the model; a large p-value is only the absence of evidence. The calibration plot draws the same groups: observed fraction against mean fitted probability, on the diagonal when the probabilities are calibrated.

Classification, ROC and cross-validation

To turn probabilities into decisions an observation is classified as an event when p is at least the cutoff (0.5 = the more probable class). The confusion matrix counts true and false positives and negatives; accuracy is the fraction correct, to be compared with the no-information rate of always predicting the majority; sensitivity (recall) is the fraction of events caught, specificity the fraction of non-events cleared, precision the fraction of predicted events that are events, F1 the harmonic mean of precision and sensitivity. All of them move with the cutoff. The ROC curve plots sensitivity against 1 − specificity for every possible cutoff, and its area (AUC, by the trapezoid rule over the distinct fitted values) is the probability that a random event receives a higher fitted probability than a random non-event: 0.5 is no discrimination, above 0.8 is good. Youden's index, sensitivity + specificity − 1, marks the cutoff that balances the two error rates with equal weight; the right cutoff in practice follows from the costs of the two errors and the prevalence. In-sample accuracy flatters a model that has been fitted to the same data; k-fold cross-validation shuffles the observations (with a seeded generator, so the folds are reproducible), deals them into k folds, fits the model to k − 1 of them and classifies the held-out fold, and reports the accuracy of those out-of-fold predictions together with the out-of-fold log-likelihood.

Predictions for new rows

For a new observation x the linear predictor is η = xᵀb with standard error √(xᵀ (XᵀWX)⁻¹ x), and the fitted probability is the inverse link of η. The confidence limits are the inverse link of η ± z SE: computed on the link scale and transformed, so the interval is asymmetric on the probability scale and never leaves 0 to 1 (R's predict with se.fit and the plogis of the limits, as on the UCLA page the default rows reproduce). Figure 3 draws the same construction along the first numeric predictor, the other numeric predictors at their means and the categorical ones at their reference levels.

Multinomial logistic regression

A response with K unordered categories is modelled against a baseline category (the first level, or the first named in Level Order): ln[P(Y = k) / P(Y = baseline)] = xᵀβ_k for each of the K − 1 other categories, so P(Y = k) = e^xᵀβ_k / Σ_j e^xᵀβ_j with β_baseline = 0. The (K − 1)p parameters are found by Newton–Raphson with the analytic Hessian, step-halved on the log-likelihood; the standard errors come from the inverse observed information. A coefficient is the change of the log relative risk of its category against the baseline per unit of the predictor, and its exponential the relative-risk ratio (R's multinom calls them relative risk ratios; they are odds ratios of the two categories). Classification is to the most probable category, and the confusion matrix is K × K.

Ordinal logistic regression (proportional odds)

For K ordered categories (listed from lowest to highest in Level Order) the cumulative model of McCullagh is ln[P(Y ≤ k) / P(Y > k)] = ζ_k − xᵀβ for k = 1 … K − 1, with K − 1 increasing thresholds ζ_k and one coefficient vector β shared by every threshold (the proportional-odds assumption). This is R's polr parametrization: a positive coefficient moves the response toward the higher categories, and e^β is the odds ratio of being above any threshold per unit of the predictor, the same for every threshold. The parameters are fitted by Newton–Raphson with the analytic Hessian and the thresholds kept ordered; the category probabilities are the differences of successive cumulative probabilities. The assumption can be checked by fitting the separate binary logistic regressions of Y ≥ k on the same predictors (the page's ordinal example does so) and comparing their slopes; when they differ much, the multinomial model is the alternative.

Linear discriminant analysis

Fisher's discriminant for two groups assumes that the predictors are normally distributed in each group with the same covariance matrix S and different means m₀, m₁. The discriminant coefficients a = S⁻¹(m₁ − m₀) define the score s = aᵀx − aᵀ(m₀ + m₁)/2 + ln(π₁/π₀), with priors π the group proportions, and the posterior probability of group 1 is e^s / (1 + e^s), a logistic function like the logit model's but with coefficients estimated from the means and the pooled covariance rather than by maximum likelihood. The Mahalanobis distance² (m₁ − m₀)ᵀ S⁻¹ (m₁ − m₀) measures how far apart the groups are in the metric of S. Discriminant analysis is more efficient than logistic regression when its normality assumption holds and worse when it does not (categorical or skewed predictors); the tool offers it as the classical alternative for the same data and reports the same classification measures, the ROC curve of the posterior and the cross-validated accuracy.

Reading the results page

The fitted model is written out with its numbers; the coefficient table gives the estimate, standard error, z, p-value and both kinds of interval; the odds-ratio table exponentiates them; the fit block gives the log-likelihoods, deviances, the likelihood-ratio test, AIC, BIC and the pseudo R²; the Hosmer–Lemeshow table and test follow; the classification block gives the confusion matrix, the five measures, the AUC with the ROC curve, the calibration plot and the probability curve; cross-validation and the predictions for new rows close the page. For the multinomial and ordinal models the coefficient tables are per comparison (or coefficients and thresholds), classification is by the most probable category, and the figure draws every category's probability against the first numeric predictor.

Limits

No interaction or polynomial terms (add the product or square as a column of the data), no offsets or weights, no grouped (binomial counts) data, no stepwise selection, no exact or conditional logistic regression, no mixed models; up to 5,000 observations, 40 coefficients, 20 levels per categorical predictor and 6 response categories; profile intervals up to 15 coefficients and not for the Firth fit; linear discriminant analysis for two groups only.

Verification

The library is checked by test-scripts/test_logistic_regression.php (776 checks) against the R output printed on the UCLA OARC "R Data Analysis Examples" pages for the open graduate-admissions data (glm with the logit and probit links: every coefficient, standard error, z and p-value, the deviances, AIC, log-likelihood, the Wald and profile-likelihood intervals, the odds ratios, the Wald tests of the rank terms, the predicted probabilities by rank and the GRE 200 row with its se.fit), the ordinal apply data (polr: coefficients, thresholds, standard errors, deviance, AIC, the two check regressions) and the hsbdemo program data (multinom: both coefficient vectors, standard errors, deviance, AIC, the fitted and predicted probabilities), and against an independent numpy / scipy re-implementation (20 seeded binary, 5 multinomial and 5 ordinal data sets at 10⁻⁸, Firth's estimates, the profile roots, the Hosmer–Lemeshow grouping, the AUC, the cross-validation folds, the discriminant coefficients) and closed forms (the Mann–Whitney form of the AUC, Firth's 2 × 2 identity, the intercept-only fit).