Classification

Graduate Quantitative Economics and Datascience

Jesse Perla
jesse.perla@ubc.ca

University of British Columbia

Overview

Motivation

  • Least squares predicts a continuous \(y\). Classification predicts a discrete one: a category, a yes/no, a code
  • Supervised learning: every training observation carries a label \(y_n\). The next lecture is the unsupervised case, features with no labels
  • Economists use classifiers two ways: the outcome itself is discrete (labor force participation, loan approval), or the classifier builds a variable that then enters a regression (an occupation code from a job title, a topic from text)
  • The second use has an econometric consequence: the classifier’s mistakes are measurement error in the generated variable. Two sections near the end say what that does to a regression, the second of them optional with the formulas

Extra Materials

Packages

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.special import expit
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split
from sklearn.metrics import (confusion_matrix, ConfusionMatrixDisplay,
                             precision_score, recall_score, RocCurveDisplay)
from sklearn.feature_extraction.text import CountVectorizer
import statsmodels.api as sm
import statsmodels.formula.api as smf

From Least Squares to Logistic Regression

The Linear Probability Model

  • With \(y_n \in \{0, 1\}\), \({\mathbb{E}}[y_n \,|\, x_n] = \mathbb{P}(y_n = 1 \,|\, x_n)\), so why not least squares?

\[ \mathbb{P}(y_n = 1 \,|\, x_n) = x_n^{\top} \beta \]

  • Fine for a quick look at average effects, and \(\beta_j\) is directly the change in the probability
  • But fitted values can be below \(0\) or above \(1\), and a unit change in \(x_j\) has the same effect whether the probability starts at \(0.5\) or \(0.99\)
  • We want a model whose output is a probability, then a rule that turns it into a class

The Sigmoid

  • The sigmoid maps \({\mathbb{R}}\) to \((0, 1)\); in scipy it is expit

\[ \sigma(z) = \frac{1}{1 + e^{-z}} \]

  • Its inverse is the log-odds, \(\log \frac{p}{1 - p} = z\)
  • Logistic regression puts the linear index inside it

\[ \mathbb{P}(y_n = 1 \,|\, x_n) = \sigma(x_n^{\top} \beta) \]

  • \(\beta_j\) is the effect of \(x_j\) on the log-odds

Question What Probability Does the Model Assign?

A logistic regression gives \(x^{\top}\beta = 0\) for an observation.

What probability does it assign to \(y = 1\)?

A. \(0\)
B. \(1/2\)
C. \(1\)
D. It depends on the classification threshold

The Likelihood

  • Least squares finds \(\beta \in {\mathbb{R}}^M\) by \(\min_{\beta} \sum_{n=1}^N (y_n - x_n^{\top} \beta)^2\)
  • Logistic regression keeps the same \(x_n^{\top} \beta\) but changes the loss: each \(y_n\) is Bernoulli with \(p_n \equiv \sigma(x_n^{\top} \beta)\), so the MLE minimizes the negative log-likelihood \(-\ell(\beta)\), which divided by \(N\) ML calls the cross-entropy or log loss

\[ \min_{\beta} -\sum_{n=1}^N \left[ y_n \log \sigma(x_n^{\top} \beta) + (1 - y_n) \log\left(1 - \sigma(x_n^{\top} \beta)\right) \right] \]

  • It is convex for any \(X\) and \(y\), so any local minimum is the global one, and strictly convex when \(X\) has full column rank, so that minimum is unique if it exists. There is no closed form

