Logistic Regression and Generalized Linear Models

Logistic regression models a binary outcome’s conditional probability and can turn that probability into a class decision using a threshold. Viewing it as a generalized linear model (GLM) connects its likelihood, loss, and fitting procedure with regression models for other outcomes. We will follow that connection through binary labels, counts, and counts observed over different exposure times.

The three pieces of a GLM

For observation i, a GLM models the conditional mean \(\mu_i=E[Y_i\mid x_i]\). The linear predictor \(\eta_i=x_i^\top\beta\) combines the input columns and coefficients; an intercept is included through a column of ones. A link function connects the mean to that predictor: \(g(\mu_i)=\eta_i\), or \(\mu_i=g^{-1}(\eta_i)\). The response distribution describes variation around the mean and is drawn from an exponential family, which includes the examples in the table. The design can still contain transformed inputs and interactions; “linear” refers to the coefficients in the predictor.

For binary data, the logit is \(\log(p/(1-p))\) and its inverse is the sigmoid. A predictor value of 0 gives probability 0.5; \(\eta=\log2\) gives odds 2 and probability 2/3. A log link has inverse \(\mu=\exp(\eta)\), so \(\eta=\log3\) gives expected count 3. A predicted count mean can be noninteger even though each observed count is an integer. A canonical link ties the linear predictor to the distribution’s natural parameter in its exponential-family form. The table lists common choices: log is canonical for Poisson, but not for Gamma or the usual fixed-dispersion negative-binomial family. A finite log predictor gives a positive mean; numerical overflow and underflow still require care.

Where the loss comes from

For conditionally independent binary observations with modeled probabilities \(p_i\), the log-likelihood is \(\sum_i[y_i\log p_i+(1-y_i)\log(1-p_i)]\). Negating it gives the summed binary cross-entropy. Using the mean instead of the sum has the same unpenalized minimizer. For one observed positive label with predicted probability 0.8, the loss is \(-\log0.8\approx0.223\); predicting 0.2 gives about 1.609. This likelihood interpretation does not guarantee that a fitted probability model is correctly specified.

A Gaussian response with constant variance gives squared error up to a positive scale and a term independent of the mean parameters. Poisson deviance is twice the difference between the saturated-model log-likelihood (allowing a separate mean for each observation) and the fitted-model log-likelihood, so minimizing it gives the same unpenalized coefficient fit as maximizing the Poisson likelihood. Deviance and negative log-likelihood are not numerically identical. Losses can also be chosen for their prediction target or decision costs; Loss Functions and Surrogates discusses those interpretations.

How it is actually fitted

General logistic regression requires numerical optimization. One method is iteratively reweighted least squares (IRLS). For the Bernoulli logit model, set \(w_i=p_i(1-p_i)\) and \(z_i=\eta_i+(y_i-p_i)/w_i\), then fit z by weighted least squares. This gives a Newton step for the unpenalized likelihood when the weights are positive and the design has full rank. It is one fitting method; scikit-learn’s default solver uses L-BFGS. The small example below rejects tiny weights and a rank-deficient solve rather than hiding them with clipping. A production solver also needs safeguards for difficult steps.

In the code, Xd is the design matrix including the intercept, and b is the current coefficient estimate. At the initial predictor value 0, p=0.5 and w=0.25, so a positive observation has working response z=0+(1−0.5)/0.25=2. Each row gets weight w and working response z; multiplying the design and response by the square root of the weight turns weighted least squares into an ordinary least-squares solve. NumPy, SciPy, and scikit-learn are used here; the last exercise also uses statsmodels.

import numpy as np
from sklearn.linear_model import LogisticRegression
from scipy.special import expit

rng = np.random.default_rng(0)
n = 2000
X = rng.normal(size=(n, 3))
Xd = np.column_stack([np.ones(n), X])
beta_true = np.array([-0.5, 1.2, -0.8, 0.4])
y = (rng.random(n) < 1 / (1 + np.exp(-Xd @ beta_true))).astype(float)

