Lab 5 — Logistic regression on a messy credit file

Author

Hicham Janati — FSML

A bank wants to rank its borrowers by risk: who will be in financial distress in the next two years. It gives you a file of 40 000 borrowers. The file is inspired by a public dataset, the Kaggle competition Give Me Some Credit, but it is not identical to it.

This is our first classifier in code: a logistic regression, scored with the metrics of the lecture. We keep a quarter of the borrowers aside to test the models, and every model gets two scores: on its training rows and on the test rows.

Boxes as always: 🔮 predict before running (no AI) · 🧠 explain in your own words · ✍️ short task, boilerplate given.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (confusion_matrix, ConfusionMatrixDisplay, precision_score,
                             recall_score, f1_score, roc_curve, roc_auc_score,
                             precision_recall_curve, average_precision_score)

BLUE, RED, GREEN, ORANGE, GRAY = "#1f6fb2", "#c92a1a", "#1c7a43", "#e8871a", "#8a8a8a"
plt.rcParams["figure.figsize"] = (6.5, 3.6)
plt.rcParams["axes.spines.top"] = plt.rcParams["axes.spines.right"] = False

Part 1 — The file and its missing values

1.1 The data and a first fit

df = pd.read_csv("https://hichamjanati.github.io/media/teaching/SMV/borrowers.csv")
print(df.shape)
df.head()

Each row is a borrower. What the columns mean:

column	meaning
`FinancialDistressNextTwoYears`	the target: 1 if the borrower was 90 days or more late on a payment in the next two years, 0 otherwise
`BorrowedToLimitRatio`	total amount borrowed (still owed) on credit cards and personal lines of credit, divided by the sum of their credit limits
`Age`	age of the borrower, in years
`NumberOfTimes90DaysLate`	number of times the borrower was 90 days or more late on a payment
`DebtToIncomeRatio`	monthly debt payments, alimony and living costs, divided by the monthly income; missing when the income is missing
`MonthlyIncome`	monthly income
`NumberOfOpenCreditLinesAndLoans`	number of open loans (car loan, mortgage) and lines of credit (credit cards)
`NumberOfRealEstateLoansOrLines`	number of mortgage and real estate loans, home equity lines of credit included

y = df["FinancialDistressNextTwoYears"]
X = df.drop(columns="FinancialDistressNextTwoYears")
print(f"{y.mean():.1%} of the borrowers are in distress: {y.sum()} of {len(y)}")

The two classes are imbalanced: fewer than 7 borrowers in 100 are in distress. As in Lab 4, train_test_split keeps a test set aside; with stratify=y, both parts get the same share of borrowers in distress.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25,
                                                    stratify=y, random_state=0)
print(f"train {len(X_train)} borrowers ({y_train.sum()} in distress) · "
      f"test {len(X_test)} borrowers ({y_test.sum()} in distress)")

LogisticRegression is the logistic regression of the lecture. It has the usual .fit, then:

.predict_proba(X) returns the predicted probabilities (we look at them in 1.2);
.predict(X) answers 1 when the probability of distress is above 0.5.

By default it carries the L2 penalty of the lecture, with C=1. C is the inverse of the penalty strength: a smaller C shrinks the coefficients more, and C=np.inf removes the penalty.

LogisticRegression().fit(X_train, y_train)

It refuses missing values. isna() marks them, and the mean of each column gives their share:

X_train.isna().mean().round(3)

The monthly income is missing for 33.0 % of the training borrowers, and with it the debt-to-income ratio. Two simple ways out: drop the borrowers with a missing value, or fill the missing values in.

1.2 Dropping the incomplete rows

X_train.notna() is True where a value is present. Run it:

X_train.notna().head()

.all(axis=0) gives, for each column, whether all its values are True; .all(axis=1) gives, for each row, whether all its values are True.

✍️ Your turn. Which axis gives the borrowers with no missing value? Use it to build complete, True for those training borrowers.

# TODO: complete = ...
complete.head()

✍️ Your turn. Fit drop_model on the complete borrowers only: a pipeline that standardizes the columns, then fits a logistic regression, with the steps named “scale” and “model”. Then complete the print: how many training borrowers are kept, and how many of them are in distress.

# TODO: drop_model = ...
# TODO: drop_model.fit(...)

