3.2. Logistic regression¶

Although having “regression” in its name, logistic regression is a classification (rather than regression) algorithm. Instead of modelling the qualitative response \(y\) directly, logistic regression models the probability that \(y\) belongs to a particular class (category).

Watch the 8-minute video below for a visual explanation of logistic regression:

Video

Explaining Logistic Regression by StatQuest, embedded according to YouTube’s Terms of Service.

3.2.1. Import libraries and load data¶

Get ready by importing the APIs needed from respective libraries.

import pandas as pd
import numpy as np
from matplotlib.gridspec import GridSpec
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, roc_auc_score, RocCurveDisplay

from statsmodels.formula.api import logit

%matplotlib inline

Load the Default dataset and inspect the first three rows.

data_url = "https://github.com/pykale/transparentML/raw/main/data/Default.csv"
df = pd.read_csv(data_url)

# Note: factorize() returns two objects: a label array and an array with the unique values. We are only interested in the first object.
df["default2"] = df.default.factorize()[0]
df["student2"] = df.student.factorize()[0]
df.head(3)

	default	student	balance	income	student2
0	No	No	729.526495	44361.625074	0
1	No	Yes	817.180407	12106.134700	1
2	No	No	1073.549164	31767.138947	0

3.2.2. Logistic model¶

Logistic regression models the probability that the response \(y\) belongs to a particular class (category) using linear regression via a logistic function, hence the name logistic regression. For binary classification, let us use 1 to denote the positive class (e.g. “Yes”) and 0 to denote the negative class (e.g. “No”).

3.2.2.1. Logit function transformation¶

Let \(\pi\) be the probability of a positive outcome

\[ \pi = \mathbb{P}(y=1| x). \]

Then the probability of a negative outcome \( \mathbb{P}(y=0| x) = 1-\pi \).

The odds of a positive outcome is the ratio of the probability of a positive outcome to the probability of a negative outcome, i.e.

\[ \frac{\pi}{1-\pi}. \]

The log of the odds is called the logit and the logit function is defined as :

\[ \mathrm{logit}(\pi)=\log \frac{\pi}{1-\pi} \]

The logit function is a monotonic transformation of the probability \(\pi\). It maps the probability \(\pi\) ranging from 0 to 1 to the real line ranging from \(-\infty\) to \(\infty\). The logit function is also known as the log-odds function.

3.2.2.2. Logistic function¶

In simple linear regression, we have modelled the relationship between outcome \(y\) and feature \(x\) with a linear equation:

\[ y = \beta_0 + \beta_1 x. \]

This linear regression model is not able to produce a binary outcome for binary classification or a probability (in the range from 0 to 1 only). With the logit function above, we can transform the probability \(\pi\) ranging from 0 to 1 to the logit function \(\mathrm{logit}(\pi)\) ranging from \(-\infty\) to \(\infty\). Now we can use the linear regression model to model the logit function \(\mathrm{logit}(\pi)\) instead of the probability \(\pi\) or \(y\) directly. The logit function is defined as:

\[ \mathrm{logit}(\pi)=\log \frac{\pi}{1-\pi} = \beta_0 + \beta_1 x. \]

The logistic function, also known as the sigmoid function, is the inverse of the logit function. It maps the logit function \(\mathrm{logit}(\pi)\) ranging from \(-\infty\) to \(\infty\) to the probability \(\pi\) ranging from 0 to 1:

\[ \pi = \mathbb{P}(y = 1 | x) = \mathrm{logit}^{-1}(\beta_0 + \beta_1 x) = \mathrm{logistic}(\beta_0 + \beta_1 x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1)}}. \]

3.2.3. Estimating the coefficients to make predictions¶

The coefficients of a logistic regression model can be estimated by maximum likelihood estimation (MLE). The likelihood function is

\[ \mathcal{L}(\beta_0, \beta_1) = \prod_{i:y_i=1} \mathbb{P}(y_i = 1) \prod_{i:y_i=0} (1 - \mathbb{P}(y_i = 1)). \]

The mathematical details are beyond the scope of this course. If interested, please read Section 4.3 of the Pattern Recognition and Machine Learning book for more details of the optimization for logistic regression.

To make predictions, taking the classification problem on the Default dataset, the probability of default given balance predicted by a logistic regression model is

\[ \mathbb{P}(\text{default} = \text{Yes} \mid \text{balance}, \beta_0, \beta_1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 \times \text{balance})}}. \]

3.2.3.1. Example explanation of system transparency¶

Let us to see the curve of the logistic function on the Default dataset.

sns.regplot(data=df, x="balance", y="default2", logistic=True)
plt.ylabel("Probability of default")
plt.xlabel("Balance")
plt.show()

