Case Study 2 — Meridian Financial Credit Scoring

MLP on Tabular Data with Gradient-Based Feature Importance

Context

Meridian Financial, a consumer lending institution, processes approximately 2 million credit applications annually. Their existing credit scoring model is a logistic regression with handcrafted features — the industry standard, chosen for its interpretability and regulatory compliance under the Equal Credit Opportunity Act (ECOA) and the Fair Credit Reporting Act (FCRA).

The data science team wants to evaluate whether an MLP can improve predictive accuracy. But accuracy alone is insufficient: Meridian must explain every adverse action (loan denial) to the applicant, citing the specific factors that contributed to the decision. This creates a tension that pervades production ML in regulated industries: nonlinear models are more accurate, but their decisions are harder to explain.

This case study builds the MLP, evaluates its accuracy advantage, and uses gradient-based analysis to understand which features drive predictions — laying the groundwork for the fairness audit in Chapter 31.

Data and Features

Meridian's credit scoring features span three categories: credit history, financial capacity, and application context.

import numpy as np
from typing import Dict, List, Tuple


def generate_credit_data(
    n_samples: int = 50_000,
    seed: int = 42,
) -> Dict[str, np.ndarray]:
    """Generate synthetic credit scoring data.

    The data-generating process includes nonlinear effects and
    feature interactions that motivate the MLP over logistic regression.

    Features:
        0: credit_score (normalized FICO-like, 0-1)
        1: debt_to_income (DTI ratio, 0-1)
        2: years_of_credit_history (log-transformed)
        3: num_open_accounts (count, normalized)
        4: num_derogatory_marks (count, normalized)
        5: total_credit_utilization (ratio, 0-1)
        6: annual_income (log-transformed, normalized)
        7: employment_length_years (normalized)
        8: loan_amount_requested (log-transformed, normalized)
        9: loan_to_income_ratio (derived)
        10: num_recent_inquiries (count, normalized)
        11: months_since_last_delinquency (normalized, -1 if never)
        12: total_revolving_balance (log-transformed, normalized)
        13: mortgage_indicator (binary)
        14: public_record_indicator (binary)

    Target: 1 = default within 24 months, 0 = no default.
    Base default rate: ~12%.

    Args:
        n_samples: Number of credit applications.
        seed: Random seed.

    Returns:
        Dictionary with train/val/test splits.
    """
    rng = np.random.default_rng(seed)
    n_features = 15

    # Generate correlated features (credit variables are correlated)
    mean = np.zeros(n_features)
    cov = np.eye(n_features)
    # Credit score correlates negatively with derogatory marks
    cov[0, 4] = cov[4, 0] = -0.5
    # DTI correlates with credit utilization
    cov[1, 5] = cov[5, 1] = 0.4
    # Income correlates with credit score
    cov[0, 6] = cov[6, 0] = 0.3
    # Employment length correlates with credit history
    cov[2, 7] = cov[7, 2] = 0.5

    X = rng.multivariate_normal(mean, cov, size=n_samples)

    # Clip binary features
    X[:, 13] = (X[:, 13] > 0.5).astype(float)
    X[:, 14] = (X[:, 14] > 1.0).astype(float)

    # Default probability: nonlinear function of features
    logit = (
        -2.0                                   # Base rate intercept
        - 1.5 * X[:, 0]                        # Credit score (strong linear)
        + 1.2 * X[:, 1]                        # DTI (strong linear)
        + 0.8 * X[:, 4]                        # Derogatory marks
        + 0.6 * X[:, 1] * X[:, 5]             # DTI x utilization interaction
        - 0.4 * X[:, 0] * X[:, 2]             # Credit score x history interaction
        + 0.5 * np.maximum(X[:, 10] - 0.5, 0) # Threshold: many inquiries = risky
        - 0.3 * X[:, 7] * (1 - X[:, 14])      # Employment helps if no public record
        + 0.3 * X[:, 8]                        # Larger loans = more risk
        + 0.7 * X[:, 14]                       # Public records
    )

    prob = 1 / (1 + np.exp(-logit))
    y = rng.binomial(1, prob)

    # Split: 70/15/15
    n_train = int(0.7 * n_samples)
    n_val = int(0.15 * n_samples)

    return {
        "X_train": X[:n_train],
        "y_train": y[:n_train],
        "X_val": X[n_train:n_train + n_val],
        "y_val": y[n_train:n_train + n_val],
        "X_test": X[n_train + n_val:],
        "y_test": y[n_train + n_val:],
        "feature_names": [
            "credit_score", "debt_to_income", "credit_history_years",
            "num_open_accounts", "num_derogatory_marks",
            "credit_utilization", "annual_income", "employment_years",
            "loan_amount", "loan_to_income", "recent_inquiries",
            "months_since_delinquency", "revolving_balance",
            "has_mortgage", "has_public_record",
        ],
    }