# TODO: n_kept = the number of training borrowers kept
# TODO: n_kept_distress = the number of kept borrowers in distress
print(f"kept {n_kept} of {len(X_train)} training borrowers, "
      f"{n_kept_distress} of the {y_train.sum()} in distress")

The bank cannot drop an applicant: it must answer every test borrower, including the 3 321 with a missing value, and drop_model has never seen one. fillna fills them with the medians of the rows the model was fitted on:

X_test_filled = X_test.fillna(X_train[complete].median())
print(drop_model.classes_)
drop_model.predict_proba(X_test_filled)[:5].round(3)       # the first five test borrowers

🧠 Explain. What do the two numbers of each row add up to? Which column is the probability of distress?

✍️ Your answer:

…

We compare models with the ROC-AUC of the lecture. roc_auc_score(y_true, scores) computes it from the true labels and the scores, here on the rows drop_model was fitted on:

scores_fit = drop_model.predict_proba(X_train[complete])[:, 1]
roc_auc_score(y_train[complete], scores_fit)

✍️ Your turn. Complete report: auc_fit is the ROC-AUC on the rows the model was fitted on, auc_test the ROC-AUC on the test rows. Then run it on drop_model, whose test rows are X_test_filled.

fitted = {}          # every model of the lab, kept for the final table

def report(name, model, X_fit, y_fit, X_te):
    """ROC-AUC on the rows the model was fitted on, and on the test rows X_te."""
    fitted[name] = (model, X_fit, y_fit, X_te)
    # TODO: auc_fit = ...
    # TODO: auc_test = ...
    print(f"{name}: train ROC-AUC {auc_fit:.3f}   test ROC-AUC {auc_test:.3f}")

# TODO: report(...)

1.3 Filling the missing values

SimpleImputer(strategy="median") is a transformer, like the scaler of Lab 4: .fit learns the median of each column on the rows it is given, .transform fills each missing value with the median of its column.

imputer = SimpleImputer(strategy="median")
imputer.fit(X_train)
pd.Series(imputer.statistics_, index=X_train.columns)      # the medians learned on the training rows

X_train_filled = pd.DataFrame(imputer.transform(X_train), columns=X_train.columns, index=X_train.index)
print(f"missing incomes before: {X_train['MonthlyIncome'].isna().sum()}   "
      f"after: {X_train_filled['MonthlyIncome'].isna().sum()}")
X_train_filled[X_train["MonthlyIncome"].isna()].head(3)

In a pipeline it goes first, and the pipeline then takes the test rows as they are.

✍️ Your turn. Build impute_model: the pipeline of drop_model with a median SimpleImputer first, in a step named “impute”. Fit it on all the training rows.

# TODO: impute_model = ...
# TODO: impute_model.fit(...)

🔮 Predict. How do you expect impute_model to compare with drop_model on the test rows? Why?

✍️ Your answer:

…

report("median imputation", impute_model, X_train, y_train, X_test)

Part 2 — Scoring a classifier with scikit-learn

The lecture defined the scores; here is how scikit-learn computes them, for impute_model on the test rows.

2.1 Thresholds, precision, recall and the two curves

confusion_matrix(y_true, y_pred) counts the four cases: the rows are the truth (0, then 1), the columns the model’s answer (0, then 1). ConfusionMatrixDisplay draws it. precision_score, recall_score and f1_score take the same two arguments.

✍️ Your turn. Compute the confusion matrix cm of impute_model on the test rows, then its precision, recall and F1.

pred = impute_model.predict(X_test)
# TODO: cm = ...
# TODO: precision = ...
# TODO: recall = ...
# TODO: f1 = ...

ConfusionMatrixDisplay(cm, display_labels=["no distress", "distress"]).plot(colorbar=False)
plt.show()
print(f"precision {precision:.3f}   recall {recall:.3f}   F1 {f1:.3f}")

With the probabilities proba, proba >= t flags the borrowers whose probability is at least t, for any threshold t:

proba = impute_model.predict_proba(X_test)[:, 1]
for t in [0.5, 0.3, 0.2, 0.1, 0.05]:
    flagged = (proba >= t).astype(int)
    print(f"t = {t:4.2f}: {flagged.sum():5d} flagged   precision {precision_score(y_test, flagged):.3f}   "
          f"recall {recall_score(y_test, flagged):.3f}   F1 {f1_score(y_test, flagged):.3f}")