We can see an S-shaped curve with the vertical axis ranging from 0 to 1, indicating the probability of a positive outcome (default) with respect to the horizontal axis \(x\) (balance).

In practice, we can use the scikit-learn or statsmodels package to fit a logistic regression model and make predictions.

Let us fit a logistic regression model using the scikit-learn API for logistic regression first.

clf = LogisticRegression(solver="newton-cg")
X_train = df.balance.values.reshape(-1, 1)
X_test = np.arange(df.balance.min(), df.balance.max()).reshape(-1, 1)
y = df.default2
clf.fit(X_train, y)
print(clf)
print("classes: ", clf.classes_)
print("coefficients: ", clf.coef_)
print("intercept :", clf.intercept_)

LogisticRegression(solver='newton-cg')
classes:  [0 1]
coefficients:  [[0.00549892]]
intercept : [-10.65132978]

Let us plot the fitted logistic regression model and study the system transparency.

prob = clf.predict_proba(X_test)

plt.scatter(X_train, y, color="orange")
plt.plot(X_test, prob[:, 1], color="lightblue")

plt.hlines(
    1,
    xmin=plt.gca().xaxis.get_data_interval()[0],
    xmax=plt.gca().xaxis.get_data_interval()[1],
    linestyles="dashed",
    lw=1,
)
plt.hlines(
    0,
    xmin=plt.gca().xaxis.get_data_interval()[0],
    xmax=plt.gca().xaxis.get_data_interval()[1],
    linestyles="dashed",
    lw=1,
)

plt.hlines(
    0.5,
    xmin=plt.gca().xaxis.get_data_interval()[0],
    xmax=plt.gca().xaxis.get_data_interval()[1],
    linestyles="dashed",
    lw=1,
)

plt.vlines(
    -clf.intercept_ / clf.coef_[0],
    ymin=plt.gca().yaxis.get_data_interval()[0],
    ymax=plt.gca().yaxis.get_data_interval()[1],
    linestyles="dashed",
    lw=1,
)

plt.hlines(
    (clf.predict_proba(np.asarray(1500).reshape(-1, 1)))[0][1],
    xmin=plt.gca().xaxis.get_data_interval()[0],
    xmax=1500,
    linestyles="dashed",
    lw=1,
)

plt.vlines(
    1500,
    ymin=plt.gca().yaxis.get_data_interval()[0],
    ymax=(clf.predict_proba(np.asarray(1500).reshape(-1, 1)))[0][1],
    linestyles="dashed",
    lw=1,
)

plt.ylabel("Probability of default")
plt.xlabel("Balance")
plt.yticks([0, 0.25, 0.5, 0.75, 1.0])
plt.xlim(xmin=-100)
plt.show()

In the above example, we learnt a logistic regression model \(f(x)\) with two parameters, \(\beta_0\) and \(\beta_1\), from the data, where

\(\beta_0 = -10.65133019 \)
\(\beta_1 = 0.00549892 \)

Using these two estimated parameters, we can examine the system logic of the logistic regression model to reveal its system transparency.

System transparency

When balance=1500 , the predicted probability of default is

\[ f(1500) = \frac{1}{1 + e^{\beta_0 + \beta_1 \times 1500}} = \frac{1}{1 + e^{-10.65133019 + 0.00549892 \times 1500}} = 0.08294763. \]

When the probability of default 0.5 (the commonly used threshold for binary classification), the corresponding balance can be calculated using the log odds equation above:

(3.1)¶\[\begin{align} \begin{aligned} \ln \left(\frac{\pi}{1-\pi}\right) &= \beta_0 + \beta_1 \times x \\ \ln \left(\frac{0.5}{1-0.5}\right) &= \beta_0 + \beta_1 \times x \\ x & = \frac{- \beta_0}{\beta_1} \\ x & = \frac{- (-10.65133019)}{0.00549892} \\ x & = 1936.9858426745614. \end{aligned} \end{align}\]

Therefore, to produce result Yes, i.e. for the probability of default \(>\) 0.5, the balance should be greater than 1936.9858426745614. Likewise, to produce result No, i.e. the probability of default \(<\) 0.5, the balance should be less than 1936.9858426745614.

Let us fit a logistic regression model using the [statsmodels API for logistic regression] next.

est = logit("default2 ~ balance", df).fit()
est.summary2().tables[1]

Optimization terminated successfully.
         Current function value: 0.079823
         Iterations 10

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	-10.651331	0.361169	-29.491287	3.723665e-191	-11.359208	-9.943453
balance	0.005499	0.000220	24.952404	2.010855e-137	0.005067	0.005931

We can see that the estimated parameters are the same as those from the scikit-learn API. Note the two last columns, “[0.025” and “0.975]”. These are Confidence Intervals, the very concept we covered in the last chapter on Linear Regression! Now, let us fit a logistic regression model using student status as the input variable.

