Logistic regression

Examples

A tiny example

import numpy as np
import pandas as pd

rng = np.random.default_rng(24)

def logistic(x):
    return 1 / (1 + np.exp(-x))

def data_generator(x, rng=rng):
    x = np.asarray(x, dtype=float)
    y = np.where(logistic(2 * x - 1) > rng.random(len(x)), "A", "B")
    return pd.DataFrame({"x": x, "y": y})

data = data_generator([0., 2., 3.])
data

	x	y
0	0.0	B
1	2.0	A
2	3.0	A

from sklearn.linear_model import LogisticRegression

X = data[["x"]].values
m = LogisticRegression(penalty=None).fit(X, data["y"])
print("intercept:", round(float(m.intercept_[0]), 4))
print("slope:    ", round(float(m.coef_[0][0]), 4))

intercept: 8.158
slope:     -8.5716

The model gives a probability for each class.

grid = np.arange(-1, 2.5, 0.5).reshape(-1, 1)
pd.DataFrame({"x": grid[:, 0],
              "P(A|x)": m.predict_proba(grid)[:, 0].round(3),
              "predicted": m.predict(grid)})

	x	P(A|x)	predicted
0	-1.0	0.000	B
1	-0.5	0.000	B
2	0.0	0.000	B
3	0.5	0.020	B
4	1.0	0.602	A
5	1.5	0.991	A
6	2.0	1.000	A

The prediction is the class with the highest probability, which with two classes means the threshold is 0.5.

The log-likelihood

We can compute the log-likelihood by hand and compare it to log_loss.

from sklearn.metrics import log_loss

def loglik(theta):
    b0, b1 = theta
    return (np.log(logistic(-b0))
            + np.log(logistic(b0 + 2 * b1))
            + np.log(logistic(-b0 - 3 * b1)))

theta_hat = [float(m.intercept_[0]), float(m.coef_[0][0])]
by_hand = loglik(theta_hat)
by_sklearn = -log_loss(data["y"], m.predict_proba(X), normalize=False)

print(f"by hand:  {by_hand:.6f}")
print(f"sklearn:  {by_sklearn:.6f}")

1

The three terms are \(P(\text{B}|0)\), \(P(\text{A}|2)\) and \(P(\text{B}|3)\), which are the observed labels.

by hand:  -17.143619
sklearn:  -0.000412

log_loss with normalize=False is the negative log-likelihood. The two agree.

Move the two parameters and watch the log-likelihood. The dashed line is the generator that produced the data. Maximum likelihood picks the parameters that make the observed labels as probable as possible.

With three data points the fit is far from the generator. This is not a bug. It is what fitting three points looks like.

With two predictors the probability is a surface and the decision boundary is a line. Move the threshold and watch the boundary shift without the probability changing at all.

Spam classification

The input is text, so we need a representation. Bag of words counts how often each word occurs in each email.

import numpy as np
import pandas as pd
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.linear_model import LogisticRegression

spam = pd.read_csv("https://go.epfl.ch/bio322-spam.csv", nrows=6000).dropna()
train, test = spam.iloc[:2000], spam.iloc[2000:4000]

vectorizer = CountVectorizer()
counts = vectorizer.fit_transform(train["text"].values).toarray()
freq = counts / np.maximum(counts.sum(axis=1, keepdims=True), 1)

print("emails:", counts.shape[0], " vocabulary:", counts.shape[1])

1

Bag of words. Each column counts how often one word occurs in one email.

2

Divided by the length of the email, so that a long email does not dominate.

emails: 2000  vocabulary: 28128

Each word is one predictor, so \(p\) is in the tens of thousands while \(n\) is 2000. With \(p > n\) the training data can be fitted perfectly, and it will be. This is why week 4 is about regularization.

m = LogisticRegression(penalty=None, max_iter=400)
m.fit(freq, train["label"].values)

test_counts = vectorizer.transform(test["text"].values).toarray()
test_freq = test_counts / np.maximum(test_counts.sum(axis=1, keepdims=True), 1)

print(f"training accuracy: {m.score(freq, train['label'].values):.3f}")
print(f"test accuracy:     {m.score(test_freq, test['label'].values):.3f}")

1

transform, not fit_transform. The vocabulary has to be the one the model was trained on, or the columns would mean different things.

training accuracy: 0.996
test accuracy:     0.968

The training accuracy is far above the test accuracy. That gap is the subject of the next two weeks.

words = np.array(vectorizer.get_feature_names_out())
order = np.argsort(m.coef_[0])
pd.DataFrame({"most ham-like": words[order[:10]],
              "most spam-like": words[order[-10:]][::-1]})

	most ham-like	most spam-like
0	schedul	your
1	enron	http
2	hourahead	more
3	louis	click
4	start	www
5	hour	here
6	messag	you
7	date	remov
8	if	onlin
9	origin	net