Logistic Regression: Full Implementation

You've worked through the sigmoid, the cross-entropy loss, the gradient update, and the evaluation metrics. Now it's time to run logistic regression on a real dataset — Breast Cancer Wisconsin — from start to finish. This is the moment where every formula you've traced by hand gets exercised on 569 patient records with 30 features each.

The goal is not just to call LogisticRegression().fit() and get a score. It's to understand what each component does: why stratify=y matters for the split, how StandardScaler changes coefficient interpretation, what the confusion matrix tells you about clinical risk, and how regularization strength (C) affects cross-validated performance. Every code block produces a number you can verify against your own run.

What an End-to-End Implementation Covers

A complete logistic regression pipeline has six stages: load and explore the data (check class balance, feature ranges), split with stratification to preserve the class ratio, scale features to zero mean and unit variance, fit the model with appropriate solver and regularization, evaluate using multiple metrics (not just accuracy), and tune the key hyperparameter (C). The scaler and model must be saved together — you can't deploy one without the other.

The Plan — Six Steps from Data to Deployment

We'll load and inspect the Breast Cancer dataset, split and scale it properly, fit a logistic regression with default settings, evaluate the predictions in clinical context, tune the regularization strength with cross-validation, and finally save the model and scaler for deployment.

Anchor dataset: Breast Cancer Wisconsin — predict malignant (0) vs benign (1) from 30 tumor features.

from sklearn.datasets import load_breast_cancer
import pandas as pd
import numpy as np

data = load_breast_cancer()
df = pd.DataFrame(data.data, columns=data.feature_names)
df['target'] = data.target  # 0=malignant, 1=benign

print(f"Shape: {df.shape}")
print(df['target'].value_counts())
print(df[['mean radius', 'mean texture', 'mean area']].describe().round(2))

Shape: (569, 31)
1    357
0    212
Name: target, dtype: int64

       mean radius  mean texture  mean area
count       569.00        569.00     569.00
mean         14.13         19.29     654.89
std           3.52          4.30     351.91
min           6.98          9.71     143.50
max          28.11         39.28    2501.00

Class distribution: 212 malignant (37.3%) and 357 benign (62.7%) — mild imbalance, but not severe enough to require special treatment for this dataset. Feature ranges span orders of magnitude (radius ~7–28, area ~143–2501), making StandardScaler essential.

Step 1: Train/Test Split and Preprocessing

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X = data.data
y = data.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)
print(f"Train: {X_train.shape}, Test: {X_test.shape}")
print(f"Train class ratio: {y_train.mean():.3f}")
print(f"Test class ratio:  {y_test.mean():.3f}")

Train: (455, 30), Test: (114, 30)
Train class ratio: 0.627
Test class ratio:  0.632

stratify=y ensures both splits preserve the 63/37 class ratio. Without stratify, a random split could produce a test set with a different ratio — inflating or deflating all class-conditional metrics.

scaler = StandardScaler()
X_train_sc = scaler.fit_transform(X_train)
X_test_sc  = scaler.transform(X_test)

The scaler is fit on training data only — test set is transformed using training mean and std.

✓ Step 1 complete. Train set: 455 samples (63% benign). Test set: 114 samples (63% benign) — stratification preserved the exact class ratio. All features are now zero-mean, unit-variance.

Step 2: Fitting the Model — Coefficients and Feature Importance

from sklearn.linear_model import LogisticRegression

model = LogisticRegression(C=1.0, solver='lbfgs', max_iter=10000, random_state=42)
model.fit(X_train_sc, y_train)

print(f"Intercept: {model.intercept_[0]:.4f}")
coef_series = pd.Series(np.abs(model.coef_[0]), index=data.feature_names)
print("\nTop 5 coefficients by magnitude:")
print(coef_series.nlargest(5).round(4))

Intercept: 0.3524

Top 5 coefficients by magnitude:
worst concave points     2.8731
worst perimeter          2.3142
mean concave points      1.9823
worst radius             1.7156
mean area                1.5601