est = logit("default2 ~ student", df).fit()
est.summary2().tables[1]

Optimization terminated successfully.
         Current function value: 0.145434
         Iterations 7

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	-3.504128	0.070713	-49.554094	0.000000	-3.642723	-3.365532
student[T.Yes]	0.404887	0.115019	3.520177	0.000431	0.179454	0.630320

The coefficient associated with student is positive, indicating that the probability of default is higher for students than non-students.

3.2.4. Multiple logistic regression¶

The logistic regression model can be extended to multiple features. A generalised logistic regression model for an instance \(\mathbf{x} = [x_1, x_2, \dots, x_D]^\top\) is

\[ \ln \left( \frac{\pi}{1-\pi} \right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_D x_D. \]

Let us fit a multiple logistic regression model using the scikit-learn API first.

Let us plot the fitted multiple logistic regression model and study the system transparency.

import warnings

warnings.filterwarnings("ignore")

X_train = df.loc[:, ["balance", "income", "student2"]]
y = df.default2

clf = LogisticRegression(solver="newton-cg", penalty="none", max_iter=1000)
clf.fit(X_train, y)
print(clf)
print("classes: ", clf.classes_)
print("coefficients: ", clf.coef_)
print("intercept :", clf.intercept_)

LogisticRegression(max_iter=1000, penalty='none', solver='newton-cg')
classes:  [0 1]
coefficients:  [[ 5.73139387e-03  2.83615167e-06 -6.51143222e-01]]
intercept : [-10.85258456]

Let us fit a multiple logistic regression model using the statsmodels API next.

est = logit("default2 ~ balance + income + student", df).fit()
est.summary2().tables[1]

Optimization terminated successfully.
         Current function value: 0.078577
         Iterations 10

	Coef.	Std.Err.	z	P>\|z\|	[0.025	0.975]
Intercept	-10.869045	0.492273	-22.079320	4.995499e-108	-11.833882	-9.904209
student[T.Yes]	-0.646776	0.236257	-2.737595	6.189022e-03	-1.109831	-0.183721
balance	0.005737	0.000232	24.736506	4.331521e-135	0.005282	0.006191
income	0.000003	0.000008	0.369808	7.115254e-01	-0.000013	0.000019

Similar parameters are estimated as using the scikit-learn API.

3.2.5. Evaluation of logistic regression¶

The evaluation of logistic regression is similar to the evaluation of linear regression. The main difference is that the evaluation metrics are different. For example, the mean squared error (MSE) is not a good metric for classification. Here we introduce two common metrics for classification: accuracy and area under receiver operating characteristic (ROC) curve, i.e. AUC (or more precisely, AUROC).

3.2.5.1. Accuracy¶

The accuracy is the fraction of correct predictions. It is defined as

\[ \text{accuracy} = \frac{\text{number of correct predictions}}{\text{total number of predictions}}. \]

Let us compute the training accuracy of the multiple logistic regression model trained using scikit-learn API above.

y_pred = clf.predict(X_train)
accuarcy = accuracy_score(y, y_pred)
print("Accuracy: ", accuarcy)

Accuracy:  0.9732

3.2.5.2. Area under the ROC curve (AUC)¶

The receiver operating characteristic (ROC) curve is a plot of the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings. The TPR and FPR are defined as

\[\begin{split} \begin{aligned} \text{TPR} &= \frac{\text{number of true positives}}{\text{number of positives}} \\ \text{FPR} &= \frac{\text{number of false positives}}{\text{number of negatives}} \end{aligned} \end{split}\]

The area under the ROC curve (AUC) is a popular measure of classifier performance. The AUC is defined as the area under the ROC curve. It is a number between 0 and 1, where 0.5 is the random guess and 1 is the perfect prediction. The AUC is a good metric for evaluating the classification since it considers all possible classification thresholds.

Let us plot the ROC curve of the multiple logistic regression model trained using scikit-learn API above.

y_prob = clf.predict_proba(X_train)[:, 1]

roc_auc = roc_auc_score(y, y_prob)
print("AUC: ", roc_auc)

clf_disp = RocCurveDisplay.from_estimator(clf, X_train, y)
plt.show()

AUC:  0.949552531739353

We can see clearly the trade-off between TPR and FPR. The AUC is 0.95, indicating that the model is good.

3.2.6. Exercise¶

1. Logistic regression is a classification algorithm.

    a. True
    
    b. False

Compare your answer with the solution below

a. True

2. All the following exercises will use the Weekly dataset.

Load the dataset, convert the categories to numerical values, and fit a logistic regression model using Volume variable as a predictor and Direction variable as the response using scikit-learn or statsmodel. Find out the fitted model’s weight and bias.

# Write your code below to answer the question