import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import matplotlib as mp
import sklearn
from IPython.display import Image, HTML
import statsmodels.api as sm
from sklearn import model_selection
from sklearn import metrics
%matplotlib inlineSo far we have seen linear regression: a continuous valued observation is estimated as a linear (or affine) function of the independent variables.
Today the observation is a 0/1 variable and we want a probability:
Regularization – a penalty on the size of \(\boldsymbol{\beta}\), for linear and logistic regression alike – follows in Lecture 11B.
Imagine that you are observing a binary variable – value 0 or 1.
That is, these could be pass/fail, admit/reject, Democrat/Republican, etc.
Assume there is some probability of observing a 1, and that probability is a function of certain independent variables.
So the key properties of a problem that make it appropriate for logistic regression are:
The following example was adapted from this URL which seems to be no longer available. There is an archive of the page and an archive of the dataset.
Let’s consider this question:
What is the probability I will be admitted to Grad School?
Let’s see how variables, such as,
affect admission into graduate school.
The response variable, admit/don’t admit, is a binary variable.
So there are three predictor variables: gre, gpa and rank.
Treating rank as a number assumes that each step down in prestige has the same effect on the log-odds. The alternative is one indicator column per rank (C(rank) in formula syntax), which costs two more coefficients. We keep the numeric version here; the categorical version is shown in Appendix: Regression Formulas with statsmodels.
The first rows of the data (400 applicants), and some summary statistics:
import os, urllib.request
def data_file(path):
"""Course data file: use a local copy (here, or in the repo root above a
notebook in jupyter_notebooks/), else fetch it from the Materials repo (Colab)."""
for prefix in ('', '../'):
if os.path.exists(prefix + path):
return prefix + path
os.makedirs(os.path.dirname(path) or '.', exist_ok=True)
urllib.request.urlretrieve(
'https://raw.githubusercontent.com/tools4ds/DS701-Materials-FA26/main/' + path, path)
return path
# original data source: http://www.ats.ucla.edu/stat/data/binary.csv
df = pd.read_csv(data_file('data/ats-admissions.csv'))
df.head(8)| admit | gre | gpa | rank | |
|---|---|---|---|---|
| 0 | 0 | 380 | 3.61 | 3 |
| 1 | 1 | 660 | 3.67 | 3 |
| 2 | 1 | 800 | 4.00 | 1 |
| 3 | 1 | 640 | 3.19 | 4 |
| 4 | 0 | 520 | 2.93 | 4 |
| 5 | 1 | 760 | 3.00 | 2 |
| 6 | 1 | 560 | 2.98 | 1 |
| 7 | 0 | 400 | 3.08 | 2 |
df.describe().round(3)| admit | gre | gpa | rank | |
|---|---|---|---|---|
| count | 400.000 | 400.000 | 400.000 | 400.000 |
| mean | 0.318 | 587.700 | 3.390 | 2.485 |
| std | 0.466 | 115.517 | 0.381 | 0.944 |
| min | 0.000 | 220.000 | 2.260 | 1.000 |
| 25% | 0.000 | 520.000 | 3.130 | 2.000 |
| 50% | 0.000 | 580.000 | 3.395 | 2.000 |
| 75% | 1.000 | 660.000 | 3.670 | 3.000 |
| max | 1.000 | 800.000 | 4.000 | 4.000 |
We can also plot histograms of the variables:
df.hist(figsize = (10, 4));
Note how
df.hist()automatically plots a histogram for each column in the dataframe as a subplot.
Mean admission rate against each feature, with standard-error bars:
For each value of a feature we plot the fraction of applicants admitted, and add error bars for the standard error of the mean. GPA takes almost as many values as there are applicants, so we aggregate it into 7 equal-width bins and plot each bin’s mean at the bin centre.
Note that there’s a greatly expanded
groupbysection in the Pandas refresher.
def admit_rate(groups):
"""Mean admission rate per group, with the standard error of the mean."""
grouped = groups['admit']
return grouped.mean(), grouped.std() / np.sqrt(grouped.count())
bins = np.linspace(df.gpa.min(), df.gpa.max(), 8)
gpa_bin = pd.cut(df.gpa, bins, include_lowest=True)
panels = [('Rank', admit_rate(df.groupby('rank'))),
('GRE', admit_rate(df.groupby('gre'))),
('GPA (bin centre)', admit_rate(df.groupby(gpa_bin, observed=True)))]
fig, axs = plt.subplots(1, 3, figsize=(13, 3.8), sharey=True)
for ax, (name, (means, errors)) in zip(axs, panels):
x = [v.mid for v in means.index] if name.startswith('GPA') else means.index
ax.errorbar(x, means, yerr=errors, marker='o', capsize=4)
ax.set_xlabel(name, fontsize=14)
ax.tick_params(labelsize=11)
axs[0].set_xticks([1, 2, 3, 4])
axs[0].set_ylabel('P[admit]', fontsize=14)
plt.tight_layout()
Finally, we plot admission status versus GRE score for each data point for each of the four ranks:
fig, axs = plt.subplots(2, 2, figsize=(10, 5), sharex=True, sharey=True)
for rank, ax in zip(range(1, 5), axs.flat):
df[df['rank'] == rank].plot.scatter('gre', 'admit', ax=ax)
ax.set_title(f'Rank {rank} Institutions')
plt.tight_layout()
What we want to do is to fit a model that predicts the probability of admission as a function of these independent variables.
Logistic regression is concerned with estimating a probability.
However, all that is available are categorical observations, which we will code as 0/1.
That is, these could be pass/fail, admit/reject, Democrat/Republican, etc.
Now, a linear function like \(\beta_0 + \beta_1 x\) cannot be used to predict probability directly, because
However, there is a transformation of probability that works: it is called log-odds.
For any probability \(p\), the odds is defined as \(p/(1-p)\), which is the ratio of the probability of an event to the probability of the non-event.
Notice that odds vary from 0 to \(\infty\), and odds < 1 indicates that \(p < 1/2\).
Odds are still not a good target for a linear function, for two reasons.
It stretches the interval \((0, 1)\) onto the whole real line, with \(p = 1/2\) at \(0\): exactly the range a line can cover.
pvec = np.linspace(0.01, 0.99, 100)
ax = plt.figure(figsize = (6, 4)).add_subplot()
ax.plot(pvec, np.log(pvec / (1-pvec)))
ax.tick_params(labelsize=12)
ax.set_xlabel('Probability', fontsize = 14)
ax.set_ylabel('Log-Odds', fontsize = 14)
ax.set_title(r'Logit Function: $\log (p/(1-p))$', fontsize = 16);
So the logistic regression model says that the log-odds of \(y = 1\) is a linear function of \(x\), a straight line \(\beta_0 + \beta_1 x\).
Solving for \(p(y{=}1 \mid x)\):
\[ \begin{aligned} \beta_0 + \beta_1 x &= \log \frac{p(y{=}1 \mid x)}{1-p(y{=}1 \mid x)} \\ e^{\beta_0 + \beta_1 x} &= \frac{p(y{=}1 \mid x)}{1-p(y{=}1 \mid x)} \quad \text{(exponentiate both sides)} \\ e^{\beta_0 + \beta_1 x} (1-p(y{=}1 \mid x)) &= p(y{=}1 \mid x) \quad \text{(clear the denominator)} \\ e^{\beta_0 + \beta_1 x} &= p(y{=}1 \mid x) + p(y{=}1 \mid x)e^{\beta_0 + \beta_1 x} \quad \text{(expand)} \\ \frac{e^{\beta_0 + \beta_1 x}}{1 +e^{\beta_0 + \beta_1 x}} &= p(y{=}1 \mid x) \quad \text{(collect terms)} \end{aligned} \]
The second line of the derivation says the odds are an exponential function of \(x\):
\[ \frac{p(y{=}1 \mid x)}{1-p(y{=}1 \mid x)} = e^{\beta_0 + \beta_1 x} = e^{\beta_0}\,\big(e^{\beta_1}\big)^{x}. \]
So each unit increase in \(x\) multiplies the odds by the same factor \(e^{\beta_1}\), whatever the starting point:
\[ \frac{\text{odds}(x+1)}{\text{odds}(x)} = e^{\beta_1}. \]
So, logistic regression fits a probability of the following form:
\[ p(y{=}1 \mid x) = \frac{e^{\beta_0+\beta_1 x}}{1+e^{\beta_0+\beta_1 x}}. \]
This is a sigmoid function; when \(\beta_1 > 0\),
We plot this on the next slide.
This type of sigmoid comes up so often that it gets its own symbol:
\[ \sigma(z) = \frac{e^{z}}{1+e^{z}} = \frac{1}{1+e^{-z}}, \qquad\text{so}\qquad p(y{=}1 \mid x) = \sigma(\beta_0 + \beta_1 x). \]
and is sometimes called the logistic function.
With \(d\) features we collect them, together with a leading \(1\) for the intercept, into \(\mathbf{x} = (1, x_1, \ldots, x_d)^T\) and the coefficients into \(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_d)^T\). The model is then
\[ p(y{=}1 \mid \mathbf{x}) = \sigma(\boldsymbol{\beta}^T\mathbf{x}), \qquad \text{logit}\, p(y{=}1 \mid \mathbf{x}) = \boldsymbol{\beta}^T\mathbf{x}. \]
Everything below is written for one feature; replace \(\beta_0 + \beta_1 x\) by \(\boldsymbol{\beta}^T\mathbf{x}\) and it holds for any number of features.
def sigmoid(z):
return 1 / (1 + np.exp(-z))
x = np.linspace(0, 50, 200)
fig, axs = plt.subplots(1, 2, figsize=(11, 3.5), sharey=True)
# left: beta_0 fixed, beta_1 varies
for b1 in [0.2, 0.4, 0.6, 1.0]:
axs[0].plot(x, sigmoid(-8 + b1 * x), label=rf'$\beta_1 = {b1}$')
axs[0].set_title(r'$\beta_0 = -8$ fixed: $\beta_1$ sets the steepness', fontsize=13)
# right: beta_1 fixed, beta_0 varies; mark where p = 1/2
for b0 in [-4, -8, -12, -16]:
line, = axs[1].plot(x, sigmoid(b0 + 0.4 * x), label=rf'$\beta_0 = {b0}$')
axs[1].plot(-b0 / 0.4, 0.5, 'o', color=line.get_color())
axs[1].axhline(0.5, color='gray', linestyle=':', linewidth=1)
axs[1].set_title(r'$\beta_1 = 0.4$ fixed: $\beta_0$ shifts the curve', fontsize=13)
for ax in axs:
ax.set_xlabel('$x$', fontsize=13)
ax.tick_params(labelsize=11)
ax.legend(fontsize=11)
axs[0].set_ylabel(r'$p(y{=}1 \mid x)$', fontsize=13)
plt.tight_layout()
So we had the logistic sigmoid function: \[ \sigma(z) = \frac{e^{z}}{1+e^{z}} = \frac{1}{1+e^{-z}}, \qquad\text{where}\qquad p(y{=}1 \mid x) = \sigma(\beta_0 + \beta_1 x). \]
which maps \((-\infty, \infty)\) to \((0, 1)\).
The inverse of that is confusingly called the logit function:
\[ \text{logit}(p) = \log \left(\frac{p}{1-p}\right) \]
that maps \(p \in (0, 1)\) to the real line \((-\infty, \infty)\).
So a compact way to describe logistic regression is that it finds regression coefficients \(\beta_0, \beta_1\) to fit:
\[ \text{logit}\left(p(y{=}1 \mid x)\right)=\log\left(\frac{p(y{=}1 \mid x)}{1-p(y{=}1 \mid x)} \right) = \beta_0 + \beta_1 x. \]
So, the best way to think of logistic regression is that we compute a linear function:
\[ \beta_0 + \beta_1 x, \]
and then map that to a probability using the \(\text{logistic}\) function:
\[ \sigma(\beta_0 + \beta_1 x) = \frac{e^{\beta_0+\beta_1 x}}{1+e^{\beta_0+\beta_1 x}}. \]
And again, in machine learning we just usually call this the sigmoid function, and denote it \(\sigma(z)\).
Let’s take a moment to compare linear and logistic regression.
In Linear regression we fit
\[ y_i = \beta_0 +\beta_1 x_i + \epsilon_i. \]
We do the fitting by minimizing the sum of squared errors \(\Vert\epsilon\Vert^2\). This can be done in closed form using either geometric arguments or by calculus.
Now, if \(\epsilon_i\) comes from a normal distribution with mean zero and some fixed variance, then minimizing the sum of squared errors is exactly the same as finding the maximum likelihood of the data with respect to the probability of the errors.
So, in the case of linear regression, it is a lucky fact that the MLE of \(\beta_0\) and \(\beta_1\) can be found by a closed-form calculation.
In Logistic regression we fit
\[ \text{logit}(p(y_i{=}1 \mid x_i)) = \beta_0 + \beta_1 x_i. \]
How should we choose parameters?
Here too, we use Maximum Likelihood Estimation of the parameters.
That is, we choose the parameter values that maximize the likelihood of the data given the model.
\[ \text{P}(y_i \mid x_i) = \left\{\begin{array}{lr}\sigma(\beta_0 + \beta_1 x_i)& \text{if } y_i = 1\\ 1 - \sigma(\beta_0 + \beta_1 x_i)& \text{if } y_i = 0\end{array}\right. \]
We can write this as a single expression. Writing \(p_i = p(y_i{=}1 \mid x_i) = \sigma(\beta_0 + \beta_1 x_i)\),
\[ \text{P}(y_i \mid x_i) = p_i^{\,y_i} (1-p_i)^{1-y_i}, \]
where we assume the parameters are fixed.
We can reinterpret this to express the likelihood of parameters \(\beta_0\), \(\beta_1\) given that we have observed the data \((x_i, y_i)\):
\[ L(\beta_0, \beta_1 \mid x_i, y_i) = p_i^{\,y_i} (1-p_i)^{1-y_i}. \]
That is one observation. We have \(n\) of them.
The observations are independent, so the likelihood of the data set is the product and its logarithm – the log-likelihood – is the sum:
\[ L(\boldsymbol{\beta}) = \prod_{i=1}^n p_i^{\,y_i}(1-p_i)^{1-y_i} \qquad\Longrightarrow\qquad \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \big[\, y_i \log p_i + (1-y_i)\log(1-p_i) \,\big]. \]
sklearn and every neural-network library report. Maximizing \(\ell\) and minimizing log-loss are the same thing.Just like before, we can create our design matrix of input observations with a leading column of 1s for the intercept, and the vector of output observations1:
\[ \mathbf{X} = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_n \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \]
And our initialized \(\boldsymbol{\beta}^{(0)} = [0, 0]^T\).
When computing the full vector of predicted probabilities \(\mathbf{p}\):
\[\mathbf{p} = \sigma(X\boldsymbol{\beta})\]
\[\mathbf{p} = \begin{bmatrix} p_1 \\ p_2 \\ \vdots \\ p_n \end{bmatrix} = \begin{bmatrix} \sigma(\boldsymbol{\beta}^T \mathbf{x}_1) \\ \sigma(\boldsymbol{\beta}^T \mathbf{x}_2) \\ \vdots \\ \sigma(\boldsymbol{\beta}^T \mathbf{x}_n) \end{bmatrix} = \begin{bmatrix} \frac{1}{1 + e^{-\boldsymbol{\beta}^T \mathbf{x}_1}} \\ \frac{1}{1 + e^{-\boldsymbol{\beta}^T \mathbf{x}_2}} \\ \vdots \\ \frac{1}{1 + e^{-\boldsymbol{\beta}^T \mathbf{x}_n}} \end{bmatrix}\]
This maps every linear prediction \(z_i\) into a valid probability \(p_i \in (0, 1)\).
To calculate the gradient of the log-likelihood, we need the derivative of the sigmoid function.
Fortunately, the sigmoid has a derivative that is easy to remember:
\[ \sigma'(z) = \sigma(z)\,\big(1 - \sigma(z)\big). \]
Using it, with \(p_i = \sigma(\boldsymbol{\beta}^T\mathbf{x}_i)\), where \(\mathbf{x}_i\) is the \(i^{th}\) row of the design matrix, and using the chain rule to break apart the derivative:
\[ \frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = \frac{\partial \ell_i}{\partial p_i} \cdot \frac{\partial p_i}{\partial z_i} \cdot \frac{\partial z_i}{\partial \boldsymbol{\beta}} \]
the partial derivative of the log-likelihood for a single observation is:
\[ \nabla \ell(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - p_i)\,\mathbf{x}_i = X^T(\mathbf{y} - \mathbf{p}), \]
where \(X\) is the \(n \times (p+1)\) matrix with rows \(\mathbf{x}_i^T\) and \(\mathbf{p} = \sigma(X\boldsymbol{\beta})\) applied entrywise.
Residual times feature – the same shape as the normal equations \(X^T(\mathbf{y} - X\boldsymbol{\beta}) = 0\) from linear regression, with a sigmoid inside. The difference is that \(\mathbf{p}\) depends on \(\boldsymbol{\beta}\) nonlinearly, so we cannot solve \(\nabla\ell = 0\) directly.
The notes has the full derivation of the gradient.
Here is the step-by-step derivation of the log-likelihood gradient for logistic regression.
The log-likelihood \(\ell(\boldsymbol{\beta})\) for \(n\) independent observations is:
\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \log(p_i) + (1 - y_i) \log(1 - p_i) \right]\]
where \(p_i = \sigma(\boldsymbol{\beta}^T \mathbf{x}_i) = \frac{1}{1 + e^{-\boldsymbol{\beta}^T \mathbf{x}_i}}\).
First, note the useful property of the sigmoid function \(\sigma(z)\):
\[\sigma'(z) = \sigma(z)(1 - \sigma(z))\]
Using this property along with the chain rule, the derivative of \(p_i\) with respect to \(\boldsymbol{\beta}\) is:
\[\frac{\partial p_i}{\partial \boldsymbol{\beta}} = \sigma'(\boldsymbol{\beta}^T \mathbf{x}_i) \cdot \mathbf{x}_i = p_i (1 - p_i) \mathbf{x}_i\]
Now differentiate the single-observation log-likelihood \(\ell_i(\boldsymbol{\beta}) = y_i \log(p_i) + (1 - y_i) \log(1 - p_i)\) with respect to \(\boldsymbol{\beta}\):
\[\frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = y_i \cdot \frac{1}{p_i} \frac{\partial p_i}{\partial \boldsymbol{\beta}} + (1 - y_i) \cdot \frac{1}{1 - p_i} \left( -\frac{\partial p_i}{\partial \boldsymbol{\beta}} \right)\]
Substitute \(\frac{\partial p_i}{\partial \boldsymbol{\beta}} = p_i (1 - p_i) \mathbf{x}_i\) into the equation:
\[\frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = y_i \cdot \frac{1}{p_i} \left[ p_i (1 - p_i) \mathbf{x}_i \right] - (1 - y_i) \cdot \frac{1}{1 - p_i} \left[ p_i (1 - p_i) \mathbf{x}_i \right]\]
Cancel \(p_i\) in the first term and \((1 - p_i)\) in the second term:
\[\frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = y_i (1 - p_i) \mathbf{x}_i - (1 - y_i) p_i \mathbf{x}_i\]
Factor out \(\mathbf{x}_i\):
\[\frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = \left( y_i - y_i p_i - p_i + y_i p_i \right) \mathbf{x}_i\]
The \(y_i p_i\) terms cancel out:
\[\frac{\partial \ell_i}{\partial \boldsymbol{\beta}} = (y_i - p_i) \mathbf{x}_i\]
Summing over all \(n\) observations yields the gradient of the full log-likelihood:
\[\nabla \ell(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - p_i) \mathbf{x}_i\]
In matrix form, this can be written compactly as:
\[\nabla \ell(\boldsymbol{\beta}) = X^T (\mathbf{y} - \mathbf{p})\]
where:
Algorithm:
Initialize parameters: Start with random values \(\boldsymbol{\beta}^{(0)}\)
Compute gradient: Calculate \(\nabla \ell(\boldsymbol{\beta}^{(t)}) = X^T(\mathbf{y} - \mathbf{p}^{(t)})\) – the direction of steepest increase
Update parameters: Take a step in that direction: \[\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + \eta \nabla \ell(\boldsymbol{\beta}^{(t)})\] where \(\eta\) is the learning rate (step size; the same \(\eta\) as in Lecture 16)
Repeat steps 2-3 until convergence (gradient \(\approx 0\) or max iterations reached)
Result: Parameters that maximize the likelihood of the observed data. The log-likelihood of logistic regression is concave, so this maximum is global.
We are climbing \(\ell\), so this is gradient ascent; the libraries minimize the log-loss instead, which is gradient descent on \(-\ell\) – the same steps. Gradient descent gets a full treatment in Neural Networks I.
Three observations, one feature: \(\mathbf{x} = [-1, 0, 2]\), \(\mathbf{y} = [0, 0, 1]\), so adding a column of 1s for the intercept gives
\[ \mathbf{X} = \begin{bmatrix} 1 & -1 \\ 1 & 0 \\ 1 & 2 \end{bmatrix} \]
Start at \(\boldsymbol{\beta}^{(0)} = [0, 0]^T\) with \(\eta = 0.1\).
We can evaluate \(\mathbf{p} = \sigma(\mathbf{X}\boldsymbol{\beta}^{(0)})\) and the gradient \(\nabla\ell(\boldsymbol{\beta}^{(0)})\) by hand.
The intercept moves down (two of three labels are 0) and the slope moves up (the 1 sits at the largest \(x\)). In today’s activity you will code exactly this loop and check that it converges to what sklearn reports.
So, in summary, we have:
Input pairs \((\mathbf{x}_i,y_i)\)
Output parameters \(\widehat{\boldsymbol{\beta}}\) that maximize the likelihood of the data given these parameters for the logistic regression model.
Method Maximum likelihood estimation, obtained by gradient ascent on \(\ell\) (the packages use a faster second-order variant of the same idea).
The standard package will give us a coefficient \(\beta_j\) for each independent variable (feature).
If we want to include a constant (i.e., \(\beta_0\)) we need to add a column of 1s (just like in linear regression).
df['intercept'] = 1.0
train_cols = df.columns[1:]
train_colsIndex(['gre', 'gpa', 'rank', 'intercept'], dtype='str')logit = sm.Logit(df['admit'], df[train_cols])
# fit the model
result = logit.fit()Optimization terminated successfully.
Current function value: 0.574302
Iterations 6Statsmodels gives us a summary of the model fit.
result.summary()| Dep. Variable: | admit | No. Observations: | 400 |
|---|---|---|---|
| Model: | Logit | Df Residuals: | 396 |
| Method: | MLE | Df Model: | 3 |
| Date: | Sat, 17 Oct 2026 | Pseudo R-squ.: | 0.08107 |
| Time: | 12:44:29 | Log-Likelihood: | -229.72 |
| converged: | True | LL-Null: | -249.99 |
| Covariance Type: | nonrobust | LLR p-value: | 8.207e-09 |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| gre | 0.0023 | 0.001 | 2.101 | 0.036 | 0.000 | 0.004 |
| gpa | 0.7770 | 0.327 | 2.373 | 0.018 | 0.135 | 1.419 |
| rank | -0.5600 | 0.127 | -4.405 | 0.000 | -0.809 | -0.311 |
| intercept | -3.4495 | 1.133 | -3.045 | 0.002 | -5.670 | -1.229 |
All three predictors are significant at the 5% level: no 95% confidence interval contains zero. The GRE interval is displayed as [0.000, 0.004] only because of rounding; unrounded it is
result.conf_int().loc['gre']0 0.000154
1 0.004434
Name: gre, dtype: float64The coefficients live on the log-odds scale, and \(e^{\beta_j}\) is the factor by which the odds are multiplied when \(x_j\) increases by one unit, holding the other features fixed.
coef_table = pd.DataFrame({
'coefficient': result.params,
'odds ratio e^beta': np.exp(result.params),
})
coef_table.loc['gre (per 100 points)'] = [100 * result.params['gre'],
np.exp(100 * result.params['gre'])]
coef_table.round(3)| coefficient | odds ratio e^beta | |
|---|---|---|
| gre | 0.002 | 1.002 |
| gpa | 0.777 | 2.175 |
| rank | -0.560 | 0.571 |
| intercept | -3.450 | 0.032 |
| gre (per 100 points) | 0.229 | 1.258 |
b0 = result.params['intercept']
mean_applicant = [df.gre.mean(), df.gpa.mean(), df['rank'].mean(), 1.0]
print(f'intercept beta_0 = {b0:.3f}; sigma(beta_0) = {1/(1+np.exp(-b0)):.3f}')
print(f'fraction admitted in the data = {df.admit.mean():.4f}')
print(f'p(y=1 | x) at the mean GRE, GPA, rank = {result.predict([mean_applicant])[0]:.3f}')intercept beta_0 = -3.450; sigma(beta_0) = 0.031
fraction admitted in the data = 0.3175
p(y=1 | x) at the mean GRE, GPA, rank = 0.298Note that by fitting a model to the data, we can make predictions for inputs that were not in the training data.
Furthermore, we can make a prediction of a probability for cases where we don’t have enough data to estimate the probability directly – e.g., for specific parameter values.
Let’s see how well the model fits the data.
We have three independent variables, so in each case we’ll use average values for the two that we aren’t evaluating.
GPA (GRE = 600, Rank = 2.5):
gpa_grid = np.linspace(df.gpa.min(), df.gpa.max(), 50)
prob = [result.predict([600, g, 2.5, 1.0]) for g in gpa_grid]
bins = np.linspace(df.gpa.min(), df.gpa.max(), 10)
groups = df.groupby(pd.cut(df.gpa, bins, include_lowest=True), observed=True)['admit'].mean()
centers = [interval.mid for interval in groups.index]
ax = plt.figure(figsize = (7, 4)).add_subplot()
ax.plot(gpa_grid, prob)
ax.plot(centers, groups, 'o')
ax.tick_params(labelsize=12)
ax.set_xlabel('gpa', fontsize = 14)
ax.set_ylabel('P[admit]', fontsize = 14)
ax.set_title('Marginal Effect of GPA', fontsize = 16);
GRE Score (GPA = 3.4, Rank = 2.5):
prob = [result.predict([b, 3.4, 2.5, 1.0]) for b in sorted(df.gre.unique())]
ax = plt.figure(figsize = (7, 4)).add_subplot()
ax.plot(sorted(df.gre.unique()), prob)
ax.plot(df.groupby('gre').mean()['admit'],'o')
ax.tick_params(labelsize=12)
ax.set_xlabel('gre', fontsize = 14)
ax.set_ylabel('P[admit]', fontsize = 14)
ax.set_title('Marginal Effect of GRE', fontsize = 16);
Institution Rank (GRE = 600, GPA = 3.4):
prob = [result.predict([600, 3.4, b, 1.0]) for b in range(1,5)]
ax = plt.figure(figsize = (7, 4)).add_subplot()
ax.plot(range(1,5), prob)
ax.plot(df.groupby('rank').mean()['admit'],'o')
ax.tick_params(labelsize=12)
ax.set_xlabel('Rank', fontsize = 14)
ax.set_xlim([0.5,4.5])
ax.set_ylabel('P[admit]', fontsize = 14)
ax.set_title('Marginal Effect of Rank', fontsize = 16);
At the start of lecture we emphasized that logistic regression is concerned with estimating a probability model for discrete (0/1) data.
However, it may well be the case that we want to do something with the probability that amounts to classification.
For example, we may classify data items using a rule such as “Assign item \(\mathbf{x}_i\) to Class 1 if \(p(y_i{=}1 \mid \mathbf{x}_i) > 0.5\)”.
For this reason, logistic regression could be considered a classification method.
Any other threshold \(t\) slides the hyperplane parallel to itself, to \(\boldsymbol{\beta}^T\mathbf{x} = \text{logit}(t)\).
Rank-2 applicants, with the boundary \(\hat{\boldsymbol{\beta}}^T\mathbf{x} = \text{logit}(t)\) drawn for \(t = 0.5\) and \(t = 0.3\) (solved for GPA as a function of GRE):
b = result.params
rank_shown = 2
gre_line = np.array([df.gre.min(), df.gre.max()])
def boundary_gpa(t):
# beta^T x = logit(t), solved for gpa with rank fixed
logit_t = np.log(t / (1 - t))
return (logit_t - b['intercept'] - b['rank'] * rank_shown - b['gre'] * gre_line) / b['gpa']
shown = df[df['rank'] == rank_shown]
colors = {0: '#56B4E9', 1: '#D55E00'} # colour-blind-safe (Okabe-Ito)
fig, ax = plt.subplots(figsize=(7, 4.5))
for label, grp in shown.groupby('admit'):
ax.scatter(grp.gre, grp.gpa, c=colors[label], s=35, alpha=0.85,
edgecolor='k', linewidth=0.3,
label='admitted' if label == 1 else 'rejected')
ax.plot(gre_line, boundary_gpa(0.5), 'k-', linewidth=2, label=r'$p(y{=}1 \mid \mathbf{x}) = 0.5$')
ax.plot(gre_line, boundary_gpa(0.3), 'k--', linewidth=2, label=r'$p(y{=}1 \mid \mathbf{x}) = 0.3$')
ax.set_ylim(df.gpa.min() - 0.1, df.gpa.max() + 0.1)
ax.set_xlabel('GRE', fontsize=13)
ax.set_ylabel('GPA', fontsize=13)
ax.set_title('Rank-2 applicants: the boundary is a line', fontsize=14)
ax.legend(loc='upper left', fontsize=10)
ax.tick_params(labelsize=11)
plt.tight_layout()
Above the line the model says “admit”. At \(t = 0.5\) almost no rank-2 applicant qualifies; lowering the threshold to \(0.3\) slides the same line down into the cloud. Changing the threshold never bends the boundary.
Let’s use our logistic regression as a classifier.
We want to ask whether we can correctly predict whether a student gets admitted to graduate school.
Let’s separate our training and test data:
X_train, X_test, y_train, y_test = model_selection.train_test_split(
df[train_cols], df['admit'],
test_size=0.4, random_state=1)Now, remember from precision and recall, we can consider the trade-offs.
Let’s say our classifier is outputting “yes” when it thinks the student will be admitted.
There are four cases:
Accuracy is the fraction of all classifications that are correct: \((\text{TP} + \text{TN}) / n\).
Accuracy has a floor set by the prevalence (the base rate).
Only 31.75% of these applicants were admitted, so a “classifier” that says no to everyone is 68% accurate without looking at a single feature.
Under class imbalance, accuracy mostly measures the imbalance – which is why we need the next two metrics.
Precision is the fraction of “yes” classifications that are correct:
\[ \mbox{Precision} = \frac{\mbox{True Positives}}{\mbox{True Positives + False Positives}}. \]
Recall is the fraction of admits that we say “yes” to:
\[ \mbox{Recall} = \frac{\mbox{True Positives}}{\mbox{True Positives + False Negatives}}. \]
def evaluate(y_train, X_train, y_test, X_test, threshold):
# learn model on training data
logit = sm.Logit(y_train, X_train)
result = logit.fit(disp=False)
# make probability predictions on test data
y_pred = result.predict(X_test)
# threshold probabilities to create classifications
y_pred = y_pred > threshold
# report metrics
precision = metrics.precision_score(y_test, y_pred, zero_division=0)
recall = metrics.recall_score(y_test, y_pred)
accuracy = metrics.accuracy_score(y_test, y_pred)
return precision, recall, accuracy
precision, recall, accuracy = evaluate(y_train, X_train, y_test, X_test, 0.5)
print(f'Precision: {precision:0.3f}, Recall: {recall:0.3f}, Accuracy: {accuracy:0.3f}')
print(f'Accuracy of "always reject" on this test set: {1 - y_test.mean():0.3f}')Precision: 0.586, Recall: 0.340, Accuracy: 0.719
Accuracy of "always reject" on this test set: 0.688The accuracy barely beats “always reject”, and the recall says why: at a 0.5 threshold the model says yes to only a small fraction of the real admits.
A single random split is noisy. Let’s average precision and recall over 20 different train/test splits:
PR = []
for i in range(20):
X_train, X_test, y_train, y_test = model_selection.train_test_split(
df[train_cols], df['admit'],
test_size=0.4, random_state=i)
PR.append(evaluate(y_train, X_train, y_test, X_test, 0.5))avgPrec = np.mean([f[0] for f in PR])
avgRec = np.mean([f[1] for f in PR])
avgAcc = np.mean([f[2] for f in PR])
print(f'Average Precision: {avgPrec:0.3f}, Average Recall: {avgRec:0.3f}, Average Accuracy: {avgAcc:0.3f}')Average Precision: 0.612, Average Recall: 0.258, Average Accuracy: 0.712Sometimes we would like a single value that describes the overall performance of the classifier.
For this, we take the harmonic mean of precision and recall, called F1 Score:
\[ \mbox{F1 Score} = \frac{2}{\mbox{Recall}^{-1} + \mbox{Precision}^{-1}}. = 2 \;\;\frac{\mbox{Precision} \cdot \mbox{Recall}}{\mbox{Precision} + \mbox{Recall}}. \]
Using this, we can evaluate other settings for the threshold.
def evalThresh(df, thresh):
PR = []
for i in range(20):
X_train, X_test, y_train, y_test = model_selection.train_test_split(
df[train_cols], df['admit'],
test_size=0.4, random_state=i)
PR.append(evaluate(y_train, X_train, y_test, X_test, thresh))
avgPrec = np.mean([f[0] for f in PR])
avgRec = np.mean([f[1] for f in PR])
# at high thresholds the model predicts no positives at all: F1 is 0, not 0/0
f1 = 0.0 if avgPrec + avgRec == 0 else 2 * (avgPrec * avgRec) / (avgPrec + avgRec)
return f1, avgPrec, avgRec
tvals = np.linspace(0.05, 0.8, 50)
f1vals = [evalThresh(df, tval)[0] for tval in tvals]plt.figure(figsize=(6, 3))
plt.plot(tvals,f1vals)
plt.ylabel('F1 Score')
plt.xlabel('Threshold for Classification')
plt.title('F1 as a function of Threshold');
Based on this plot, we can say that the best classification threshold appears to be around 0.35 – well below 0.5 – where precision and recall are:
F1, Prec, Rec = evalThresh(df, 0.35)
print('Precision: {:0.3f}, Recall: {:0.3f}, F1: {:0.3f}'.format(Prec, Rec, F1))Precision: 0.464, Recall: 0.618, F1: 0.530Compared with the 0.5 threshold we traded precision for recall: the model now says yes to most of the real admits, at the cost of more false alarms.
\(0.5\) is not a law of nature: it is where the model thinks the two classes are equally likely, which says nothing about what the two kinds of mistake cost.
If a false positive costs \(C_{FP}\) and a false negative costs \(C_{FN}\), saying “yes” is worth it exactly when \(p\,C_{FN} > (1-p)\,C_{FP}\), i.e. when
\[ p > t^* = \frac{C_{FP}}{C_{FP} + C_{FN}}. \]
A miss that is four times as costly as a false alarm puts the threshold at \(t^* = 0.2\). Choosing \(t\) is a decision about costs, and F1 is just one cost-free way to summarize the trade-off.
The example here is based on http://blog.yhathq.com/posts/logistic-regression-and-python.html where you can find additional details.
Two rates, both out of the true classes:
True positive rate (TPR, = recall): \(\text{TP}/(\text{TP}+\text{FN})\), the fraction of real admits we say yes to.
False positive rate (FPR): \(\text{FP}/(\text{FP}+\text{TN})\), the fraction of real rejects we wrongly say yes to.
The ROC curve plots TPR against FPR, one point per threshold. The diagonal is a random classifier.
The AUC is the area under the ROC curve: a threshold-free summary of how well the model ranks admits above rejects.
# Fit on a training split and score the held-out split -- a ROC curve on the
# training data would flatter the model.
X_train, X_test, y_train, y_test = model_selection.train_test_split(
df[train_cols], df['admit'], test_size=0.4, random_state=1)
roc_fit = sm.Logit(y_train, X_train).fit(disp=False)
y_score = roc_fit.predict(X_test)
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_score)
auc_score = metrics.roc_auc_score(y_test, y_score)
fig, ax = plt.subplots(figsize=(4, 4))
ax.plot(fpr, tpr, linewidth=2, label=f'ROC curve (AUC = {auc_score:.3f})')
ax.plot([0, 1], [0, 1], 'k--', linewidth=1, label='Random classifier')
ax.set_xlabel('False Positive Rate', fontsize=12)
ax.set_ylabel('True Positive Rate', fontsize=12)
ax.set_title('ROC curve on the held-out split', fontsize=13)
ax.legend(loc='lower right', fontsize=10)
ax.grid(True, alpha=0.3)
ax.tick_params(labelsize=11)
plt.tight_layout()
plt.show()
print(f"AUC Score: {auc_score:.4f}")
AUC Score: 0.7061class_weight="balanced".Shown here for 1 feature for simplicity.↩︎