🧠 Explain. Which threshold would you pick for the bank, and what does it cost?

✍️ Your answer:

…

roc_curve(y_true, proba) returns the false positive rates, the true positive rates and the thresholds of the ROC curve, and roc_auc_score its area. precision_recall_curve(y_true, proba) returns the precisions, the recalls and the thresholds of the PR curve, and average_precision_score(y_true, proba) its area, the PR-AUC (scikit-learn calls it the average precision).

fpr, tpr, _ = roc_curve(y_test, proba)
prec, rec, _ = precision_recall_curve(y_test, proba)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.8))
ax1.plot(fpr, tpr, color=BLUE, label=f"ROC-AUC {roc_auc_score(y_test, proba):.3f}")
ax1.plot([0, 1], [0, 1], color=GRAY, ls="--", label="random ranking")
ax1.set_xlabel("false positive rate"); ax1.set_ylabel("true positive rate (recall)"); ax1.legend()
ax2.plot(rec, prec, color=BLUE, label=f"PR-AUC {average_precision_score(y_test, proba):.3f}")
ax2.axhline(y_test.mean(), color=GRAY, ls="--", label="random ranking")
ax2.set_xlabel("recall"); ax2.set_ylabel("precision"); ax2.legend()
plt.show()

2.2 Tuning the penalty

As in Lab 4, GridSearchCV chooses a hyperparameter by 5-fold cross-validation on the training rows. The parameter is C of the step "model", so its key is "model__C"; the score to maximize is the ROC-AUC, scoring="roc_auc"; and return_train_score=True also keeps the score on the training folds.

🔮 Predict. impute_model used C=1. Do you expect the best C to improve its test ROC-AUC by much? Why?

✍️ Your answer:

…

✍️ Your turn. Tune the C of impute_model over Cs by 5-fold cross-validation, scored by the ROC-AUC and keeping the training scores.

Cs = np.logspace(-2, 3, 11)
# TODO: tuned = ...
# TODO: tuned.fit(...)

print(f"best C {tuned.best_params_['model__C']:.3g}   CV ROC-AUC {tuned.best_score_:.3f}")
report("median imputation, C tuned", tuned, X_train, y_train, X_test)

plt.semilogx(Cs, tuned.cv_results_["mean_train_score"], "o-", color=BLUE, ms=4, label="training folds")
plt.semilogx(Cs, tuned.cv_results_["mean_test_score"], "o-", color=RED, ms=4, label="validation folds (CV)")
plt.xlabel("C"); plt.ylabel("ROC-AUC"); plt.legend()
plt.show()

🧠 Explain. What does this curve say about the penalty?

✍️ Your answer:

…

Part 3 — Looking at the data

3.1 Diagnostics

Before trying another model, look at what this one learned. The columns were standardized, so each coefficient is the change in the log-odds of distress for one standard deviation of its column.

pd.Series(impute_model["model"].coef_[0], index=X_train.columns).round(3)

🧠 Explain. BorrowedToLimitRatio measures how much of their credit limits a borrower uses. What sign would you expect for its coefficient, and what does the model give?

✍️ Your answer:

…

Now look at the raw values of this ratio. s.plot.box(vert=False) draws a boxplot of the Series s: the box spans the middle half of the values, and the circles are the values far from it. s.quantile(q) returns the value below which a share q of s lies: 99 % of the values are below s.quantile(0.99).

ratio = "BorrowedToLimitRatio"
X_train[ratio].plot.box(vert=False, figsize=(8, 2))
plt.show()
X_train[ratio].quantile([0.5, 0.9, 0.99, 0.995, 0.999, 1]).round(3)

🧠 Explain. Which values of the ratio cannot be real?

✍️ Your answer:

…

The ratio is the amount a borrower owes on their credit cards and lines of credit, divided by their credit limits. It usually lies between 0, nothing used, and 1, the limits reached: 97.9 % of the training borrowers are in that range. A borrower near or above 1 has used up their credit, so a larger ratio, above 1 especially, should mean more distress. The number of borrowers and the share in distress by slice of the ratio:

edges = [0, 0.25, 0.5, 1, 2, 5, 15, np.inf]
labels = [f"{a:g}-{b:g}" for a, b in zip(edges[:-2], edges[1:-1])] + ["above 15"]
slices = pd.cut(X_train[ratio], edges, labels=labels, include_lowest=True)
by_slice = y_train.groupby(slices, observed=False).agg(["size", "mean"])

