Logistic Regression

Regression models predict a continuous quantity, and a probability is one. A probability-predicting model becomes a classifier once we add a decision rule, such as "if the probability is 50% or more, call it a cat".

Logistic regression predicts probabilities, so strictly it is a regression method, even though machine learning usually files it under classification because that is what it is mostly used for. Methods like SVMs are classifiers in the stricter sense: they predict an outcome without a probability.

In scikit-learn:

X = np.array([[1, 1], [1, 2], [1, 3], [1, 4], [1, 5]])  # first column is a constant, for the intercept
y = np.array([0, 0, 0, 1, 1])
 
model = LogisticRegression(random_state=0).fit(X, y)
model.predict_proba(X)[:, 1]  # [0.063, 0.161, 0.353, 0.608, 0.816]
model.predict(X)              # [0, 0, 0, 1, 1]

Odds

The odds of an event are the probability it happens over the probability it doesn't:

From a linear model to a probability

The outcome follows a Bernoulli distribution (a binomial with one trial): 1 with some unknown probability $p$, 0 with probability $1 - p$. The goal is to estimate $p$ for any combination of the independent variables.

A linear combination of the variables,

can take any real value, so we need a link between the real line and the interval $(0, 1)$. The logit, the natural log of the odds, maps $(0, 1)$ onto the whole real line:

That has the probabilities on the x axis; we want them as the output. So set the logit equal to $z$ and invert:

This is the sigmoid or logistic function $\sigma(z)$. It maps any real input into the open interval $(0, 1)$, never quite reaching 0 or 1.

Fitting by maximum likelihood

The coefficients are usually found by maximum likelihood estimation. The likelihood $L(\boldsymbol{\beta})$ of a coefficient vector $\boldsymbol{\beta}$ is the product of the predicted probabilities of what actually happened:

where $\mathbf{x}_i$ is the feature vector of observation $i$, $y_i \in \{0, 1\}$ its outcome, and $p(\mathbf{x}_i) = \sigma(\mathbf{x}_i \cdot \boldsymbol{\beta})$ the predicted probability of the positive class.

A product of many small numbers underflows, so we maximise the log-likelihood instead:

There is no closed-form solution, so we climb it iteratively. The sigmoid's derivative, $\sigma'(z) = \sigma(z)(1 - \sigma(z))$, makes the gradient come out neatly:

That is, the error for each observation, weighted by its features and summed. The transpose is what makes the shapes line up: $X$ is $n \times m$ and the error is $n \times 1$, so $X^\top$ (which is $m \times n$) times the error gives one $m \times 1$ gradient entry per weight.

From scratch

1. Start with all weights at zero.
2. Compute scores $X\boldsymbol{\beta}$ (one weighted sum per row) and pass them through the sigmoid.
3. Subtract the predictions from the labels to get the error.
4. Take a step up the gradient (gradient ascent, since we are maximising), scaled by the learning rate.

def sigmoid(x):
    return 1 / (1 + np.exp(-x))
 
 
def predict_proba(x, weights):
    return sigmoid(np.dot(x, weights))
 
 
def logistic_regression(X, y, num_steps: int, learning_rate: float):
    weights = np.zeros(X.shape[1])
    for step in range(num_steps):
        error = y - predict_proba(X, weights)
        gradient = np.dot(X.T, error)
        weights += learning_rate * gradient
    return weights
 
 
weights = logistic_regression(X, y, num_steps=1_000_000, learning_rate=5e-5)
predict_proba(X, weights)  # [0.002, 0.024, 0.248, 0.817, 0.984]

Both models classify all five points correctly, but ours is much more confident. That is not a bug: scikit-learn's LogisticRegression applies L2 regularisation by default (C=1.0), which shrinks the coefficients and pulls the probabilities towards 0.5. With penalty=None, or a very large C, it matches the from-scratch fit. With data this cleanly separable the unregularised optimum is at infinity, so our weights depend on how long we run.

Why not linear regression?

A linear model treats the classes as the numbers 0 and 1 and fits the best straight line through them. That line is not a probability: it goes below 0 and above 1, and a few extreme $x$ values tilt it badly. Logistic regression stays bounded and behaves at both extremes.

Interpreting coefficients: the odds ratio

An odds ratio compares two odds. A loaded coin that lands heads 70% of the time has odds $0.7 / 0.3 \approx 2.33$; a fair coin has odds 1. So the odds ratio of heads, loaded against fair, is

In logistic regression, $e^{\beta_j}$ is the odds ratio for a one-unit increase in $x_j$, holding the other variables constant. Take a fictitious model of sleep apnoea against body weight, where weight has an odds ratio of 1.07:

• Each extra pound multiplies the odds of apnoea by 1.07, a 7% increase.
• Ten extra pounds multiply them by $1.07^{10} \approx 1.97$, nearly doubling the odds.
• This holds at any point on the weight scale. It is the odds that change by a constant factor, not the probability.

Original notebook

Tags: Statistics, Machine Learning