The data-generating process includes two types of nonlinearity that logistic regression cannot capture: (1) feature interactions (DTI $\times$ utilization, credit score $\times$ history length), and (2) a threshold effect on recent inquiries (risk increases only after more than a certain number of inquiries).

Model Comparison

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (
    roc_auc_score, log_loss, average_precision_score,
    precision_recall_curve,
)


def evaluate_credit_models(
    data: Dict[str, np.ndarray],
) -> Dict[str, Dict[str, float]]:
    """Train and compare logistic regression vs. MLP on credit data.

    Args:
        data: Dictionary with train/val/test splits.

    Returns:
        Dictionary of model results.
    """
    # Standardize
    scaler = StandardScaler()
    X_train_s = scaler.fit_transform(data["X_train"])
    X_val_s = scaler.transform(data["X_val"])
    X_test_s = scaler.transform(data["X_test"])

    results = {}

    # --- Logistic Regression ---
    lr = LogisticRegression(C=1.0, max_iter=1000, random_state=42)
    lr.fit(X_train_s, data["y_train"])
    y_prob_lr = lr.predict_proba(X_test_s)[:, 1]

    results["logistic_regression"] = {
        "auc": roc_auc_score(data["y_test"], y_prob_lr),
        "logloss": log_loss(data["y_test"], y_prob_lr),
        "ap": average_precision_score(data["y_test"], y_prob_lr),
    }

    # --- MLP (numpy) ---
    config = MLPConfig(
        layer_dims=[15, 64, 32, 16, 1],
        activation="relu",
        output_activation="sigmoid",
        seed=42,
    )
    mlp = NumpyMLP(config)
    mlp.train(
        X_train=X_train_s,
        y_train=data["y_train"],
        X_val=X_val_s,
        y_val=data["y_val"],
        epochs=300,
        batch_size=256,
        lr=0.003,
        verbose=False,
    )
    y_prob_mlp = mlp.forward(X_test_s).ravel()

    results["mlp"] = {
        "auc": roc_auc_score(data["y_test"], y_prob_mlp),
        "logloss": binary_cross_entropy(y_prob_mlp, data["y_test"]),
        "ap": average_precision_score(data["y_test"], y_prob_mlp),
    }

    return results


data = generate_credit_data()
results = evaluate_credit_models(data)

print(f"{'Model':<25} {'AUC':>8} {'LogLoss':>10} {'AP':>8}")
print("-" * 55)
for model_name, metrics in results.items():
    print(f"{model_name:<25} {metrics['auc']:>8.4f} "
          f"{metrics['logloss']:>10.4f} {metrics['ap']:>8.4f}")

Expected results:

The MLP provides a 2.3-point AUC improvement. In credit scoring, where models process millions of decisions annually and each percentage point of AUC translates to significant revenue (reduced defaults) and reduced harm (fewer incorrect denials), this improvement is substantial.

Gradient-Based Feature Importance

Regulatory requirements demand that Meridian explain its credit decisions. For logistic regression, the coefficients directly indicate feature importance. For the MLP, we use gradient-based analysis: the gradient of the output with respect to each input feature indicates how much a small change in that feature would change the prediction.

def compute_gradient_importance(
    model: NumpyMLP,
    X: np.ndarray,
    feature_names: List[str],
) -> Dict[str, float]:
    """Compute feature importance via average absolute input gradients.

    For each test example, compute the gradient of the output with
    respect to the input features. The average absolute gradient
    across examples gives a measure of each feature's importance.

    Args:
        model: Trained NumpyMLP.
        X: Input data of shape (n, d).
        feature_names: Names of input features.

    Returns:
        Dictionary mapping feature names to importance scores.
    """
    n, d = X.shape
    importances = np.zeros(d)

    for i in range(n):
        x_i = X[i:i+1]  # (1, d)
        model.forward(x_i)

        # Backpropagate to get gradient w.r.t. input
        L = len(model.config.layer_dims) - 1
        y_pred = model.cache[f"A{L}"]

        # Start with output gradient (for sigmoid output, dL/dz = 1 for
        # the gradient of the output itself, not the loss)
        dA = np.ones_like(y_pred)

        # Propagate backward through all layers
        for l in range(L, 0, -1):
            if l == L:
                dZ = dA * sigmoid_gradient(model.cache[f"Z{l}"])
            else:
                dZ = dA * relu_gradient(model.cache[f"Z{l}"])

            if l > 1:
                dA = dZ @ model.params[f"W{l}"]
            else:
                # Gradient w.r.t. input
                input_grad = dZ @ model.params[f"W{l}"]  # (1, d)
                importances += np.abs(input_grad.ravel())

    importances /= n

    # Normalize to sum to 1
    importances /= importances.sum()

    return {
        name: float(imp)
        for name, imp in zip(feature_names, importances)
    }


