Counterfactual Explanations in XAI: Python Guide

The Story That Explains Counterfactual Explanations

A counterfactual explanation answers the question: "What is the smallest change to the input that would flip the model's output?" It is the XAI technique closest to how humans naturally think about decisions — not "why did this happen?" but "what would I need to do differently?"

Where Counterfactuals Fit in the XAI Landscape

Explainable AI (XAI) techniques are broadly divided by scope (global vs local) and approach (feature attribution vs example-based vs rule-based). Counterfactual explanations are local, example-based methods — they explain a single prediction by generating a contrast example.

The Four Properties of a Good Counterfactual

Not all counterfactuals are useful. A counterfactual that says "if you were 20 years younger and had a PhD, you'd be approved" is technically valid but practically useless. Good counterfactuals satisfy four properties.

The Mathematical Formulation

Counterfactual search is formulated as a constrained optimisation problem. The original formulation by Wachter et al. (2017) is:

Animated Diagram — How Counterfactual Search Works

The animation below shows a binary classifier's decision boundary in 2D feature space. Watch how the search algorithm nudges the original point (red) across the boundary to find the nearest valid counterfactual (green).

Major Counterfactual Algorithms Compared

The field has produced several distinct algorithms, each with different trade-offs between plausibility, diversity, and computational cost.

Animated Diagram — DICE Diversity vs Single CF

This animation contrasts finding a single counterfactual (Wachter-style) versus a diverse set (DICE-style). Each path represents a different recourse strategy.

Python Implementation — From Scratch (Wachter Method)

Let us build a minimal counterfactual explainer from scratch to understand what libraries like DICE do under the hood. We will use a loan approval dataset and a trained Random Forest classifier.

# ── Step 1: Build dataset and train model ───────────────────
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

np.random.seed(42)
n = 2000

# Synthetic loan data
df = pd.DataFrame({
    'income':       np.random.normal(35000, 12000, n).clip(10000, 150000),
    'credit_score': np.random.normal(640, 80, n).clip(300, 850),
    'debt_ratio':   np.random.uniform(0.1, 0.9, n),
    'employment_yr': np.random.exponential(5, n).clip(0, 30),
    'loan_amount':  np.random.normal(18000, 8000, n).clip(1000, 60000),
})

# Label: approved if strong financials
score = (
    (df['income'] / 50000) * 0.35 +
    (df['credit_score'] / 850) * 0.30 +
    (1 - df['debt_ratio']) * 0.20 +
    (df['employment_yr'] / 30) * 0.15
)
df['approved'] = (score + np.random.normal(0, 0.08, n) > 0.48).astype(int)

X = df.drop('approved', axis=1)
y = df['approved']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Train Random Forest (no scaling needed, but we show pipeline for completeness)
model = RandomForestClassifier(n_estimators=200, max_depth=8, random_state=42)
model.fit(X_train, y_train)
print(f"Test accuracy: {model.score(X_test, y_test):.3f}")

# ── Step 2: Wachter-style Counterfactual from scratch ───────
from scipy.optimize import minimize
from scipy.stats import median_abs_deviation

def wachter_cf(model, x_orig, y_desired, feature_names,
               immutable=None, lam=0.5, max_iter=500):
    """
    Wachter-style counterfactual finder.
    model       : sklearn classifier with predict_proba
    x_orig      : pd.Series — the original input row
    y_desired   : int — desired class (e.g. 1 = approved)
    immutable   : list of feature names that cannot change
    lam         : regularisation weight (higher = prioritise flipping)
    """
    x0 = x_orig.values.copy().astype(float)
    names = feature_names
    immut_idx = [names.index(f) for f in (immutable or [])]

    # MAD normalisation factors (from training data stats)
    mad = X_train.apply(lambda col: median_abs_deviation(col)).values + 1e-8

    def objective(x_cf):
        # Prediction loss
        proba = model.predict_proba(pd.DataFrame([x_cf], columns=names))[:, y_desired][0]
        pred_loss = (proba - 1.0) ** 2           # want proba → 1

        # L1 proximity (MAD-normalised)
        proximity = np.sum(np.abs(x_cf - x0) / mad)

        return lam * pred_loss + proximity

    def grad_constraint(x_cf):
        # Ensure immutable features don't change
        return x_cf[immut_idx] - x0[immut_idx]

    constraints = [{'type': 'eq', 'fun': grad_constraint}] if immut_idx else []

    result = minimize(
        objective,
        x0=x0,
        method='SLSQP',
        constraints=constraints,
        options={'maxiter': max_iter, 'ftol': 1e-6}
    )

    x_cf = pd.Series(result.x, index=names)
    proba = model.predict_proba(pd.DataFrame([result.x], columns=names))[:, y_desired][0]
    valid = proba > 0.5

    return x_cf, valid, proba


