Defining a custom risk for binary classification¶

MAPIE ships with several predefined risks and performance metrics for binary classification (precision, recall, accuracy, false positive rate...). When the metric you care about is not one of them, you can define your own with BinaryRisk and control it exactly like a predefined one.

In this example, we show how to define specificity (the true negative rate) as a custom risk and use it with the BinaryClassificationController.

Specificity is the proportion of actual negatives that are correctly predicted as negative. It is a natural target in screening problems where we want to avoid raising too many alarms on healthy/benign cases, while keeping recall (the ability to catch the positive cases) as high as possible.

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import make_circles
from sklearn.metrics import recall_score
from sklearn.neural_network import MLPClassifier

from mapie.risk_control import BinaryClassificationController, BinaryRisk
from mapie.utils import train_conformalize_test_split

RANDOM_STATE = 42

First, load the dataset and split it into training, calibration (for conformalization), and test sets, then fit a classifier on the training data.

X, y = make_circles(n_samples=5000, noise=0.3, factor=0.3, random_state=RANDOM_STATE)
(X_train, X_calib, X_test, y_train, y_calib, y_test) = train_conformalize_test_split(
    X,
    y,
    train_size=0.7,
    conformalize_size=0.1,
    test_size=0.2,
    random_state=RANDOM_STATE,
)

clf = MLPClassifier(max_iter=150, random_state=RANDOM_STATE)
clf.fit(X_train, y_train)