The largest absolute coefficients identify the most influential features. worst concave points = 2.87 means a 1 standard deviation increase in this feature changes the log-odds of benign by 2.87 — the strongest signal in the model.

✓ Step 2 complete. Model fitted with intercept = 0.3524. Top features by coefficient magnitude: worst concave points (2.87), worst perimeter (2.31), mean concave points (1.98). Each represents a 1‑std increase in log-odds of benign.

Step 3: Predictions and Probability Outputs

y_pred = model.predict(X_test_sc)
y_prob = model.predict_proba(X_test_sc)[:, 1]  # P(benign)

for i in range(5):
    print(f"Sample {i}: P(benign)={y_prob[i]:.4f}, predicted={y_pred[i]}, true={y_test[i]}")

Sample 0: P(benign)=0.9823, predicted=1, true=1
Sample 1: P(benign)=0.0041, predicted=0, true=0
Sample 2: P(benign)=0.8916, predicted=1, true=1
Sample 3: P(benign)=0.9991, predicted=1, true=1
Sample 4: P(benign)=0.0018, predicted=0, true=0

The model is confident on all 5 shown samples — probabilities near 0 or 1. predict_proba()[:,1] extracts the probability of the positive class (benign = 1).

✓ Step 3 complete. The model assigns extreme probabilities — near 0 or 1 — to the first 5 test samples, indicating high confidence in its predictions.

Step 4: Full Evaluation — Confusion Matrix and AUC

from sklearn.metrics import (confusion_matrix, classification_report,
                              roc_auc_score, average_precision_score)

cm = confusion_matrix(y_test, y_pred)
print("Confusion Matrix:")
print(cm)
print()
print(classification_report(y_test, y_pred, target_names=['Malignant', 'Benign']))
print(f"AUC-ROC: {roc_auc_score(y_test, y_prob):.4f}")
print(f"AUC-PR:  {average_precision_score(y_test, y_prob):.4f}")

Confusion Matrix:
[[40  2]
 [ 1 71]]

              precision    recall  f1-score   support

   Malignant       0.98      0.95      0.96        42
      Benign       0.97      0.99      0.98        72

    accuracy                           0.97       114

AUC-ROC: 0.9981
AUC-PR:  0.9989

Interpreting the confusion matrix:

• TP = 71: benign tumors correctly identified as benign
• TN = 40: malignant tumors correctly identified as malignant
• FP = 2: benign tumors misclassified as malignant (unnecessary biopsy)
• FN = 1: malignant tumor called benign — the most serious error

FN = 1 is the critical failure: one cancer patient would receive a false all-clear. Medical screening systems tune the threshold toward higher recall (lower threshold), accepting more FP to ensure near-zero FN.

✓ Step 4 complete. Confusion matrix: TP=71, TN=40, FP=2, FN=1. AUC-ROC = 0.9981, AUC-PR = 0.9989. The single FN is the critical error — one malignant tumor called benign.

Step 5: Hyperparameter Tuning — Regularization Strength C

C = 1/λ in sklearn. Smaller C = stronger regularization. Note the inverse convention from Ridge/Lasso:

from sklearn.model_selection import cross_val_score

C_values = [0.001, 0.01, 0.1, 1, 10, 100]
for C in C_values:
    m = LogisticRegression(C=C, solver='lbfgs', max_iter=10000)
    scores = cross_val_score(m, X_train_sc, y_train, cv=5, scoring='roc_auc')
    print(f"C={C:>7}: CV AUC = {scores.mean():.4f} ± {scores.std():.4f}")

C=  0.001: CV AUC = 0.9841 ± 0.0082
C=   0.01: CV AUC = 0.9944 ± 0.0049
C=    0.1: CV AUC = 0.9970 ± 0.0030
C=      1: CV AUC = 0.9975 ± 0.0031
C=     10: CV AUC = 0.9974 ± 0.0030
C=    100: CV AUC = 0.9973 ± 0.0030

At C=0.001 (strong regularization): AUC drops to 0.984 — the model is over-constrained. At C=1–100: performance plateaus — the data is well-conditioned and very little additional regularization helps. Optimal C≈1 (sklearn default).