# ── Apply to Maria's application ────────────────────────────
maria = pd.Series({
    'income': 26500,
    'credit_score': 580,
    'debt_ratio': 0.52,
    'employment_yr': 1.5,
    'loan_amount': 20000,
})

feature_names = list(X.columns)
cf, valid, proba = wachter_cf(model, maria, y_desired=1,
                                feature_names=feature_names,
                                immutable=[],       # all mutable
                                lam=1.0)

print(f"\nOriginal prediction: {model.predict_proba(maria.to_frame().T)[0,1]:.3f}")
print(f"CF prediction:       {proba:.3f}  (Valid: {valid})")
print("\nFeature Changes:")
print(f"{'Feature':20s} {'Original':>12s} {'Counterfactual':>16s} {'Delta':>12s}")
print("-"*65)
for feat in feature_names:
    delta = cf[feat] - maria[feat]
    if abs(delta) > 0.01 * abs(maria[feat]):
        print(f"{feat:20s} {maria[feat]:>12.1f} {cf[feat]:>16.1f} {delta:>+12.1f}")

DICE Library — Diverse Counterfactuals in Practice

The DiCE library (Microsoft Research) is the most widely used production-grade counterfactual library. It supports sklearn, TensorFlow, PyTorch, and any custom model.

# ── DICE: Diverse Counterfactual Explanations ───────────────
import dice_ml
from dice_ml import Dice

# 1. Wrap dataset
d = dice_ml.Data(
    dataframe=pd.concat([X_train, y_train.rename('approved')], axis=1),
    continuous_features=['income', 'credit_score', 'debt_ratio',
                         'employment_yr', 'loan_amount'],
    outcome_name='approved'
)

# 2. Wrap model
m = dice_ml.Model(model=model, backend='sklearn')

# 3. Create explainer
exp = Dice(d, m, method='random')      # method: 'random' | 'genetic' | 'kdtree'

# 4. Generate 4 diverse counterfactuals for Maria
query = maria.to_frame().T
cf_result = exp.generate_counterfactuals(
    query_instances=query,
    total_CFs=4,
    desired_class='opposite',             # flip the prediction
    permitted_range={
        'income':       [26500, 80000],    # can only increase
        'debt_ratio':  [0.1, 0.52],         # can only decrease debt
        'employment_yr': [1.5, 30],          # can only increase tenure
    },
    features_to_vary=['income', 'credit_score', 'debt_ratio', 'employment_yr']
    # loan_amount not varied — she already fixed the amount
)

cf_result.visualize_as_dataframe(show_only_changes=True)

Interactive — Actionability Constraints

Actionability constraints are what separate useful counterfactuals from technically valid but practically useless ones. The interactive diagram below shows how constraining mutable features changes the counterfactual space.

Real-World Applications

Animated Diagram — Proximity vs Plausibility Trade-off

A fundamental tension in counterfactual generation: the nearest point that crosses the decision boundary may be off the data manifold — a combination of feature values that never occurs in reality. The animation below shows this trade-off.

How to Evaluate Counterfactual Explanations

Generating counterfactuals is only half the challenge. Evaluating their quality requires a multi-dimensional scorecard. No single metric captures everything — you need all five.

# ── Evaluation helper ───────────────────────────────────────
from scipy.spatial.distance import cdist

def evaluate_cfs(x_orig, cfs_df, model, X_train, y_desired=1):
    """Compute all 5 CF quality metrics."""
    x0 = x_orig.values
    cf_vals = cfs_df.values
    mad = X_train.apply(lambda c: median_abs_deviation(c)).values + 1e-8

    # 1. Validity
    preds = model.predict(cf_vals)
    validity = np.mean(preds == y_desired)

    # 2. Proximity (L1 MAD-normalised)
    proximity = np.mean([np.sum(np.abs(cf - x0) / mad) for cf in cf_vals])

    # 3. Sparsity (number of changed features)
    sparsity = np.mean([np.sum(np.abs(cf - x0) > 0.01) for cf in cf_vals])

    # 4. Diversity (pairwise distances)
    if len(cf_vals) > 1:
        pw = cdist(cf_vals, cf_vals, metric='minkowski', p=1)
        diversity = pw[np.triu_indices_from(pw, k=1)].mean()
    else:
        diversity = 0.0

    # 5. Plausibility: distance to nearest training point
    dists = cdist([cf for cf in cf_vals],
                  X_train.values, metric='minkowski', p=1)
    plausibility = np.mean(dists.min(axis=1))

    return {
        'validity':    round(validity, 3),
        'proximity':   round(proximity, 3),
        'sparsity':    round(sparsity, 2),
        'diversity':   round(diversity, 3),
        'plausibility': round(plausibility, 3),
    }

Common Pitfalls and How to Avoid Them

Counterfactuals vs SHAP — When to Use Each

Complete XAI Pipeline with Counterfactuals

Golden Rules