MLPClassifier(max_iter=150, random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

MLPClassifier

?Documentation for MLPClassifieriFitted

Parameters

	max_iter max_iter: int, default=200 Maximum number of iterations. The solver iterates until convergence (determined by 'tol') or this number of iterations. For stochastic solvers ('sgd', 'adam'), note that this determines the number of epochs (how many times each data point will be used), not the number of gradient steps.	150
	random_state random_state: int, RandomState instance, default=None Determines random number generation for weights and bias initialization, train-test split if early stopping is used, and batch sampling when solver='sgd' or 'adam'. Pass an int for reproducible results across multiple function calls. See :term:`Glossary <random_state>`.	42
	hidden_layer_sizes hidden_layer_sizes: array-like of shape(n_layers - 2,), default=(100,) The ith element represents the number of neurons in the ith hidden layer.	(100,)
	activation activation: {'identity', 'logistic', 'tanh', 'relu'}, default='relu' Activation function for the hidden layer. - 'identity', no-op activation, useful to implement linear bottleneck, returns f(x) = x - 'logistic', the logistic sigmoid function, returns f(x) = 1 / (1 + exp(-x)). - 'tanh', the hyperbolic tan function, returns f(x) = tanh(x). - 'relu', the rectified linear unit function, returns f(x) = max(0, x)	'relu'
	solver solver: {'lbfgs', 'sgd', 'adam'}, default='adam' The solver for weight optimization. - 'lbfgs' is an optimizer in the family of quasi-Newton methods. - 'sgd' refers to stochastic gradient descent. - 'adam' refers to a stochastic gradient-based optimizer proposed by Kingma, Diederik, and Jimmy Ba For a comparison between Adam optimizer and SGD, see :ref:`sphx_glr_auto_examples_neural_networks_plot_mlp_training_curves.py`. Note: The default solver 'adam' works pretty well on relatively large datasets (with thousands of training samples or more) in terms of both training time and validation score. For small datasets, however, 'lbfgs' can converge faster and perform better.	'adam'
	alpha alpha: float, default=0.0001 Strength of the L2 regularization term. The L2 regularization term is divided by the sample size when added to the loss. For an example usage and visualization of varying regularization, see :ref:`sphx_glr_auto_examples_neural_networks_plot_mlp_alpha.py`.	0.0001
	batch_size batch_size: int, default='auto' Size of minibatches for stochastic optimizers. If the solver is 'lbfgs', the classifier will not use minibatch. When set to "auto", `batch_size=min(200, n_samples)`.	'auto'
	learning_rate learning_rate: {'constant', 'invscaling', 'adaptive'}, default='constant' Learning rate schedule for weight updates. - 'constant' is a constant learning rate given by 'learning_rate_init'. - 'invscaling' gradually decreases the learning rate at each time step 't' using an inverse scaling exponent of 'power_t'. effective_learning_rate = learning_rate_init / pow(t, power_t) - 'adaptive' keeps the learning rate constant to 'learning_rate_init' as long as training loss keeps decreasing. Each time two consecutive epochs fail to decrease training loss by at least tol, or fail to increase validation score by at least tol if 'early_stopping' is on, the current learning rate is divided by 5. Only used when ``solver='sgd'``.	'constant'
	learning_rate_init learning_rate_init: float, default=0.001 The initial learning rate used. It controls the step-size in updating the weights. Only used when solver='sgd' or 'adam'.	0.001
	power_t power_t: float, default=0.5 The exponent for inverse scaling learning rate. It is used in updating effective learning rate when the learning_rate is set to 'invscaling'. Only used when solver='sgd'.	0.5
	shuffle shuffle: bool, default=True Whether to shuffle samples in each iteration. Only used when solver='sgd' or 'adam'.	True
	tol tol: float, default=1e-4 Tolerance for the optimization. When the loss or score is not improving by at least ``tol`` for ``n_iter_no_change`` consecutive iterations, unless ``learning_rate`` is set to 'adaptive', convergence is considered to be reached and training stops.	0.0001
	verbose verbose: bool, default=False Whether to print progress messages to stdout.	False
	warm_start warm_start: bool, default=False When set to True, reuse the solution of the previous call to fit as initialization, otherwise, just erase the previous solution. See :term:`the Glossary <warm_start>`.	False
	momentum momentum: float, default=0.9 Momentum for gradient descent update. Should be between 0 and 1. Only used when solver='sgd'.	0.9
	nesterovs_momentum nesterovs_momentum: bool, default=True Whether to use Nesterov's momentum. Only used when solver='sgd' and momentum > 0.	True
	early_stopping early_stopping: bool, default=False Whether to use early stopping to terminate training when validation score is not improving. If set to True, it will automatically set aside ``validation_fraction`` of training data as validation and terminate training when validation score is not improving by at least ``tol`` for ``n_iter_no_change`` consecutive epochs. The split is stratified, except in a multilabel setting. If early stopping is False, then the training stops when the training loss does not improve by more than ``tol`` for ``n_iter_no_change`` consecutive passes over the training set. Only effective when solver='sgd' or 'adam'.	False
	validation_fraction validation_fraction: float, default=0.1 The proportion of training data to set aside as validation set for early stopping. Must be between 0 and 1. Only used if early_stopping is True.	0.1
	beta_1 beta_1: float, default=0.9 Exponential decay rate for estimates of first moment vector in adam, should be in [0, 1). Only used when solver='adam'.	0.9
	beta_2 beta_2: float, default=0.999 Exponential decay rate for estimates of second moment vector in adam, should be in [0, 1). Only used when solver='adam'.	0.999
	epsilon epsilon: float, default=1e-8 Value for numerical stability in adam. Only used when solver='adam'.	1e-08
	n_iter_no_change n_iter_no_change: int, default=10 Maximum number of epochs to not meet ``tol`` improvement. Only effective when solver='sgd' or 'adam'. .. versionadded:: 0.20	10
	max_fun max_fun: int, default=15000 Only used when solver='lbfgs'. Maximum number of loss function calls. The solver iterates until convergence (determined by 'tol'), number of iterations reaches max_iter, or this number of loss function calls. Note that number of loss function calls will be greater than or equal to the number of iterations for the `MLPClassifier`. .. versionadded:: 0.22	15000

Fitted attributes

Name	Type	Value
best_loss_ best_loss_: float or None The minimum loss reached by the solver throughout fitting. If `early_stopping=True`, this attribute is set to `None`. Refer to the `best_validation_score_` fitted attribute instead.	float64	0.3084
best_validation_score_ best_validation_score_: float or None The best validation score (i.e. accuracy score) that triggered the early stopping. Only available if `early_stopping=True`, otherwise the attribute is set to `None`.	NoneType	None
classes_ classes_: ndarray or list of ndarray of shape (n_classes,) Class labels for each output.	ndarray[int64](2,)	[0,1]
coefs_ coefs_: list of shape (n_layers - 1,) The ith element in the list represents the weight matrix corresponding to layer i.	list	[array([[-9.47...7182985e-02]]), array([[-8.04...4735727e-02]])]
intercepts_ intercepts_: list of shape (n_layers - 1,) The ith element in the list represents the bias vector corresponding to layer i + 1.	list	[array([-0.021... -0.25089974]), array([0.04181039])]
loss_ loss_: float The current loss computed with the loss function.	float64	0.3086
loss_curve_ loss_curve_: list of shape (`n_iter_`,) The ith element in the list represents the loss at the ith iteration.	list	[np.float64(0.6887701418956516), np.float64(0.6599292705502692), np.float64(0.6365389999524382), np.float64(0.6155448000404404), ...]
n_features_in_ n_features_in_: int Number of features seen during :term:`fit`. .. versionadded:: 0.24	int	2
n_iter_ n_iter_: int The number of iterations the solver has run.	int	145
n_layers_ n_layers_: int Number of layers.	int	3
n_outputs_ n_outputs_: int Number of outputs.	int	1
out_activation_ out_activation_: str Name of the output activation function.	str	'lo...ic'
t_ t_: int The number of training samples seen by the solver during fitting.	int	507500
validation_scores_ validation_scores_: list of shape (`n_iter_`,) or None The score at each iteration on a held-out validation set. The score reported is the accuracy score. Only available if `early_stopping=True`, otherwise the attribute is set to `None`.	NoneType	None

Defining the custom risk¶

Any metric that can be written as sum(occurrence if condition) / sum(condition) can be controlled by the BinaryClassificationController, thanks to the Learn Then Test framework implemented in MAPIE.

A BinaryRisk is therefore defined by two per-sample functions, both taking the ground-truth labels y_true and the predictions y_pred and returning boolean arrays:

risk_condition: which samples count towards the metric (the denominator),
risk_occurrence: among those, which ones are a "success" (the numerator).

For specificity (true negative rate), we look at the actual negatives (y_true == 0, the condition) and count those that are correctly predicted as negative (y_pred == 0, the occurrence). In other words, specificity = count(y_pred == 0 and y_true == 0) / count(y_true == 0).

Because a higher specificity is better, we set higher_is_better=True so that MAPIE treats it as a performance metric rather than a risk.

specificity = BinaryRisk(
    risk_occurrence=lambda y_true, y_pred: y_pred == 0,
    risk_condition=lambda y_true, y_pred: y_true == 0,
    higher_is_better=True,
)

Controlling the custom risk¶

We now initialize a BinaryClassificationController using the probability estimation function of the fitted estimator (clf.predict_proba), our custom specificity risk, a target level, and a confidence level. The controller thresholds the predicted probability of the positive class: a higher threshold predicts fewer positives, which mechanically increases specificity but lowers recall.

When a custom risk is used, best_predict_param_choice="auto" cannot guess a sensible secondary objective, so we specify one explicitly. Here we maximize recall: among all thresholds that guarantee the target specificity, the controller selects the one that catches the most positive cases.

target_specificity = 0.8
confidence_level = 0.9
bcc = BinaryClassificationController(
    clf.predict_proba,
    specificity,
    target_level=target_specificity,
    confidence_level=confidence_level,
    best_predict_param_choice="recall",
)
bcc.calibrate(X_calib, y_calib)

print(
    f"{len(bcc.valid_predict_params)} thresholds found that guarantee a specificity "
    f"of at least {target_specificity} with a confidence of {confidence_level}.\n"
    "Among those, the one that maximizes the secondary objective (recall here) is: "
    f"{bcc.best_predict_param:.2f}."
)

Out:

60 thresholds found that guarantee a specificity of at least 0.8 with a confidence of 0.9.
Among those, the one that maximizes the secondary objective (recall here) is: 0.40.

Just like with a predefined risk, we can visualize how the threshold values impact specificity and which thresholds are statistically guaranteed.

proba_positive_class = clf.predict_proba(X_calib)[:, 1]

tested_thresholds = bcc._predict_params
specificities = np.full(len(tested_thresholds), np.inf)
for i, threshold in enumerate(tested_thresholds):
    y_pred = (proba_positive_class >= threshold).astype(int)
    # specificity is the recall of the negative class
    specificities[i] = recall_score(y_calib, y_pred, pos_label=0)

valid_thresholds_indices = np.array(
    [t in bcc.valid_predict_params for t in tested_thresholds]
)
best_threshold_index = np.where(tested_thresholds == bcc.best_predict_param)[0][0]

plt.figure()
plt.scatter(
    tested_thresholds[valid_thresholds_indices],
    specificities[valid_thresholds_indices],
    c="tab:green",
    label="Valid thresholds",
)
plt.scatter(
    tested_thresholds[~valid_thresholds_indices],
    specificities[~valid_thresholds_indices],
    c="tab:red",
    label="Invalid thresholds",
)
plt.scatter(
    tested_thresholds[best_threshold_index],
    specificities[best_threshold_index],
    c="tab:green",
    label="Best threshold",
    marker="*",
    edgecolors="k",
    s=300,
)
plt.axhline(target_specificity, color="tab:gray", linestyle="--")
plt.text(
    0.05,
    target_specificity + 0.02,
    "Target specificity",
    color="tab:gray",
    fontstyle="italic",
)
plt.xlabel("Threshold")
plt.ylabel("Specificity")
plt.legend()
plt.show()

plot risk control custom risk

Like in the quickstart, low thresholds correspond to specificity values below the target and are therefore invalid. Some thresholds whose empirical specificity is above the target are also rejected: risk control takes a margin to account for the uncertainty due to the finite size of the calibration set, so that the guarantee holds on unseen data with high probability.

We can confirm this on the test set by using the predict method of the controller, which applies the best threshold.

y_pred_test = bcc.predict(X_test)
test_specificity = recall_score(y_test, y_pred_test, pos_label=0)
test_recall = recall_score(y_test, y_pred_test, pos_label=1)

print(
    "With the best threshold, on the test set:\n"
    f"- specificity is {test_specificity:.3f} (target was {target_specificity}),\n"
    f"- recall is {test_recall:.3f}."
)

Out:

With the best threshold, on the test set:
- specificity is 0.809 (target was 0.8),
- recall is 0.881.

The specificity measured on the test set is close to (and above) the target, illustrating that the guarantee obtained on the calibration set transfers to unseen data. The same recipe applies to any binary metric that can be written as sum(occurrence if condition) / sum(condition): define it once with BinaryRisk, then control it with the BinaryClassificationController.

Total running time of the script: ( 0 minutes 1.203 seconds)

Download Python source code: plot_risk_control_custom_risk.py

Download Jupyter notebook: plot_risk_control_custom_risk.ipynb

Gallery generated by mkdocs-gallery