✓ Step 5 complete. Cross-validated AUC peaks at C=1 (0.9975) and plateaus across C=1–100. Strong regularization (C=0.001) degrades AUC to 0.984.

Step 6: Solver Comparison

Different solvers implement the same logistic regression objective but differ in numerical method and supported penalties:

for solver in ['lbfgs', 'liblinear', 'saga']:
    m = LogisticRegression(solver=solver, max_iter=10000, random_state=42)
    m.fit(X_train_sc, y_train)
    acc = m.score(X_test_sc, y_test)
    print(f"solver={solver:12s}: accuracy={acc:.4f}")

solver=lbfgs       : accuracy=0.9737
solver=liblinear   : accuracy=0.9737
solver=saga        : accuracy=0.9737

All three converge to the same solution on this well-conditioned dataset. Choose based on dataset size: lbfgs for small-medium (default), saga for large datasets or ElasticNet penalty, liblinear for L1 penalty on small datasets.

✓ Step 6 complete. All three solvers (lbfgs, liblinear, saga) converge to 97.37% accuracy. The choice matters for dataset size and supported penalties, not for final accuracy on this well-conditioned data.

Step 7: Saving and Loading for Deployment

import joblib

# Save both model and scaler — they must travel together
joblib.dump({'model': model, 'scaler': scaler}, 'breast_cancer_lr.pkl')

# Load and infer on a new sample
bundle = joblib.load('breast_cancer_lr.pkl')
new_sample = X_test[:1]                                        # use one test sample for demo
new_scaled = bundle['scaler'].transform(new_sample)
prob = bundle['model'].predict_proba(new_scaled)[0, 1]
print(f"P(benign) = {prob:.4f}")

P(benign) = 0.9823

Always save the scaler alongside the model. If you load the model and apply a new or unfitted scaler, predictions will be wrong — the model was trained on features with specific mean and std, and those must be reproduced at inference time.

✓ Step 7 complete. Model and scaler saved together as breast_cancer_lr.pkl. Inference on a new sample gives P(benign) = 0.9823 — matching the original test prediction.

Pipeline walkthrough

Related Concepts

This implementation applies everything from Posts 01–04: the sigmoid, cross-entropy loss, gradient descent, precision/recall, and AUC evaluation on a real dataset. The pipeline pattern — split, scale, fit, evaluate, tune, deploy — carries to every sklearn model. The hyperparameter tuning here (manual C sweep) is a preview of GridSearchCV and RandomizedSearchCV covered in Post 06.

Honest Limitations

AUC-ROC = 0.9981 and accuracy = 97% on this dataset are exceptionally high — the 30 features include many strong predictors of cancer (concave points, perimeter, area). Real clinical datasets are harder. This implementation also skips one important step for clinical use: calibration. The model's predict_proba output of 0.98 should genuinely mean 98% of samples with that score are benign — sklearn.calibration.CalibratedClassifierCV checks and corrects this.

The FN=1 error highlights the gap between ML metrics and deployment risk. In a 10,000-sample clinical screen, 1% FNR means 100 missed cancers. The decision threshold should be set by the acceptable FNR for the use case — not by whatever maximizes accuracy or F1.

Test Your Understanding

1. The model uses C=1.0. The AUC at C=0.001 is 0.9841 and at C=1 is 0.9975. The difference is 0.013 AUC points. Is this clinically meaningful for a cancer screening model? What metric would you focus on instead?

2. stratify=y was used in the split. Suppose the original dataset has 212 malignant and 357 benign. After an 80/20 split with stratify, exactly how many malignant and benign samples are in the test set?

3. FN=1 in the test set (1 malignant predicted as benign). If you lower the threshold from 0.5 to 0.2, what happens to FN and FP? Which confusion matrix entries could increase and which could only decrease?

4. The saga solver and lbfgs solver both achieve 97.37% accuracy but use different optimization algorithms. If training data doubled to 1,000 samples, which solver would you prefer and why?

5. A colleague saves only the model (without the scaler) to model.pkl. At inference time, they apply the raw features directly. Qualitatively, what happens to the model's predictions?