b = np.zeros(4)
for it in range(30):
    eta = Xd @ b
    p = expit(eta)
    w = p * (1 - p)
    if np.min(w) < 1e-12:
        raise RuntimeError("Tiny IRLS weights: inspect separation or scaling")
    z = eta + (y - p) / w
    root_w = np.sqrt(w)
    b_new, _, rank, _ = np.linalg.lstsq(root_w[:, None] * Xd, root_w * z, rcond=None)
    if rank != Xd.shape[1]:
        raise RuntimeError("IRLS design is rank deficient")
    if np.abs(b_new - b).max() < 1e-10:
        b = b_new
        break
    b = b_new

else:
    raise RuntimeError("IRLS did not converge within 30 iterations")

assert np.max(np.abs(Xd.T @ (expit(Xd @ b) - y))) < 1e-7
sk = LogisticRegression(C=np.inf, tol=1e-10, max_iter=5000).fit(X, y)
print(f"converged in {it + 1} iterations")
print("IRLS   ", np.round(b, 6))
print("sklearn", np.round(np.r_[sk.intercept_, sk.coef_[0]], 6))
print("true   ", beta_true)
np.testing.assert_allclose(b, np.r_[sk.intercept_, sk.coef_[0]], atol=1e-6, rtol=0)
# converged in 6 iterations
# IRLS    [-0.484414  1.112975 -0.746429  0.458875]
# sklearn [-0.484414  1.112975 -0.746429  0.458875]
# true    [-0.5  1.2 -0.8  0.4]

The loop reports convergence only after its coefficient-change criterion is met. It also checks the likelihood gradient and compares coefficients with a tightly converged scikit-learn fit. Both coefficient vectors agree to the displayed six decimal places here and pass an absolute-tolerance check of 1e-6. This is a numerical check for this dataset. Newton’s local quadratic convergence needs a sufficiently close starting point and a smooth objective with a nonsingular Hessian; convexity alone does not guarantee rapid or global convergence. Separation and ill-conditioning can make this iteration fail.

The fitted coefficients differ from the generating coefficients because a finite sample has random outcomes. The solver comparison helps separate numerical disagreement from sampling variation here. Standard errors quantify sampling uncertainty under assumptions; a coefficient close to its generating value in one run does not by itself validate those uncertainty estimates.

Separation

Under complete separation, the unpenalized logistic likelihood has a supremum but no finite coefficient vector attaining it. Scaling a separating direction makes every fitted class probability approach its observed label. Not every coefficient must diverge along every fitting path. A numerical optimizer can stop at large finite values because its tolerance is satisfied, even though a finite unpenalized maximum does not exist.

import numpy as np
from sklearn.linear_model import LogisticRegression

X = np.array([[-2.0], [-1.0], [1.0], [2.0]])
y = np.array([0, 0, 1, 1])                      # perfectly separable

for C in (1e0, 1e2, 1e4, 1e8):
    m = LogisticRegression(C=C, tol=1e-10, max_iter=100_000).fit(X, y)
    print(f"C={C:8.0e}  coef {m.coef_[0, 0]:9.4f}"
          f"  P(y=1 | x=1) {m.predict_proba([[1.0]])[0, 1]:.8f}")
# C=   1e+00  coef    1.0066  P(y=1 | x=1) 0.73235312
# C=   1e+02  coef    3.9453  P(y=1 | x=1) 0.98102172
# C=   1e+04  coef    7.8441  P(y=1 | x=1) 0.99960810
# C=   1e+08  coef   16.3200  P(y=1 | x=1) 0.99999992

Each finite C in this table still specifies a penalized fit. In this two-class, full-rank example, positive L2 regularization gives a finite optimum; increasing C weakens that penalty and allows a larger coefficient. The limiting unpenalized problem has no finite solution. The code uses a tighter tolerance than the default so a very weak penalty is less likely to be confused with early numerical stopping. Report the penalty and convergence settings along with a coefficient of this size.

Separation can occur in small samples, high-dimensional designs, or rare categorical levels, and it can also arise from a genuinely deterministic relationship. Investigate prediction-time availability if a feature perfectly predicts the outcome, but separation alone is not evidence of leakage. Large coefficients, convergence warnings, and very extreme probabilities are clues; probabilities need not round to exactly 0 or 1 for separation to be present.

Scikit-learn applies regularization by default, supporting finite and more stable fits in many settings. The unpenalized comparisons here use C=np.inf with the current API. Removing the penalty does not ensure that a finite maximum exists or that the optimizer found it. Bias-reducing approaches such as Firth logistic regression provide another option under appropriate design conditions; they do not resolve leakage or identify a causal effect.