Gradient, Hessian, and Newton’s Method

  • With \(\sigma'(z) = \sigma(z)(1 - \sigma(z))\),

\[ \nabla (-\ell)(\beta) = -\sum_{n=1}^N (y_n - p_n)\, x_n, \qquad \nabla^2 (-\ell)(\beta) = \sum_{n=1}^N p_n (1 - p_n)\, x_n x_n^{\top} \]

  • As in least squares the gradient is residual times \(x_n\), but \(p_n\) is nonlinear in \(\beta\)
  • The Hessian is \(X^{\top} W X\) with \(W\) diagonal with entries \(p_n (1 - p_n)\), positive definite when \(X\) has full column rank, so \(-\ell\) is strictly convex and any minimum is the only one (recall the definiteness lecture)
  • No closed form, but Newton’s method is a sequence of weighted least squares problems

\[ \beta_{k+1} = \beta_k + \left(X^{\top} W_k X\right)^{-1} X^{\top} (y - p_k) \]

Perfect Separation is Interpolation

  • Perfect separation: some \(\beta\) has \(x_n^{\top} \beta > 0\) for every \(y_n = 1\) and \(x_n^{\top} \beta < 0\) for every \(y_n = 0\), so the training data are classified with zero errors
  • Scaling it by \(s \to \infty\) sends \(\sigma(s\, x_n^{\top} \beta) \to y_n\) for every \(n\), so the loss falls toward \(0\) but no finite \(\beta\) attains it: the MLE does not exist
  • The classification version of interpolation: as with least squares when \(M > N\), the data are fit exactly and do not pick a unique \(\beta\). There we chose the min-norm solution
  • Common with many features and few observations

Question Double the Coefficients

Finite coefficients perfectly separate the training data. Double the intercept and slopes:

\[ \begin{aligned} p_n &= \sigma(\beta_0 + x_n^{\top}\beta),\\ \widetilde p_n &= \sigma\!\left(2\beta_0 + x_n^{\top}(2\beta)\right). \end{aligned} \]

At threshold \(1/2\), what changes?

A. Same classes; probabilities nearer \(0\) or \(1\).
B. Some classes change.
C. Same classes and probabilities.
D. Probabilities nearer \(1/2\).

Perfect Separation: Which Line?

  • Every line between the clusters separates the data, so the candidates are a whole cone of directions, not a point
  • The solid line is the furthest from both clusters, the widest margin. The limit of vanishing regularization picks it

Regularization and scikit-learn’s C

  • LogisticRegression does not fit the MLE by default. It adds a ridge penalty to the same objective

\[ \min_{\beta_0, \beta} \left\{ -\sum_{n=1}^N \left[ y_n \log \sigma(\beta_0 + x_n^{\top} \beta) + (1 - y_n) \log\left(1 - \sigma(\beta_0 + x_n^{\top} \beta)\right) \right] + \frac{1}{2C} \|\beta\|_2^2 \right\} \]

  • The penalty is on the slopes only (the intercept \(\beta_0\) is not penalized): smaller C, more shrinkage (default C=1.0). It fixes perfect separation and reduces overfitting at the cost of biased coefficients
  • C=np.inf turns the penalty off and matches statsmodels.Logit; l1_ratio=1.0, solver="saga" swaps in the lasso penalty \(\|\beta\|_1\), which sets coefficients exactly to zero (the default lbfgs solver only handles ridge)

Regularization Picks the Widest Margin

  • With separated data every separating direction drives the loss to \(0\), so the penalty decides: the fit wants the one with the smallest \(\|\beta\|_2\)
  • Scaling \((\beta_0, \beta)\) does not move the boundary \(\beta_0 + x^{\top} \beta = 0\), so fix the scale by requiring every point to clear it by at least \(1\)

\[ \min_{\beta_0, \beta} \|\beta\|_2 \,\text{s.t.}\,\beta_0 + x_n^{\top} \beta \geq 1 \text{ if } y_n = 1, \qquad \beta_0 + x_n^{\top} \beta \leq -1 \text{ if } y_n = 0 \]

  • The distance from \(x_n\) to the boundary is \(|\beta_0 + x_n^{\top} \beta| / \|\beta\|_2\), so the nearest point sits at \(1/\|\beta\|_2\) and the smallest norm is the widest margin: the min-norm solution again
  • As \(C \to \infty\) the ridge-penalized logistic direction converges to it (Rosset, Zhu, and Hastie 2004)
  • Any finite C lands in between, with shrunk coefficients as in ridge regression

Simulated Data and Two Fits

rng = np.random.default_rng(526)
N = 500
X = rng.normal(size=(N, 2))
beta_true = np.array([1.5, -1.0])
p = expit(0.3 + X @ beta_true)
y = (rng.random(N) < p).astype(int)
clf = LogisticRegression(C=np.inf).fit(X, y)
print("sklearn:    ", clf.intercept_, clf.coef_[0])
logit = sm.Logit(y, sm.add_constant(X)).fit(disp=0)
print("statsmodels:", logit.params)
sklearn:     [0.33244846] [ 1.43336144 -1.09263345]
statsmodels: [ 0.33279073  1.43333638 -1.09242626]
  • Identical coefficients because C=np.inf turns the penalty off. Econometrics reads \(\beta\) and marginal effects (logit.get_margeff()); ML reads predicted classes

Predicted Classes and Probabilities

print(clf.predict_proba(X[:3]).round(3))
print(clf.predict(X[:10]))
print(f"accuracy on the training data: {clf.score(X, y):.3f}")
[[0.846 0.154]
 [0.113 0.887]
 [0.817 0.183]]
[0 1 0 0 1 1 1 0 1 0]
accuracy on the training data: 0.752
  • predict_proba gives \(\left[\mathbb{P}(y = 0 \,|\, x), \mathbb{P}(y = 1 \,|\, x)\right]\); predict applies a threshold at \(1/2\)

The Decision Boundary is Linear

  • \(\sigma(x^{\top} \beta) > 1/2\) exactly when \(x^{\top} \beta > 0\), so the boundary \(\{x : x^{\top} \beta = 0\}\) is a hyperplane
  • Nonlinear boundaries come from nonlinear features (polynomials, interactions, embeddings), as in least squares
  • Points near the boundary have \(p \approx 1/2\): that is where the mistakes are

Decision Rules and Evaluation

The Confusion Matrix

  • Cross-tabulate the true class against the predicted class, on held-out data in practice; the examples below reuse the training sample
  • Rows are the truth, columns the prediction (the sklearn convention)
  • Two kinds of mistakes with different costs: a false positive (predicted \(1\), truly \(0\)) and a false negative (predicted \(0\), truly \(1\))
  • The same matrix returns for k-means next lecture, where cluster labels are arbitrary

Accuracy, Precision, and Recall

Each metric answers a different question about the confusion matrix.

Metric Calculation Interpretation
Accuracy \((TP + TN) / N\) Among all observations, what share did we classify correctly?
Precision \(TP / (TP + FP)\) Among predicted positives (\(\hat y = 1\)), what share are actually positive?
Recall \(TP / (TP + FN)\) Among actual positives (\(y = 1\)), what share did we find?
  • Precision and recall count the same true positives, but divide by different groups.
  • Precision conditions on the prediction; recall conditions on the truth.

Error Rates and Metrics in Code

  • False positive rate: among actual negatives (\(y = 0\)), what share did we flag?

\[ p_0 = \frac{FP}{FP + TN} \]

  • False negative rate: among actual positives (\(y = 1\)), what share did we miss?

\[ p_1 = \frac{FN}{TP + FN} = 1 - \text{recall} \]

y_hat = clf.predict(X)
print(confusion_matrix(y, y_hat))
print(f"accuracy  {clf.score(X, y):.3f}")
print(f"precision {precision_score(y, y_hat):.3f}")
print(f"recall    {recall_score(y, y_hat):.3f}")
[[158  67]
 [ 57 218]]
accuracy  0.752
precision 0.765
recall    0.793

Rare Classes Break Accuracy

  • With \(5\%\) positives, always predicting \(0\) is \(95\%\) accurate and useless
p_rare = expit(-4.0 + X @ beta_true)
y_rare = (rng.random(N) < p_rare).astype(int)
clf_rare = LogisticRegression(C=np.inf).fit(X, y_rare)
y_hat_rare = clf_rare.predict(X)
print(f"share y = 1: {y_rare.mean():.3f}, always-zero accuracy: "
      f"{1 - y_rare.mean():.3f}, classifier accuracy: "
      f"{clf_rare.score(X, y_rare):.3f}")
print(confusion_matrix(y_rare, y_hat_rare))
print(f"recall {recall_score(y_rare, y_hat_rare):.3f}")
share y = 1: 0.048, always-zero accuracy: 0.952, classifier accuracy: 0.950
[[473   3]
 [ 22   2]]
recall 0.083
  • Report precision and recall whenever a class is rare, which is most variables economists want to code: a strike, a discrimination complaint, a recession

The Threshold is a Cost Decision

  • predict uses \(\hat{y} = \mathbf{1}\left\{{\hat{p}(x) > 1/2}\right\}\). With a false positive costing \(c_{FP}\) and a false negative \(c_{FN}\), the cheaper prediction is \(\hat{y} = 1\) when

\[ \hat{p}(x) > \frac{c_{FP}}{c_{FP} + c_{FN}} \equiv \tau \]

  • Equal costs give \(\tau = 1/2\); if a miss is nine times worse than a false alarm, \(\tau = 0.1\)
  • \(\hat{p}(x)\) estimates the posterior \(\mathbb{P}(y = 1 \,|\, x)\) of the probability lecture. Logistic regression models it directly; generative classifiers such as naive Bayes model \(\mathbb{P}(x \,|\, y)\) and \(\mathbb{P}(y)\) and use Bayes’ law
  • Lowering \(\tau\) raises recall at the expense of false alarms; class_weight="balanced" pushes the fit the same way

Question Which Prediction Minimizes Expected Cost?

For one observation, \(\mathbb{P}(y = 1 \,|\, x) = 0.2\).

A false positive costs 1, a false negative costs 5, and correct predictions cost nothing.

Which prediction minimizes expected cost?

A. Predict \(0\), because it is more likely.
B. Predict \(1\), because missing a positive is sufficiently costly.
C. Both predictions have the same expected cost.
D. We must refit the model before deciding.

ROC Curves

The ROC Curve: Intuition

  • Imagine a radar system: \(y = 1\) means an aircraft is present; \(y = 0\) means none is present. Predicting \(1\) means raise an alarm
  • The model estimates the probability an aircraft is present, given the radar readings
  • A low cutoff is cautious about missing an aircraft: raise alarms more readily, detecting more aircraft but accepting more false alarms
  • A high cutoff is cautious about false alarms: require stronger evidence before sounding an alarm, missing more aircraft
  • The receiver operating characteristic (ROC) curve shows this tradeoff for the same fitted model across all cutoffs:
    • Horizontal: false alarm rate (false positive rate)
    • Vertical: detection rate (true positive rate, or recall)

The ROC Curve: What Exactly Is \(\tau\)?

Fix the model’s estimated probabilities \(\hat p_n = \hat p(x_n)\) of an aircraft being present. The threshold \(\tau \in [0,1]\) sets the decision:

\[ \hat y_n(\tau) = \begin{cases} 1 \text{ (alarm)} & \text{if } \hat p_n > \tau,\\ 0 \text{ (no alarm)} & \text{if } \hat p_n \leq \tau. \end{cases} \]

  • If \(\hat p_n = 0.2\), we raise an alarm at \(\tau = 0.1\), but not at \(\tau = 0.5\)
  • Each cutoff gives one point: FPR horizontally and TPR vertically

\[ \text{FPR}(\tau) = \frac{FP(\tau)}{FP(\tau) + TN(\tau)}, \qquad \text{TPR}(\tau) = \frac{TP(\tau)}{TP(\tau) + FN(\tau)}. \]

  • \(\tau\) selects a point on the curve; it is not an axis. The probabilities stay fixed

The ROC Curve

  • Marked point: \(\tau = 1/2\) means raise an alarm if \(\hat p_n > 1/2\). On this schematic curve, the false alarm rate is \(12\%\) and the detection rate is \(72\%\)
  • Lower \(\tau\): both rates weakly increase. With sigmoid scores in \((0,1)\), \(\tau = 1\) never raises an alarm at \((0,0)\); \(\tau = 0\) always raises an alarm at \((1,1)\)
  • The area under the curve (AUC) summarizes the ranking across thresholds: the probability a random positive scores above a random negative, with half credit for ties

Question Which Model Is “Better”?

Financial transactions:
\(y = 1\): fraudulent; \(y = 0\): legitimate.

Which model is better? Explain using AUC.

ROC in Code: Ranking Versus Threshold

Share of \(y = 1\) in each dataset: balanced \(55\%\); rare \(4.8\%\).

fig, ax = plt.subplots(figsize=(4.5, 3.4))
s_bal = clf.predict_proba(X)[:, 1]
s_rare = clf_rare.predict_proba(X)[:, 1]
RocCurveDisplay.from_predictions(
    y, s_bal, ax=ax, name="balanced")
RocCurveDisplay.from_predictions(
    y_rare, s_rare, ax=ax, name="rare class")
plt.show()
for tau in [0.5, 0.1]:
    y_hat = (s_rare > tau).astype(int)
    print(f"tau = {tau}: recall "
          f"{recall_score(y_rare, y_hat):.2f}")

tau = 0.5: recall 0.08
tau = 0.1: recall 0.71
  • Lower \(\tau\): recall rises from \(8\%\) to \(71\%\); scores and ROC stay fixed
  • Rare positives: low FPR need not mean high precision

Generalization

LASSO versus Ridge

Minimize the negative log-likelihood plus a penalty on the slopes. The intercept is not penalized.

LASSO: \(L^1\) penalty

\[ \lambda \|\beta\|_1 = \lambda \sum_{j=1}^M |\beta_j| \]

  • Encourages many exact zeros: selects features

Ridge: squared \(L^2\) penalty

\[ \frac{\lambda}{2} \|\beta\|_2^2 = \frac{\lambda}{2} \sum_{j=1}^M \beta_j^2 \]

  • Encourages small coefficients, usually nonzero
  • Larger \(\lambda\): stronger penalty; \(\lambda = 0\): no penalty
  • With the summed loss above, \(\lambda = 1/C\): smaller C means more regularization

Train, Test, and Overfitting

  • Fit on training data; evaluate on untouched test data. Tune C on separate validation data
  • All features: \(2\) true features plus \(200\) noise features
X_noise = np.hstack([X, rng.normal(size=(N, 200))])
X_tr, X_te, y_tr, y_te = train_test_split(X_noise, y, test_size=0.5,
                                          random_state=526)
mle_all = LogisticRegression(C=np.inf, max_iter=5000).fit(X_tr, y_tr)
mle_true = LogisticRegression(C=np.inf).fit(X_tr[:, :2], y_tr)
ridge = LogisticRegression(l1_ratio=0, C=0.1).fit(X_tr, y_tr)
Unpenalized, all features  train 1.000  test 0.640
Unpenalized, true features train 0.732  test 0.740
Ridge, all features        train 0.992  test 0.624
  • With \(202\) features and \(250\) training rows the data are separated, so the unpenalized fit has no MLE: its coefficients are wherever the optimizer stopped
  • Ridge penalizes coefficient size; it does not guarantee higher test accuracy

More Than Two Classes

Multinomial Logit, the Softmax, and Random Utility

  • With \(K\) classes, give each its own \(\beta_k \in {\mathbb{R}}^M\) and apply the softmax

\[ \mathbb{P}(y = k \,|\, x) = \frac{e^{x^{\top} \beta_k}}{\sum_{j=1}^K e^{x^{\top} \beta_j}} \]

  • Only differences \(\beta_k - \beta_j\) matter, so one class can be normalized to \(\beta_1 = 0\) (scikit-learn keeps all \(K\) vectors and lets the penalty pin them down); with \(K = 2\) this is the sigmoid. The log-likelihood is still concave, and LogisticRegression fits it automatically when y has more than two values
  • Economists reached the same formula through discrete choice (McFadden 1974): an agent picks the option with the highest utility \(U_k = x^{\top} \beta_k + \epsilon_k\) with i.i.d. Gumbel \(\epsilon_k\), and \(\mathbb{P}(\text{choose } k \,|\, x)\) is the softmax
  • The same map from a score per option to a distribution is how a language model picks the next token, with \(K\) the vocabulary size

Application: Coding Occupations from Job Titles

The Problem

  • Map free-text titles such as POLICE OFFICER and SR DATA ANALYST to Standard Occupational Classification (SOC) major groups
  • O*NET supplies \(8{,}189\) titles already labelled with SOC codes. Their first two digits identify the major group; \(22\) groups appear in these data
  • Train and evaluate: fit the word-count classifier on \(75\%\) of these labelled titles; compare predictions with the known labels on the remaining \(25\%\)
  • Apply: classify Chicago payroll titles, which have no supplied SOC labels

Bag of Words

  • The simplest feature map counts how often each vocabulary word appears, ignoring order. CountVectorizer builds the vocabulary from the training titles and returns a sparse matrix with one column per word
demo = CountVectorizer().fit(["Police Officer", "Police Sergeant",
                              "Chief Officer"])
print(demo.get_feature_names_out())
print(demo.transform(["Police Officer", "Police Sergeant"]).toarray())
['chief' 'officer' 'police' 'sergeant']
[[0 1 1 0]
 [0 0 1 1]]
  • Words never seen in training are dropped at prediction time. The text lecture adds weighting and embeddings; counting is enough here

Loading the Data

  • The data files live in slides/data/ of the course repository. The notebooks read a local copy if there is one and the course website otherwise
def data_path(name):
    url = "https://jlperla.github.io/grad_econ_datascience/slides/data/"
    return f"data/{name}" if os.path.exists(f"data/{name}") else url + name

onet = pd.read_csv(data_path("onet_titles.csv.gz"))
sample = onet.sample(4, random_state=1)
sample = sample[["title", "soc_code", "major_group_title"]]
print(sample.to_string(index=False))
print(onet["major_group_title"].value_counts().head(4))
                                      title   soc_code                      major_group_title
              MHA (Mental Health Assistant) 29-2053.00 Healthcare Practitioners and Technical
                            Line Supervisor 51-1011.00                             Production
Home Care Physical Therapist (Home Care PT) 29-1123.00 Healthcare Practitioners and Technical
                  Licensed Reactor Operator 51-8011.00                             Production
major_group_title
Production                                1016
Healthcare Practitioners and Technical     734
Educational Instruction and Library        613
Construction and Extraction                534
Name: count, dtype: int64

Fitting the Coder

  • Word counts of the title, the major group as the label (\(22\) classes), a quarter of the titles held out
titles_tr, titles_te, y_tr, y_te = train_test_split(
    onet["title"], onet["major_group_title"], test_size=0.25, random_state=526)
vec = CountVectorizer()
X_tr = vec.fit_transform(titles_tr)
X_te = vec.transform(titles_te)
coder = LogisticRegression(solver="newton-cg").fit(X_tr, y_tr)
print(f"vocabulary size {X_tr.shape[1]}, classes {len(coder.classes_)}")
print(f"train accuracy {coder.score(X_tr, y_tr):.3f}, "
      f"test accuracy {coder.score(X_te, y_te):.3f}")
vocabulary size 2954, classes 22
train accuracy 0.861, test accuracy 0.719
  • newton-cg is Newton’s method from earlier with each step solved by conjugate gradient: about ten iterations, where the default lbfgs takes sixty. Train well above test: overfitting again

Where Does it Go Wrong?

top = onet["major_group_title"].value_counts().index[:6]
ConfusionMatrixDisplay.from_predictions(
    y_te, coder.predict(X_te), labels=top,
    display_labels=[t[:22] for t in top],
    xticks_rotation=40, colorbar=False)
plt.tight_layout()
plt.show()

  • The six largest groups on held-out titles, a submatrix that drops the titles predicted outside them. The off-diagonal mistakes are between groups whose titles share words

Coding a Payroll File

  • Apply the coder to \(8{,}000\) full-time salaried positions from the City of Chicago’s payroll
chicago = pd.read_csv(data_path("chicago_payroll.csv.gz"))
X_chicago = vec.transform(chicago["job_title"])
chicago["group"] = coder.predict(X_chicago)
chicago["p_max"] = coder.predict_proba(X_chicago).max(axis=1).round(2)
counts = chicago.groupby(["job_title", "group"])["p_max"]
counts = counts.agg(n="size", p_max="first")
print(counts.sort_values("n", ascending=False).head(8).to_string())
                                                                                  n  p_max
job_title                              group                                              
POLICE OFFICER                         Protective Service                      2465   0.89
FIREFIGHTER-EMT                        Healthcare Practitioners and Technical   550   0.18
SERGEANT                               Production                               447   0.15
POLICE OFFICER (ASSIGNED AS DETECTIVE) Protective Service                       376   0.99
FIREFIGHTER-EMT (RECRUIT)              Healthcare Practitioners and Technical   171   0.18
LIEUTENANT-EMT                         Protective Service                       144   0.13
FIRE ENGINEER-EMT                      Architecture and Engineering             141   0.75
PARAMEDIC                              Healthcare Practitioners and Technical   132   0.22
  • Chicago’s titles use a different dialect: FIREFIGHTER-EMT lands in Healthcare Practitioners, since Firefighter Paramedic is one of only three O*NET titles with the word and EMT appears in none; a title with no known word gets the largest-intercept group. Check the common titles and p_max before using the codes

The Regression the Codes Were Built For

premia = smf.ols("np.log(annual_salary) ~ C(group) + C(department)",
                 data=chicago).fit()
coef = premia.params.filter(like="group").round(2)
coef.index = coef.index.str.replace(r"C\(group\)\[T\.(.*)\]", r"\1", regex=True)
print(pd.concat([coef.nsmallest(4), coef.nlargest(4)]).to_string())
Building and Grounds Cleaning and Maintenance   -0.38
Healthcare Support                              -0.29
Healthcare Practitioners and Technical          -0.20
Office and Administrative Support               -0.18
Personal Care and Service                        0.26
Educational Instruction and Library              0.19
Management                                       0.19
Production                                       0.19
  • Log-salary differences relative to Architecture and Engineering within a department, the four lowest and four highest
  • Plausible magnitudes, but the regressor is a classifier’s output with about \(28\%\) of titles miscoded on held-out data: measurement error in a regressor. What does it do to the coefficients?

Predictions as Measurements

Classification Errors in a Salary Regression

  • Two equally sized groups: \(d_n = 1\) for managers; \(d_n = 0\) for other workers

\[ \log(\text{salary}_n) = \alpha + \beta d_n + u_n, \qquad {\mathbb{E}}[u_n \,|\, d_n] = 0 \]

  • \(\beta\) is the true gap in mean log salary between the two groups
  • A classifier mislabels \(25\%\) of each group, independently of salary within each true group
  • Regress log salary on the predicted group \(\hat d_n\) instead of the true \(d_n\):

\[ \underbrace{\beta = 0.2}_{\text{true salary-gap coefficient}} \quad\longrightarrow\quad \underbrace{\hat\beta \approx 0.1}_{\text{estimated coefficient in a large sample}} \]

  • Attenuation: bias in \(\hat\beta\) toward zero, even with a large sample

What to Do About It

  • Keep a random labelled sample, report its confusion matrix with the regression, and never treat classifier output as truth
  • Check the error rates by subgroup before trusting a correction; the text lecture does this on Newswire
  • Bias corrections use labelled validation data and assumptions about the errors. The optional section gives examples and introduces ValidMLInference (Battaglia, Christensen, Hansen, and Sacher, 2024)

Measurement Error from a ClassifierOptional

A Mismeasured Regressor (Aigner 1973)

  • The true \(d_n \in \{0, 1\}\) enters \(y_n = \alpha + \beta\, d_n + u_n\), but we observe the classifier’s \(\hat{d}_n\), with the error rates of its confusion matrix

\[ p_0 \equiv \mathbb{P}(\hat{d}_n = 1 \,|\, d_n = 0), \qquad p_1 \equiv \mathbb{P}(\hat{d}_n = 0 \,|\, d_n = 1) \]

  • If the mistakes are non-differential (given \(d_n\), independent of \(u_n\) and the other regressors), OLS on \(\hat{d}_n\) converges to

\[ \hat{\beta} \approx \beta \cdot \lambda, \qquad \lambda \approx 1 - p_0 - p_1 \]

  • For any classifier better than a coin flip (\(p_0 + p_1 < 1\)), \(\lambda\) is between \(0\) and \(1\): attenuation toward zero, never a sign flip, unchanged as \(N\) grows

Monte Carlo

def attenuation(N, p0, p1, pi=0.3, beta=1.0, reps=100):
    slopes = []
    for _ in range(reps):
        d = rng.random(N) < pi
        y = 2.0 + beta * d + rng.normal(size=N)
        d_hat = np.where(d, rng.random(N) > p1, rng.random(N) < p0)
        d_hat = d_hat.astype(float)
        slopes.append(sm.OLS(y, sm.add_constant(d_hat)).fit().params[1])
    return np.mean(slopes)

for N in [1_000, 10_000, 100_000]:
    print(f"N = {N:>7,}: mean slope {attenuation(N, 0.05, 0.20):.3f}")
pi, p0, p1 = 0.3, 0.05, 0.20
q = pi * (1 - p1) + (1 - pi) * p0
factor = (1 - p0 - p1) * pi * (1 - pi) / (q * (1 - q))
print(f"1 - p0 - p1 = {1 - p0 - p1:.3f}, exact factor {factor:.3f}")
N =   1,000: mean slope 0.789
N =  10,000: mean slope 0.789
N = 100,000: mean slope 0.789
1 - p0 - p1 = 0.750, exact factor 0.790
  • True \(\beta = 1\). The estimate sits at the exact factor \((1 - p_0 - p_1)\, \pi(1 - \pi) / (q(1 - q))\), with \(\pi\) the true share of ones and \(q\) the predicted share, and does not move with \(N\)

A Mismeasured Outcome

  • Now the classifier’s output is the outcome: the regression wants \(y_n \in \{0, 1\}\) and gets \(\tilde{y}_n\), with a false positive rate \(p_0\) and a false negative rate \(p_1\) that do not depend on \(x_n\)

\[ {\mathbb{E}}[\tilde{y}_n \,|\, x_n] = p_0 + (1 - p_0 - p_1)\, {\mathbb{E}}[y_n \,|\, x_n] \]

  • In a linear regression on \(\tilde{y}_n\) every coefficient shrinks by exactly the factor \(1 - p_0 - p_1\) (the regressor case had an approximation), and the intercept also picks up \(p_0\). Dividing the slopes by the factor undoes it
  • The text lecture applies the division to a coverage regression on Newswire, where the outcome is a topic classifier’s label

Differential Errors

  • Both results assume error rates that do not depend on the regressors. Error rates that differ by subgroup break it: a region whose images look unlike the training set, a season, the wording of a prompt. A different base rate alone does not, since \(p_0\) and \(p_1\) condition on the truth
  • Dividing by \(1 - p_0 - p_1\) can then overshoot: differential errors bias in either direction, not only toward zero
  • Diagnose it with the confusion matrix by subgroup, so the labelled sample has to cover the subgroups the regression compares

Corrections Exist

  • The held-out confusion matrix estimates \(p_0\) and \(p_1\), so the attenuation can be undone. Three to know about
    • Battaglia, Christensen, Hansen, and Sacher (2024): correct the moments with the validation-sample error rates; package ValidMLInference
    • Prediction-powered inference (Angelopoulos et al. 2023): combine predictions on the full sample with the labelled sample into valid inference, exactly unbiased for a mean
    • Design-based supervised learning (Egami, Hinck, Stewart, and Wei 2023): the same idea for text, with the labelled sample drawn at random from the corpus
  • Each wants a random labelled sample from the same population as the regression (Battaglia et al. can also use error rates reported elsewhere): it is the direct evidence on the classifier’s errors, and the only place differential error can show up

Review QuestionsQuestions

Question Same Classes, Same Probabilities?

Two logistic regressions both classify every observation in the same training sample correctly, using a threshold of \(1/2\).

Neither model assigns probability exactly \(1/2\) to any training observation.

Must the models assign the same probability of \(y_n = 1\) to each \(x_n\)?

A. Yes: both models classify every training observation correctly.
B. No: their probabilities only need to fall on the correct side of \(1/2\).

Question Precision and Recall for Fraud Detection

We predict whether financial transactions are fraudulent (\(y = 1\)) or legitimate (\(y = 0\)).

In a held-out sample:

  • Actually fraudulent: 20.
  • Flagged as fraudulent: 10.
  • Fraudulent among those flagged: 8.

What are its precision and recall?

A. Precision \(40\% = 8/20\); recall \(80\% = 8/10\).
B. Precision \(80\% = 8/10\); recall \(80\% = 8/10\).
C. Precision \(80\% = 8/10\); recall \(40\% = 8/20\).
D. Precision \(40\% = 8/20\); recall \(40\% = 8/20\).

Question LASSO or Ridge?

For \(\lambda > 0\), identify \(P_1\) and what increasing \(\lambda\) especially encourages.

\[ P_1(\beta) = \lambda \sum_{j=1}^M |\beta_j|, \qquad P_2(\beta) = \frac{\lambda}{2} \sum_{j=1}^M \beta_j^2. \]

A. LASSO; many exact zeros.
B. LASSO; small, usually nonzero coefficients.
C. Ridge; many exact zeros.
D. Ridge; small, usually nonzero coefficients.