# Compute and display importance
importance = compute_gradient_importance(
    mlp, X_test_s, data["feature_names"]
)

print(f"\n{'Feature':<30} {'Importance':>12}")
print("-" * 44)
for name, score in sorted(
    importance.items(), key=lambda x: x[1], reverse=True
):
    bar = "#" * int(score * 100)
    print(f"{name:<30} {score:>12.4f}  {bar}")

Expected top features (approximate):

The gradient importance ranking aligns with domain knowledge and with the true data-generating process. Credit score and DTI are the dominant predictors, followed by derogatory marks and utilization. The MLP has learned the feature interactions (DTI $\times$ utilization, credit score $\times$ history) even though these were not explicitly engineered as input features.

Individual Decision Explanation

For adverse action notices, Meridian needs per-applicant explanations — not just global feature importance.

def explain_decision(
    model: NumpyMLP,
    x: np.ndarray,
    feature_names: List[str],
    top_k: int = 5,
) -> List[Tuple[str, float, str]]:
    """Explain a single credit decision using input gradients.

    Args:
        model: Trained NumpyMLP.
        x: Single input of shape (d,).
        feature_names: Feature names.
        top_k: Number of top factors to return.

    Returns:
        List of (feature_name, gradient_value, direction) tuples.
    """
    x_input = x.reshape(1, -1)
    model.forward(x_input)

    L = len(model.config.layer_dims) - 1
    dA = np.ones((1, 1))

    for l in range(L, 0, -1):
        if l == L:
            dZ = dA * sigmoid_gradient(model.cache[f"Z{l}"])
        else:
            dZ = dA * relu_gradient(model.cache[f"Z{l}"])

        if l > 1:
            dA = dZ @ model.params[f"W{l}"]
        else:
            input_grad = (dZ @ model.params[f"W{l}"]).ravel()

    # Feature contribution = gradient * feature value
    contributions = input_grad * x

    # Sort by absolute contribution
    sorted_indices = np.argsort(np.abs(contributions))[::-1][:top_k]

    explanations = []
    for idx in sorted_indices:
        direction = "increases risk" if contributions[idx] > 0 else "decreases risk"
        explanations.append(
            (feature_names[idx], float(contributions[idx]), direction)
        )

    return explanations


# Example: explain a high-risk applicant
high_risk_idx = np.argmax(y_prob_mlp)
prediction = y_prob_mlp[high_risk_idx]
print(f"\nApplicant #{high_risk_idx}: P(default) = {prediction:.4f}")
print(f"{'Factor':<30} {'Contribution':>14} {'Direction':<20}")
print("-" * 66)

explanations = explain_decision(
    mlp, X_test_s[high_risk_idx], data["feature_names"]
)
for name, contrib, direction in explanations:
    print(f"{name:<30} {contrib:>14.4f} {direction:<20}")

Limitations and Forward References

This gradient-based explanation is a first approximation. It has known limitations:

1. Gradients reflect local sensitivity, not global importance. The gradient at a specific input tells you the effect of an infinitesimal perturbation, not the effect of a large change. Integrated gradients (Exercise 6.25) address this by averaging gradients along the path from a baseline to the input.

2. The explanation depends on the model, not the truth. If the model has learned a spurious correlation, the gradient faithfully reports that correlation as important. The fairness audit in Chapter 31 will examine whether the model's use of proxy features creates disparate impact across demographic groups.

3. Gradient explanations are not contrastive. They do not answer "what would need to change for this applicant to be approved?" — a question that requires counterfactual analysis (Chapter 19).

4. Regulatory compliance requires more than gradients. ECOA mandates specific adverse action reason codes. Mapping continuous gradient values to discrete reason codes is an additional engineering challenge covered in Chapter 35.

Despite these limitations, gradient-based analysis provides a principled starting point for understanding neural network decisions on tabular data — a bridge between the black-box MLP and the fully interpretable logistic regression it replaces.