fig, axes = plt.subplots(2, 1, figsize=(6.5, 5), sharex=True)
by_slice["size"].plot.bar(ax=axes[0], color=GRAY, logy=True)
axes[0].set_ylim(1, None); axes[0].set_ylabel("borrowers")
by_slice["mean"].plot.bar(ax=axes[1], color=RED, rot=0, label="share in distress")
axes[1].axhline(y_train.mean(), color=GRAY, ls="--", lw=1, label="all training borrowers")
axes[1].set_xlabel(ratio); axes[1].set_ylabel("share in distress"); axes[1].legend()
plt.show()

The share in distress climbs with the ratio, from 2.1 % below 0.25 to 41.7 % between 1 and 2: a strong signal, which the coefficient of −0.016 misses. No borrower has a ratio between 5 and 15, and 53 are above 15.

A borrower can owe a little more than their credit limits, but not hundreds of times more.

✍️ Your turn. Print the 60 largest ratios of the training rows and use them to choose the cutoff above which the ratio cannot be real. First time here: np.sort(a)[::-1] sorts a from the largest to the smallest value. Then build ratio_bad, True for the training borrowers above the cutoff, and compute the share in distress among them, rate_bad, and among the training borrowers just over their limits, with a ratio between 1 and 2, rate_over.

# TODO: print the 60 largest ratios of the training rows
# TODO: cutoff = the value above which the ratio cannot be real
# TODO: ratio_bad = ...
over_limit = X_train[ratio].between(1, 2, inclusive="right")
# TODO: rate_bad = ...
# TODO: rate_over = ...
print(f"above the cutoff: {ratio_bad.sum()} training borrowers, {rate_bad:.1%} in distress")
print(f"between 1 and 2:  {over_limit.sum()} training borrowers, {rate_over:.1%} in distress")
print(f"all:              {len(y_train)} training borrowers, {y_train.mean():.1%} in distress")
print(f"test borrowers above the cutoff: {(X_test[ratio] > cutoff).sum()}")

🧠 Explain. Why remove the borrowers above the cutoff, but keep those just over their limits?

✍️ Your answer:

…

🧠 Explain. We drop the training borrowers above the cutoff. What should happen to a test borrower above it?

✍️ Your answer:

…

s.mask(cond) replaces with NaN every value of s where cond is True, so the imputer of the pipeline fills it like any missing value.

🔮 Predict. Once these borrowers are gone, what sign do you expect for the coefficient of the ratio? Why?

✍️ Your answer:

…

✍️ Your turn. Build the training rows below the cutoff and the test rows with their ratios above the cutoff masked. Then fit clean_model, a new pipeline with the same three steps as impute_model, on those training rows, and report it.

# TODO: X_train_clean, y_train_clean = ...
X_test_clean = X_test.copy()
# TODO: X_test_clean[ratio] = ...
# TODO: clean_model = ...
# TODO: clean_model.fit(...)

report("impossible ratios removed", clean_model, X_train_clean, y_train_clean, X_test_clean)

for name, m in [("impute_model", impute_model), ("clean_model", clean_model)]:
    print(f"{name}: standard deviation of the ratio in the scaler {m['scale'].scale_[0]:.3f}   "
          f"its coefficient {m['model'].coef_[0][0]:.3f}")
pd.Series(clean_model["model"].coef_[0], index=X_train.columns).round(3)

Part 4 — The missing incomes, and the final table

4.1 Why is a value missing?

A value can be missing for three kinds of reasons. Let’s consider the variable Income as an example with missing values:

MCAR (missing completely at random): the missingness is independent of everything; (for e.g: random system error at data entry)
MAR (missing at random): the missingness can be explained by some observed columns; (for e.g by how data was collected: online vs at the bank in person)
MNAR (missing not at random): the missingness depends on the missing value itself or on the outcome; (for example, borrowers with a low income are less likely to share their income).

✍️ Your turn. Compute the share of borrowers in distress among the clean training borrowers whose income is missing, rate_missing, and among the others, rate_observed.