A log-linear mean for counts

import numpy as np
from sklearn.linear_model import LinearRegression, PoissonRegressor

rng = np.random.default_rng(2)
n = 4000
x = rng.normal(size=(n, 2))
counts = rng.poisson(np.exp(0.5 + 0.8 * x[:, 0] - 0.4 * x[:, 1]))

ols = LinearRegression().fit(x, counts)
poisson = PoissonRegressor(alpha=0).fit(x, counts)
print(f"OLS      coef {np.round(ols.coef_, 4)}  intercept {ols.intercept_:.4f}")
print(f"Poisson  coef {np.round(poisson.coef_, 4)}  intercept {poisson.intercept_:.4f}")
print(f"true          [0.8, -0.4]  intercept 0.5")
print(f"OLS predicts a negative count for {np.mean(ols.predict(x) < 0):.4f} of rows"
      f", min {ols.predict(x).min():.3f}")
print(f"Poisson minimum prediction {poisson.predict(x).min():.4f}")
# OLS      coef [ 1.9393 -1.0111]  intercept 2.4533
# Poisson  coef [ 0.7939 -0.4025]  intercept 0.5010
# true          [0.8, -0.4]  intercept 0.5
# OLS predicts a negative count for 0.1320 of rows, min -6.450
# Poisson minimum prediction 0.0449

The fitted Poisson coefficients are close to the generating log-mean coefficients in this sample. OLS fits an additive mean on the original scale, so its coefficients do not estimate those same parameters. The example compares representations of a log-linear mean; it does not imply that every count prediction task requires a full Poisson distributional assumption.

OLS produces negative predicted means for 13.2% of these training rows, down to −6.45. That is incompatible with a count mean. The log-link model produces positive means by construction in exact arithmetic. Positivity alone does not establish a correct mean model, calibrated uncertainty, or good held-out performance; this code reports fitted values on the data used for estimation.

For a Poisson response, conditional variance equals the conditional mean at fixed inputs. Overall sample variance can exceed the overall mean simply because the means vary across inputs, even in this correctly generated Poisson example. Extra conditional dispersion can invalidate Poisson-based standard errors. Under a correctly specified log mean and suitable regularity conditions, Poisson quasi-likelihood can still consistently estimate the mean parameters; this is not a finite-sample unbiasedness guarantee. Alternatives include appropriate robust covariance estimates or a justified negative-binomial variance model, often \(\mu+\alpha\mu^2\). Neither repairs a wrong mean automatically.

Exercises

1. Read a coefficient as an odds ratio. Fit a logistic regression on a binary feature and confirm that \(e^{\beta}\) equals the ratio of the odds in the two groups. Then calculate how the implied risk difference changes across hypothetical baseline probabilities.

Compare numerical agreement with the empirical odds ratio, then calculate probability changes at several hypothetical baselines.

import numpy as np
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(0)
n = 200_000
treated = (rng.random(n) < 0.5).astype(float)
p = np.where(treated == 1, 0.40, 0.20)
y = (rng.random(n) < p).astype(int)

m = LogisticRegression(C=np.inf, tol=1e-10, max_iter=5000).fit(treated.reshape(-1, 1), y)
beta = m.coef_[0, 0]

p1, p0 = y[treated == 1].mean(), y[treated == 0].mean()
empirical = (p1 / (1 - p1)) / (p0 / (1 - p0))
print(f"coefficient {beta:.6f}   exp(beta) {np.exp(beta):.6f}")
print(f"empirical odds ratio      {empirical:.6f}")

print("\nsame coefficient, different baselines:")
for base in (0.01, 0.10, 0.50, 0.90):
    odds = base / (1 - base) * np.exp(beta)
    print(f"  baseline {base:.2f} -> {odds / (1 + odds):.4f}"
          f"   risk difference {odds / (1 + odds) - base:+.4f}")
np.testing.assert_allclose(np.exp(beta), empirical, rtol=1e-6, atol=0)
# coefficient 0.981409   exp(beta) 2.668213
# empirical odds ratio      2.668213
# same coefficient, different baselines:
#   baseline 0.01 -> 0.0262   risk difference +0.0162
#   baseline 0.10 -> 0.2287   risk difference +0.1287
#   baseline 0.50 -> 0.7274   risk difference +0.2274
#   baseline 0.90 -> 0.9600   risk difference +0.0600

