<jupyter_start><jupyter_text>*Studying and mitigating the effects of data drifts on ML model performance* Example notebook: Recalibration strategy illustrated at one example endpoint This notebook is part of the supporting information to the manuscript entitled: *Studying and mitigating the effects of data drifts on ML model performance at the example of chemical toxicity data** Authors: <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME>* Last updated: September 2021 In this notebook, the use of the Conformal prediction (CP) classes to study and mitigate the effects of data drifts on model calibration are demonstrated. The ChEMBL dataset, as used in the manuscript, for endpoint "CHEMBL228" is loaded as example. You can also exchange it with your own data.  Import packages and modules<jupyter_code>import pandas as pd
import numpy as np
import sys
import tarfile
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsRegressor
from nonconformist.nc import InverseProbabilityErrFunc, NcFactory
from continuous_calibration import (
StratifiedRatioSampler, BalancedStratifiedRatioSampler, CrossValidationSampler, InductiveConformalPredictor, 
    ContinuousCalibrationAggregatedConformalPredictor, CrossValidator
)
from continuous_calibration_helper_functions import draw_scatter_plot_one_endpoint, create_umap_data, plot_umap<jupyter_output><empty_output><jupyter_text>Define paths and parameters Here, we load the example dataset `ChEMBL228` , which is provided within this GitHub repository.Note that the full dataset with all 12 ChEMBL endpoints is due to its size not directly included, but can be downloaded from Zenodo ([link](https://zenodo.org/record/5167636)). Copy it to the `data` folder if you wish to work on a different endpoint (more info below).<jupyter_code># Dataset for ChEMBL228 endpoint
data_path = "../data/CHEMBL228_chembio_normalizedDesc.csv.tar.bz2"<jupyter_output><empty_output><jupyter_text>Dataset defining the splits by publication year into the CP subsets.Contains information for all 12 endpoints.<jupyter_code>time_split_threshold_path = "../data/data_size_chembio_chembl.csv"<jupyter_output><empty_output><jupyter_text>Choose an endpoint. Note that due to the size of the datasets, only one dataset for *ChEMBL228* endpoint is available from GitHub.Data for other ChEMBL endpoints can be downloaded from Zenodo ([link](https://zenodo.org/record/5167636)).<jupyter_code># Selected endpoint
endpoint = "CHEMBL228"<jupyter_output><empty_output><jupyter_text>Define parameters for machine learning and conformal prediction.* Note that the calculations take time and computational resources. If you just want to try out the notebook, we recommend using the small values provided. * To get the same plots as in the manuscript (maybe small deviations due to random splitting and random forest), replace them by the larger values (as provided in the comments).<jupyter_code>ntrees = 50  # 500
n_folds_acp = 3  # 20
# Defines how to split training data into calibration and proper
# training set using stratified spliting. Default is test_ratio=0.3,
# i.e. allocate 30% of data to calibration set
ratio_sampler = StratifiedRatioSampler(n_folds=n_folds_acp)<jupyter_output><empty_output><jupyter_text>Define UMAP parameters for compound descriptor space analysis.<jupyter_code>descriptors = "chembio"
n_neighbors = 100
min_distances = 0.8
distance_metric = "euclidean"
umap_colours = ["navy", "magenta", "orangered", "green"]<jupyter_output><empty_output><jupyter_text>Load and preprocess the data  Load data for ChEMBL228 endpoint Due to the size of the datasets, which contain the precalculated chembio descriptors, the example data for ChEMBL228 endpoint was deposited on GitHub as compressed file. The data contains:* Molecule ChEMBL ID* SMILES* Publication date* Binary activity (0 or 1)* Precalculated ChemBio descriptors<jupyter_code>data = pd.read_csv(data_path, compression='bz2', index_col=0)
data.head()<jupyter_output><empty_output><jupyter_text>Load data for another endpoint If you want to run the example notebook for any other ChEMBL endpoint, a compressed file with all ChEMBL datasets used in this work can be downloaded from [Zenodo](https://zenodo.org/record/5167636) as tar.bz2 files. * Select the correct `data_path` * Change the `endpoint` * Unpac the respective file.To do so uncomment convert the following four cells from `Raw NBConvert` to `Code`, so they will be executed..<jupyter_code># Update path
data_path = "../data/chembl_chembio_descriptors.tar.bz2"
# Update selected endpoint
endpoint = "CHEMBL220"
# Due to the size of the datasets containing the chembio descriptors, 
# they were uploaded to GitHub as tar.bz2 files and need to be unpacked first
tar = tarfile.open(data_path)
# The data is collected in the `files` dict, containing as keys the ChEMBL endpoint name and as value the dataframe
files = {}
for i, name in zip(tar, tar.getnames()):
    file = tar.extractfile(i)
    print(file)
    if file:  # Extraction also yields an empty folder (?) which will be ignored
        n = name.split("/")[1].split('_')[0]  # Get ChEMBL endpoint name from filename
        if n == endpoint:
            print(n)
            df = pd.read_csv(file)  # Read df
            files[n] = df  # Store in dict
# Select a specific example dataset
data = files[endpoint]<jupyter_output><empty_output><jupyter_text>Define arrays with labels, descriptors and publication years Drop rows without any information on the publication date.<jupyter_code># The publication date is required for temporal data splitting.
data.dropna(subset=["year"], inplace=True)<jupyter_output><empty_output><jupyter_text>Define the labels for the ML experiment, i.e., the measured binarized bioactivities.<jupyter_code>y = data[f"{endpoint}_bioactivity"].values
y<jupyter_output><empty_output><jupyter_text>Define feature columns from the ChemBio descriptor, i.e., p-values and chemical descriptors.<jupyter_code># Any other columns present in the dataset, but not part of the descriptors, are explicitly excluded.

columns = [
            col
            for col in data.columns
            if (not col.startswith("Toxicity"))  # Exclude in vivo labels used for p-values
            and (col != f"{endpoint}_bioactivity")  # Exclude ChEMBL bioactivity labels
            and (col != "SMILES (Canonical)")  # Exclude canonical smiles
            and (col != "smiles")  # Exclude smiles
            and (col != "year")  # Exclude publication year
            and (col != "molecule_chembl_id")  # Exclude molecule ChEMBL id
        ]

X = data[columns].values
X<jupyter_output><empty_output><jupyter_text>Collect the publication years in an arry. This information will be used for temporal data splitting.<jupyter_code>years = data["year"].values
years<jupyter_output><empty_output><jupyter_text>Split the data Here, we implement a time-split to separate the data in four subsets based on the publication year:* training set* update 1 set* update 2 set* holdout set  Load the dataframe defining time-split thresholds per endpointThe `time_split_threshold_path` variable holds the precalculated information which year is used per data set to create the time split data while retaining specific ratios (see manuscript for more details).<jupyter_code>splits_df = pd.read_csv(time_split_threshold_path, index_col=0, 
                        usecols=["chembl_id", "train_thresh", "update1_thresh", "update2_thresh"])
splits_df
splits_df.index<jupyter_output><empty_output><jupyter_text>Get year thresholds to split the data for the selected endpoint.<jupyter_code>thresholds = splits_df["train_thresh"][endpoint], splits_df["update1_thresh"][endpoint], splits_df["update2_thresh"][endpoint]
thresholds<jupyter_output><empty_output><jupyter_text>Do the splitting based on the threshold.<jupyter_code># Collect the respective indices based on the threshold
mask_train = years <= thresholds[0]
mask_update1 = (years > thresholds[0]) & (years <= thresholds[1])
mask_update2 = (years > thresholds[1]) & (years <= thresholds[2])
mask_holdout = years > thresholds[2]

# Split the data accordingly
X_train, y_train = X[mask_train], y[mask_train]
X_update1, y_update1 = X[mask_update1], y[mask_update1]
X_update2, y_update2 = X[mask_update2], y[mask_update2]
X_holdout, y_holdout = X[mask_holdout], y[mask_holdout]<jupyter_output><empty_output><jupyter_text>Inspect the data subsets.<jupyter_code># Check shape of the arrays
print(f"Training set shape: y={y_train.shape}, X={X_train.shape}")
print(f"Update1 set shape: y={y_update1.shape}, X={X_update1.shape}")
print(f"Update2 set shape: y={y_update2.shape}, X={X_update2.shape}")
print(f"Holdout set shape: y={y_holdout.shape}, X={X_holdout.shape}")
# Check if we have actives and inactives in all datasets
print(f"Active compounds in \n\t training set: {int(y_train.sum())}")
print(f"\t update1 set: {int(y_update1.sum())}\n\t update2 set: {int(y_update2.sum())}\n\t holdout set: {int(y_holdout.sum())}")<jupyter_output>Active compounds in 
	 training set: 858
	 update1 set: 373
	 update2 set: 235
	 holdout set: 196<jupyter_text>Conformal prediction recalibration strategy without cross-validation (CV) * The experiments in the original study were performed in a fivefold CV. * For illustration of the recalibration strategy, it is here explained, step by step, without CV first. * Note that the analysis of the models and predictions will be presented after the `crossvalidation` part.  Prepare Aggregated Conformal Predictor (ACP)<jupyter_code>clf = RandomForestClassifier(n_estimators=ntrees)
error_function = InverseProbabilityErrFunc()
nc = NcFactory.create_nc(
        clf, err_func=error_function, normalizer_model=None
    )
icp = InductiveConformalPredictor(
        nc_function=nc, smoothing=False, condition=(lambda instance: instance[1])
    )  # Mondrian as (default) condition)
acp = ContinuousCalibrationAggregatedConformalPredictor(predictor=icp, sampler=ratio_sampler,
        aggregation_func=np.median
    )<jupyter_output><empty_output><jupyter_text>Fit and calibrate the model, make predictions for holdout set Calibrate with the `original` calibration set, which is a split from X_train/y_train.<jupyter_code>acp.fit_calibrate(X_train=X_train, y_train=y_train)<jupyter_output><empty_output><jupyter_text>Make predictions with the model and original calibration set on the holdout set.<jupyter_code>original_prediction = acp.predict(X_test=X_holdout)

original_prediction[:5]<jupyter_output><empty_output><jupyter_text>Update calibration set I Calibrate the fitted (not yet calibrated) predictors with a 'new' dataset<jupyter_code>acp.calibrate_update(X_update=X_update1, y_update=y_update1)<jupyter_output><empty_output><jupyter_text>Make predictions with the model and the updated calibration set 1 on the holdout set.<jupyter_code>calibrate_update_prediction = acp.predict_calibrate_update(updated_number=1, X_test=X_holdout)
calibrate_update_prediction[:5]<jupyter_output><empty_output><jupyter_text>Update calibration set II We can update with an even 'newer' calibration set (as often as we want). The model with the updated calibration set will be saved in the acp object.<jupyter_code># Check how many updated calibration sets we already have. 
# Here, this should be 1.
print(acp.predictors_calibrated_update.keys())
acp.calibrate_update(X_update=X_update2, y_update=y_update2)<jupyter_output><empty_output><jupyter_text>Make predictions with the model and the updated calibration set 2 on the holdout set.<jupyter_code>calibrate_update_2_predcition = acp.predict_calibrate_update(updated_number=2, X_test=X_holdout)
calibrate_update_2_predcition[:5]<jupyter_output><empty_output><jupyter_text>Conformal prediction recalibration strategy within cross-validation (CV) Instead of using the individual acp's, we might want to include everything into a cross-validator. Advantages: * Splitting of CV into test and train set are handled* Fixed and reproducible pipeline (except for random splitting and random forest)* Updating of the same trained models possibleNote that we use fivefold CV by default.  Prepare ACP and cross-validator class<jupyter_code>cv_acp = ContinuousCalibrationAggregatedConformalPredictor(predictor=icp, sampler=ratio_sampler,
        aggregation_func=np.median
    )
cross_validator = CrossValidator(predictor=cv_acp, cv_splitter=CrossValidationSampler())<jupyter_output><empty_output><jupyter_text>Fit and calibrate the model, make predictions for holdout set Predictions on the holdout set.Crossvalidate with the "original" training set (no update calibration set)<jupyter_code>cv_eval_df = cross_validator.cross_validate(X_train=X_train, y_train=y_train, X_test=X_holdout, y_test=y_holdout, steps=10, endpoint=endpoint)<jupyter_output><empty_output><jupyter_text>Analyse internal cross-validation (with original calibration set) * Evaluation per CV fold.* Note that the metrics are calculated for significant levels in the range of [0,1] (stepsize 0.1). Thus, you receive 11 entries per fold.<jupyter_code>cv_eval_df[["validity", "efficiency", "accuracy", "significance_level", "fold"]].set_index("fold").head()<jupyter_output><empty_output><jupyter_text>Averaged (all folds) evaluations, again per significance level<jupyter_code>cross_validator.averaged_evaluation_df_cv[["validity mean", "efficiency mean", "accuracy mean"]]<jupyter_output><empty_output><jupyter_text>Since often the datasets are very unbalanced w.r.t. the ration of active and inactive molecules, balanced metrics are calculated as well.Balanced metrics, averaged (all folds) evaluations, again per significant level<jupyter_code>cross_validator.averaged_evaluation_df_cv[["validity_bal mean", "efficiency_bal mean", "accuracy_bal mean"]]<jupyter_output><empty_output><jupyter_text>Analyse performance using the CV-models when predicting the holdout set Note the holdout set was already annotated when calling `cross_validator.cross_validate`.  The `cross_validator` object stores the predictions in a dictionary that contains the information for the predictions within the cross-validation and the predictions for the test (holdout) set. The predictions (values) are in the form of a list, list elements are the predicted p-values and the known labels.<jupyter_code>print("keys: ", cross_validator._predictions.keys(), "\n")
print("p-values: ", cross_validator._predictions["pred_test"][0][0][:5], "\n")
# [0]: predictions, [0]: first loop of cv, [:5]: first five p-values
print("labels: ", cross_validator._predictions["pred_test"][1])
# [1]: labels<jupyter_output>keys:  dict_keys(['cv', 'pred_test']) 

p-values:  [[0.38983051 0.0410628 ]
 [0.33050847 0.04830918]
 [0.33898305 0.04589372]
 [0.27118644 0.05797101]
 [0.16949153 0.11111111]] 

labels:  [[0. 0. 1. ... 0. 0. 1.]
 [0. 0. 1. ... 0. 0. 1.]
 [0. 0. 1. ... 0. 0. 1.]
 [0. 0. 1. ... 0. 0. 1.]
 [0. 0. 1. ... 0. 0. 1.]]<jupyter_text>Evaluation on the holdout set (averaged over all folds)<jupyter_code>cross_validator.averaged_evaluation_df_pred_test[["validity mean", "efficiency mean", "accuracy mean"]]<jupyter_output><empty_output><jupyter_text>Update the calibration set within the crossvalidator I<jupyter_code>cross_validator.cross_validate_calibrate_update(X_update=X_update1, y_update=y_update1, X_test=X_holdout, y_test=y_holdout, steps=10, endpoint=endpoint)<jupyter_output><empty_output><jupyter_text>Evaluating the predictions on the holdout set using the update I calibration set (averaged over all folds)<jupyter_code>cross_validator.averaged_evaluation_df_cal_update_1[["validity mean", "efficiency mean", "accuracy mean"]]<jupyter_output><empty_output><jupyter_text>Update the calibration set within the crossvalidator II<jupyter_code>cross_validator.cross_validate_calibrate_update(X_update=X_update2, y_update=y_update2, X_test=X_holdout, y_test=y_holdout, steps=10, endpoint=endpoint)<jupyter_output><empty_output><jupyter_text>Evaluating the predictions on the holdout set using the update II calibration set (averaged over all folds)<jupyter_code>cross_validator.averaged_evaluation_df_cal_update_2[["validity mean", "efficiency mean", "accuracy mean"]]<jupyter_output><empty_output><jupyter_text>Visualisations  Scatter plot to compare different strategies To compare balanced validity, balanced efficiency and balanced accuracy for the three strategies, they are plotted in a scatter plot. In the following cells, the strategies and evaluation measures to be plotted are defined.<jupyter_code>strategies = ["cv_original", "original", "update1"]
eval_df_columns = ["validity_bal mean", "efficiency_bal mean", "accuracy_bal mean", "significance_level"]
evaluation_measures = [m.split(" ")[0] for m in eval_df_columns[:-1]]
print(evaluation_measures)
evaluation_dfs_dict = {}
evaluation_dfs_dict["cv_original"] = cross_validator.averaged_evaluation_df_cv[eval_df_columns]
evaluation_dfs_dict["original"] = cross_validator.averaged_evaluation_df_pred_test[eval_df_columns]
evaluation_dfs_dict["update1"] = cross_validator.averaged_evaluation_df_cal_update_1[eval_df_columns]
evaluation_dfs_dict["update2"] = cross_validator.averaged_evaluation_df_cal_update_2[eval_df_columns]<jupyter_output><empty_output><jupyter_text>Draw a scatter plot that compares balanced validity, efficiency and accuracy within the CV, and predicting the test set using the original and updated calibration sets. The dashed line represents the expected validity (80%) at a significance level of 0.2.<jupyter_code>draw_scatter_plot_one_endpoint(endpoint, evaluation_dfs_dict, evaluation_measures=evaluation_measures)<jupyter_output><empty_output><jupyter_text>UMAP to analyse descriptor space of subsets A UMAP can be helpful to analyse the chemical space of the subsets. Ideally, the update and holdout sets cover a similar descriptor space.* In the first step, we create an embedding. * In a second step, this will be plotted.<jupyter_code>embedding, colours_numbers = create_umap_data(
    chembl_id=endpoint, descriptors_df=data, thresholds=thresholds, n_neighbors=n_neighbors,
    min_distances=min_distances, distance_metric=distance_metric, dataset_colours=umap_colours)
umap_plot = plot_umap(embedding, endpoint, colours_datasets=colours_numbers, umap_colours=umap_colours, n_neighbor=n_neighbors, min_distance=min_distances, distance_metric=distance_metric)<jupyter_output>2111
{'fraction_non_alphanumeric': 0.08612153854523692, 'fraction_numerical': 0.016593221261084815, 'mean_word_length': 6.62106135986733, 'pattern_counts': {'":': 0, '<': 156, '<?xml version=': 0, '>': 156, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '37', 'dir': 'jupyter-structured-clean-dedup', 'id': '370918', 'license': 'MIT', 'n_tokens_mistral': 5897, 'n_tokens_neox': 5506, 'n_words': 1728, 'path': 'code/1_continuous_calibration_example.ipynb', 'repo_name': 'volkamerlab/CPrecalibration_manuscript_SI', 'size': '18421'}