income_missing = X_train_clean["MonthlyIncome"].isna()
# TODO: rate_missing = ...
# TODO: rate_observed = ...
print(f"in distress: {rate_missing:.1%} when the income is missing, {rate_observed:.1%} when it is observed")

🧠 Explain. Which of the three mechanisms do these two shares rule out?

✍️ Your answer:

…

The median imputer turns a missing income into the median income. SimpleImputer(strategy="median", add_indicator=True) also adds, for each column with missing values in the training rows (here the income and the debt-to-income ratio), a 0/1 column equal to 1 where the value was missing. The model then gets its own coefficient for “the income is missing”. Under MAR, that column should add little once the model sees the other columns.

✍️ Your turn. Fit indicator_model: the pipeline of clean_model with add_indicator=True in its imputer, on the same clean rows, and report it.

# TODO: indicator_model = ...
# TODO: indicator_model.fit(...)

report("+ missing indicator", indicator_model, X_train_clean, y_train_clean, X_test_clean)

Could this gain come from which borrowers happen to be in the test set? Another 10 000 test borrowers would give other ROC-AUCs. We have only one test set, but the bootstrap imitates drawing new ones: draw 10 000 borrowers with replacement from the test set (some come twice, others not at all), score both models on them, and repeat 1 000 times. The middle 95 % of the 1 000 scores, between their 2.5th and 97.5th percentiles, is a 95 % confidence interval for the test ROC-AUC.

✍️ Your turn. Complete the bootstrap: i, the positions of the drawn borrowers, then the interval lo, hi of each model. First time here: rng.integers(0, n, n) draws n positions between 0 and n − 1 with replacement, and np.percentile(a, [2.5, 97.5]) returns the two percentiles of a.

p_clean = clean_model.predict_proba(X_test_clean)[:, 1]
p_indicator = indicator_model.predict_proba(X_test_clean)[:, 1]
y_t = y_test.to_numpy()
n = len(y_t)

rng = np.random.default_rng(0)
auc_clean, auc_indicator = [], []
for _ in range(1000):
    # TODO: i to be defined as n positions drawn with replacement
    auc_clean.append(roc_auc_score(y_t[i], p_clean[i]))
    auc_indicator.append(roc_auc_score(y_t[i], p_indicator[i]))
auc_clean, auc_indicator = np.array(auc_clean), np.array(auc_indicator)

for name, p, aucs in [("clean_model", p_clean, auc_clean), ("indicator_model", p_indicator, auc_indicator)]:
    # TODO: lo, hi to be defined as the 2.5th and 97.5th percentiles of aucs
    print(f"{name}: test ROC-AUC {roc_auc_score(y_t, p):.3f}, 95 % interval [{lo:.3f}, {hi:.3f}]")

Both models were scored on the same drawn borrowers, so each draw also gives a gain, auc_indicator - auc_clean: a paired bootstrap.

gains = auc_indicator - auc_clean
lo, hi = np.percentile(gains, [2.5, 97.5])
print(f"gain in test ROC-AUC {roc_auc_score(y_t, p_indicator) - roc_auc_score(y_t, p_clean):+.3f}, "
      f"95 % interval [{lo:+.3f}, {hi:+.3f}]")

🧠 Explain. Could the gain of the indicator be the luck of the test draw? Why?

✍️ Your answer:

…

🧠 Explain. Which mechanism does this file show? Why does the indicator help here?

✍️ Your answer:

…

4.2 The models side by side

The PR-AUC depends on the share of positives in the rows scored: the training PR-AUC of drop_model is computed on its kept rows, with 2.2 % of positives, so it cannot be compared with its test PR-AUC.

✍️ Your turn. Complete the table of the lab: for each model, pr holds its PR-AUC on its training rows and on the test rows.

table = {}
for name, (model, X_fit, y_fit, X_te) in fitted.items():
    p_fit, p_te = model.predict_proba(X_fit)[:, 1], model.predict_proba(X_te)[:, 1]
    roc = (roc_auc_score(y_fit, p_fit), roc_auc_score(y_test, p_te))
    # TODO: pr = (train PR-AUC, test PR-AUC)
    table[name] = [*roc, *pr]

pd.DataFrame(table, index=["train ROC-AUC", "test ROC-AUC", "train PR-AUC", "test PR-AUC"]).T.round(3)

🧠 Explain. What lessons do you draw from comparing these models?

✍️ Your answer:

…