2. Logit against probit. Fit both links to the same data and compare their coefficients and their predicted probabilities. Explain what the ratio between the coefficients is.

Inspect every coefficient ratio and both the mean and maximum probability differences. Is the rescaling exact?

import numpy as np
from scipy.optimize import minimize
from scipy.special import expit, ndtr, log_ndtr

rng = np.random.default_rng(1)
n = 20_000
X = np.column_stack([np.ones(n), rng.normal(size=n), rng.normal(size=n)])
true = np.array([-0.3, 1.0, -0.6])
y = (rng.random(n) < expit(X @ true)).astype(float)
sign = 2 * y - 1

def logit_loss_grad(b):
    eta = X @ b
    return np.mean(np.logaddexp(0, eta) - y * eta), X.T @ (expit(eta) - y) / n

def probit_loss_grad(b):
    eta = X @ b
    logp = log_ndtr(sign * eta)
    logpdf = -0.5 * eta ** 2 - 0.5 * np.log(2 * np.pi)
    grad = -(X.T @ (sign * np.exp(logpdf - logp))) / n
    return -np.mean(logp), grad

fits = []
for objective in (logit_loss_grad, probit_loss_grad):
    result = minimize(objective, np.zeros(3), jac=True, method="BFGS", options={"gtol": 1e-8})
    if not result.success:
        raise RuntimeError(result.message)
    fits.append(result.x)
logit, probit = fits
print("logit ", np.round(logit, 4))
print("probit", np.round(probit, 4))
print("ratio ", np.round(logit / probit, 4), "   pi/sqrt(3) =", round(np.pi / np.sqrt(3), 4))
pl, pp = expit(X @ logit), ndtr(X @ probit)
print(f"max probability difference {np.abs(pl - pp).max():.4f}   mean {np.abs(pl - pp).mean():.5f}")
# logit  [-0.3124  0.9989 -0.5965]
# probit [-0.1865  0.5978 -0.3583]
# ratio  [1.6752 1.671  1.6649]    pi/sqrt(3) = 1.8138
# max probability difference 0.0117   mean 0.00410

3. Exposure and an offset. Model counts observed over unequal time windows, once with the window as a predictor and once as an offset. Compare estimated rate parameters and identify which fitting objectives are equivalent.

Compare an hours predictor, an unweighted rate fit, and an exposure-weighted rate fit. Check the last against an explicit offset model.

import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import PoissonRegressor

rng = np.random.default_rng(2)
n = 20_000
x = rng.normal(size=n)
hours = rng.uniform(1, 40, n)
rate = np.exp(-1.0 + 0.7 * x)
counts = rng.poisson(rate * hours)

as_feature = PoissonRegressor(alpha=0, tol=1e-10, max_iter=1000).fit(
    np.column_stack([x, hours]), counts)
print(f"raw hours feature   x coef {as_feature.coef_[0]:.4f}"
      f"  hours coef {as_feature.coef_[1]:.4f}  intercept {as_feature.intercept_:.4f}")
for name, weights in (("unweighted rate", None), ("weighted rate", hours)):
    fitted = PoissonRegressor(alpha=0, tol=1e-10, max_iter=1000).fit(
        x.reshape(-1, 1), counts / hours, sample_weight=weights)
    print(f"{name:18s} x coef {fitted.coef_[0]:.4f}  intercept {fitted.intercept_:.4f}")

offset_fit = sm.GLM(counts, sm.add_constant(x), family=sm.families.Poisson(),
                    offset=np.log(hours)).fit()
print(f"explicit offset    x coef {offset_fit.params[1]:.4f}  intercept {offset_fit.params[0]:.4f}")
print("true               x coef 0.7000  intercept -1.0000")
np.testing.assert_allclose(np.r_[fitted.intercept_, fitted.coef_], offset_fit.params, atol=1e-6, rtol=0)
# raw hours feature   x coef 0.6974  hours coef 0.0516  intercept 0.7978
# unweighted rate    x coef 0.7010  intercept -1.0034
# weighted rate      x coef 0.6979  intercept -1.0011
# explicit offset    x coef 0.6979  intercept -1.0011
# true               x coef 0.7000  intercept -1.0000