4 Model Evaluation

A model is only as useful as its ability to generalize to new, unseen data. This chapter covers the key principles of evaluating machine learning models in the context of chemistry, where naive evaluation strategies can paint a dangerously optimistic picture.

4.1 Train / Validation / Test Philosophy

The fundamental principle: never evaluate a model on the data it was trained on. To get an honest estimate of how well a model will perform in practice, we split our data into distinct subsets:

Training set (~60–70%): Used to fit the model parameters.
Validation set (~15–20%): Used to tune hyperparameters (e.g., regularization strength, polynomial degree) and select among competing models.
Test set (~15–20%): Held out until the very end. Touched exactly once to report final, unbiased performance.

The validation set acts as a proxy for unseen data during model development. If we used the test set to choose hyperparameters, we would effectively be “training” on it and our final performance estimate would be optimistic.

A common mistake is to repeatedly peek at the test set during development. Even if the test set is never directly used for gradient updates, selecting the model that performs best on it causes information leakage — the test set silently guides our modeling choices.

4.1.1 Interactive example: train/validation/test split

Adjust the split ratios and the noise level to see how the amount of training data and the split proportions affect model performance across the three subsets. Only training and validation are set directly; the test set takes whatever remains, and the plot title reports the resulting three-way split.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 600

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({
    "font.family": "serif",
    "figure.facecolor": "#fafaf7",
    "axes.facecolor": "#fafaf7",
    "axes.grid": True,
    "grid.alpha": 0.3,
})

app_ui = ui.page_fluid(
    ui.layout_columns(
        ui.input_slider("train_frac", "Training fraction", min=0.3, max=0.8, value=0.6, step=0.05),
        ui.input_slider("val_frac", "Validation fraction", min=0.05, max=0.65, value=0.2, step=0.05),
        ui.input_slider("noise", "Noise level", min=0.1, max=2.0, value=0.8, step=0.1),
        ui.input_slider("degree", "Polynomial degree", min=1, max=12, value=3, step=1),
        ui.input_slider("n_points", "Total data points", min=30, max=200, value=80, step=10),
        col_widths=(4, 4, 4, 6, 6),
    ),
    ui.output_plot("split_plot"),
)

def server(input, output, session):
    @reactive.effect
    @reactive.event(input.train_frac)
    def _limit_val_frac():
        # Cap the validation slider at what the training fraction leaves behind,
        # keeping at least 5% for the test set. Without this the slider can be
        # set to a value the plot cannot honour.
        cap = round(1.0 - input.train_frac() - 0.05, 2)
        ui.update_slider("val_frac", max=cap, value=min(input.val_frac(), cap))

    @output
    @render.plot
    def split_plot():
        np.random.seed(42)
        n = input.n_points()
        x = np.sort(np.random.uniform(-3, 3, n))
        y_true = np.sin(x) + 0.3 * x
        y = y_true + np.random.normal(0, input.noise(), n)

        # The test set is whatever is left over, so the validation fraction is
        # capped to always leave at least 5% of the points for testing. Without
        # this the test set can be empty and its score comes out as nan.
        train_frac = input.train_frac()
        val_frac = min(input.val_frac(), 1.0 - train_frac - 0.05)
        test_frac = 1.0 - train_frac - val_frac

        train_end = int(n * train_frac)
        val_end = train_end + int(n * val_frac)

        idx = np.random.permutation(n)
        train_idx = idx[:train_end]
        val_idx = idx[train_end:val_end]
        test_idx = idx[val_end:]

        x_train, y_train = x[train_idx], y[train_idx]
        x_val, y_val = x[val_idx], y[val_idx]
        x_test, y_test = x[test_idx], y[test_idx]

        deg = input.degree()
        coeffs = np.polyfit(x_train, y_train, deg)
        poly = np.poly1d(coeffs)

        x_line = np.linspace(-3.2, 3.2, 300)
        y_pred_line = poly(x_line)

        mse_train = np.mean((y_train - poly(x_train))**2)
        mse_val = np.mean((y_val - poly(x_val))**2) if len(x_val) > 0 else float("nan")
        mse_test = np.mean((y_test - poly(x_test))**2) if len(x_test) > 0 else float("nan")

        fig, ax = plt.subplots(figsize=(6, 4))
        ax.scatter(x_train, y_train, c="#2b6cb0", s=30, alpha=0.7, label=f"Train (MSE={mse_train:.3f})", zorder=3)
        ax.scatter(x_val, y_val, c="#e8833a", s=30, alpha=0.7, label=f"Val (MSE={mse_val:.3f})", zorder=3)
        ax.scatter(x_test, y_test, c="#c53030", s=30, alpha=0.7, label=f"Test (MSE={mse_test:.3f})", zorder=3)
        ax.plot(x_line, np.sin(x_line) + 0.3 * x_line, "--", color="#999", linewidth=1.2, label="True function")
        y_clipped = np.clip(y_pred_line, -6, 6)
        ax.plot(x_line, y_clipped, color="#276749", linewidth=2, label=f"Degree {deg} fit")
        ax.set_ylim(-5, 5)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_title(f"Train {train_frac:.0%} / Val {val_frac:.0%} / Test {test_frac:.0%}")
        ax.legend(fontsize=9, loc="upper left")
        fig.tight_layout()
        return fig

app = App(app_ui, server)

Notice how the test MSE gives a more honest picture of generalization than the training MSE, especially for high-degree polynomials where the model memorizes the training data.

4.2 K-Fold Cross-Validation

When data is limited (a common situation in chemistry), setting aside a fixed validation set wastes precious training examples. K-fold cross-validation offers a solution: partition the data into $k$ equally sized folds, then train $k$ separate models, each time holding out one fold as the validation set and training on the remaining $k - 1$ folds.

The cross-validated score is the average performance over the $k$ held-out folds:

\[ \text{CV}(k) = \frac{1}{k} \sum_{i=1}^{k} \mathcal{L}_i \]

where $\mathcal{L}_i$ is the loss on fold $i$. Common choices are $k = 5$ or $k = 10$. The extreme case $k = n$ is leave-one-out cross-validation (LOOCV), which has low bias but high variance and is computationally expensive. (The bias–variance trade-off is introduced in the linear models chapter.)

Advantages:

Every data point is used for both training and validation.
Yields $k$ scores rather than a single number, so the result comes with an uncertainty. We interpret the standard error of the mean , $\sigma/\sqrt{k}$, where $\sigma$ is the standard deviation across the $k$ folds, to give a sense of the variability in the cross-validated performance.

Disadvantages:

$k$ times more expensive than a single train/val split.
Folds can be correlated, especially for LOOCV.

4.3 Metrics

Choosing the right metric is essential — a model that looks excellent by one metric may be mediocre by another. The appropriate metric depends on whether the task is regression or classification.

4.3.1 Regression metrics

Root Mean Squared Error (RMSE): Penalizes large errors more heavily due to squaring.

\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \]

Mean Absolute Error (MAE): More robust to outliers; gives equal weight to all errors.

\[ \text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i| \]

Coefficient of Determination ($R^2$): Measures the proportion of variance explained. $R^2 = 1$ is perfect; $R^2 = 0$ means the model is no better than predicting the mean; $R^2 < 0$ means the model is worse than the mean.

\[ R^2 = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2} \]

4.3.2 Interactive example: regression metrics

Drag the noise slider and add an outlier to see how RMSE, MAE, and $R^2$ respond differently.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 520

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({
    "font.family": "serif",
    "figure.facecolor": "#fafaf7",
    "axes.facecolor": "#fafaf7",
    "axes.grid": True,
    "grid.alpha": 0.3,
})

app_ui = ui.page_fluid(
    ui.layout_columns(
        ui.input_slider("noise_met", "Noise level", min=0.1, max=3.0, value=0.5, step=0.1),
        ui.input_slider("outlier_size", "Outlier magnitude", min=0.0, max=8.0, value=0.0, step=0.5),
        ui.input_slider("n_pts_met", "Number of points", min=20, max=100, value=40, step=5),
    ),
    ui.output_plot("metrics_plot"),
)

def server(input, output, session):
    @reactive.calc
    def compute_data():
        np.random.seed(21)
        n = input.n_pts_met()
        x = np.sort(np.random.uniform(0, 10, n))
        y_true = 0.8 * x + 1.0
        y = y_true + np.random.normal(0, input.noise_met(), n)

        # Add outlier at the last point
        if input.outlier_size() > 0:
            y[-1] = y_true[-1] + input.outlier_size()

        # Simple linear fit
        coeffs = np.polyfit(x, y, 1)
        y_pred = np.polyval(coeffs, x)

        rmse = np.sqrt(np.mean((y - y_pred)**2))
        mae = np.mean(np.abs(y - y_pred))
        ss_res = np.sum((y - y_pred)**2)
        ss_tot = np.sum((y - np.mean(y))**2)
        r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0

        return x, y, y_true, y_pred, coeffs, rmse, mae, r2

    @output
    @render.plot
    def metrics_plot():
        x, y, y_true, y_pred, coeffs, rmse, mae, r2 = compute_data()

        fig, ax = plt.subplots(figsize=(6, 4))
        ax.scatter(x, y, c="#2b6cb0", s=30, alpha=0.7, zorder=3, label="Data")
        x_line = np.linspace(-0.5, 10.5, 100)
        ax.plot(x_line, np.polyval(coeffs, x_line), color="#e8833a", linewidth=2, label="Linear fit")
        ax.plot(x_line, 0.8 * x_line + 1.0, "--", color="#999", linewidth=1, label="True line")

        # Draw residuals
        for xi, yi, yp in zip(x, y, y_pred):
            ax.plot([xi, xi], [yi, yp], color="#c53030", alpha=0.3, linewidth=1)

        if input.outlier_size() > 0:
            ax.scatter([x[-1]], [y[-1]], c="#c53030", s=80, zorder=5, marker="x", linewidths=2, label="Outlier")

        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.set_title("Regression Metrics Comparison")
        ax.legend(fontsize=9, loc="upper left")
        ax.text(0.98, 0.02,
                f"RMSE = {rmse:.3f}\nMAE  = {mae:.3f}\nR²   = {r2:.3f}",
                transform=ax.transAxes, ha="right", va="bottom",
                fontsize=10, family="monospace",
                bbox=dict(boxstyle="round,pad=0.4", facecolor="white", edgecolor="#ddd"))
        fig.tight_layout()
        return fig

app = App(app_ui, server)

Note that RMSE is always greater than or equal to MAE, and a single large outlier inflates RMSE far more than MAE.

4.3.3 Classification metrics

For binary classification tasks (e.g., active/inactive, toxic/non-toxic), we use different metrics.

Confusion matrix terminology:

	Predicted positive	Predicted negative
Actual positive	True Positive (TP)	False Negative (FN)
Actual negative	False Positive (FP)	True Negative (TN)

Accuracy $= \frac{TP + TN}{TP + TN + FP + FN}$ — misleading for imbalanced datasets.

Precision $= \frac{TP}{TP + FP}$ — of those predicted positive, how many are correct?

Recall (Sensitivity) $= \frac{TP}{TP + FN}$ — of actual positives, how many did we catch?

ROC-AUC: The area under the Receiver Operating Characteristic curve, which plots recall vs. false positive rate at various classification thresholds. AUC = 1.0 is perfect; AUC = 0.5 is random guessing.

4.3.4 Interactive example: ROC curve and threshold

Adjust the classification threshold and class separability to see how the ROC curve, precision, and recall change.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 700

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({
    "font.family": "serif",
    "figure.facecolor": "#fafaf7",
    "axes.facecolor": "#fafaf7",
    "axes.grid": True,
    "grid.alpha": 0.3,
})

app_ui = ui.page_fluid(
    ui.layout_columns(
        ui.input_slider("separability", "Class separability", min=0.5, max=4.0, value=2.0, step=0.25),
        ui.input_slider("threshold", "Classification threshold", min=0.0, max=1.0, value=0.5, step=0.02),
        ui.input_slider("imbalance", "Class imbalance (% positive)", min=10, max=90, value=50, step=5),
    ),
    ui.output_plot("roc_plot", height="560px"),
)

def server(input, output, session):
    @reactive.calc
    def compute_class_data():
        np.random.seed(33)
        n = 200
        n_pos = int(n * input.imbalance() / 100)
        n_neg = n - n_pos

        sep = input.separability()
        scores_pos = np.random.normal(sep / 2, 1, n_pos)
        scores_neg = np.random.normal(-sep / 2, 1, n_neg)

        # Convert to probabilities via sigmoid
        def sigmoid(x):
            return 1 / (1 + np.exp(-x))

        probs_pos = sigmoid(scores_pos)
        probs_neg = sigmoid(scores_neg)

        y_true = np.concatenate([np.ones(n_pos), np.zeros(n_neg)])
        y_probs = np.concatenate([probs_pos, probs_neg])

        return y_true, y_probs

    @output
    @render.plot
    def roc_plot():
        from matplotlib.gridspec import GridSpec
        from matplotlib.colors import ListedColormap

        y_true, y_probs = compute_class_data()
        thresh = input.threshold()

        # ROC curve over all thresholds
        thresholds = np.linspace(0, 1, 200)
        tpr_list = []
        fpr_list = []
        for t in thresholds:
            y_pred = (y_probs >= t).astype(int)
            tp_t = np.sum((y_pred == 1) & (y_true == 1))
            fn_t = np.sum((y_pred == 0) & (y_true == 1))
            fp_t = np.sum((y_pred == 1) & (y_true == 0))
            tn_t = np.sum((y_pred == 0) & (y_true == 0))
            tpr_list.append(tp_t / max(tp_t + fn_t, 1))
            fpr_list.append(fp_t / max(fp_t + tn_t, 1))

        sorted_idx = np.argsort(fpr_list)
        fpr_sorted = np.array(fpr_list)[sorted_idx]
        tpr_sorted = np.array(tpr_list)[sorted_idx]
        auc = np.trapz(tpr_sorted, fpr_sorted)

        # Confusion matrix at the current threshold
        y_pred_curr = (y_probs >= thresh).astype(int)
        tp = int(np.sum((y_pred_curr == 1) & (y_true == 1)))
        fn = int(np.sum((y_pred_curr == 0) & (y_true == 1)))
        fp = int(np.sum((y_pred_curr == 1) & (y_true == 0)))
        tn = int(np.sum((y_pred_curr == 0) & (y_true == 0)))
        curr_tpr = tp / max(tp + fn, 1)
        curr_fpr = fp / max(fp + tn, 1)

        accuracy = (tp + tn) / len(y_true)
        precision = tp / max(tp + fp, 1)
        recall = tp / max(tp + fn, 1)

        fig = plt.figure(figsize=(6, 6.2))
        gs = GridSpec(2, 2, figure=fig, height_ratios=[1, 1], hspace=0.55, wspace=0.4)
        ax = fig.add_subplot(gs[0, 0])
        ax2 = fig.add_subplot(gs[0, 1])
        ax_cm = fig.add_subplot(gs[1, 0])
        ax_met = fig.add_subplot(gs[1, 1])

        # Top-left: score distributions
        pos_probs = y_probs[y_true == 1]
        neg_probs = y_probs[y_true == 0]
        ax.hist(neg_probs, bins=25, alpha=0.6, color="#2b6cb0", label="Negative", density=True)
        ax.hist(pos_probs, bins=25, alpha=0.6, color="#c53030", label="Positive", density=True)
        ax.axvline(thresh, color="#e8833a", linewidth=2, linestyle="--", label=f"Threshold={thresh:.2f}")
        ax.set_xlabel("Predicted probability")
        ax.set_ylabel("Density")
        ax.set_title("Score distributions")
        ax.legend(fontsize=8)

        # Top-right: ROC curve
        ax2.plot(fpr_sorted, tpr_sorted, color="#2b6cb0", linewidth=2, label=f"ROC (AUC={auc:.3f})")
        ax2.plot([0, 1], [0, 1], "--", color="#999", linewidth=1, label="Random")
        ax2.scatter([curr_fpr], [curr_tpr], c="#e8833a", s=80, zorder=5, label="Current threshold")
        ax2.set_xlabel("False Positive Rate")
        ax2.set_ylabel("True Positive Rate")
        ax2.set_title("ROC curve")
        ax2.legend(fontsize=8, loc="lower right")
        ax2.set_xlim(-0.02, 1.02)
        ax2.set_ylim(-0.02, 1.02)

        # Bottom-left: confusion matrix at the current threshold
        cm_cmap = ListedColormap(["#f5c6c6", "#c6e6c6"])
        grid = np.array([[1, 0], [0, 1]])
        ax_cm.imshow(grid, cmap=cm_cmap, vmin=0, vmax=1)
        ax_cm.grid(False)
        cells = [[("TP", tp), ("FN", fn)], [("FP", fp), ("TN", tn)]]
        for r in range(2):
            for c in range(2):
                lbl, val = cells[r][c]
                ax_cm.text(c, r, f"{lbl}\n{val}", ha="center", va="center",
                           fontsize=12, fontweight="bold", color="#222")
        ax_cm.set_xticks([0, 1])
        ax_cm.set_xticklabels(["Pred +", "Pred -"])
        ax_cm.set_yticks([0, 1])
        ax_cm.set_yticklabels(["Actual +", "Actual -"])
        ax_cm.set_title("Confusion matrix")

        # Bottom-right: derived metrics
        ax_met.axis("off")
        ax_met.text(0.0, 0.9,
                    f"Accuracy   = {accuracy:.3f}\n\n"
                    f"Precision  = {precision:.3f}\n\n"
                    f"Recall     = {recall:.3f}",
                    transform=ax_met.transAxes, va="top", ha="left",
                    fontsize=12, family="monospace")

        return fig

app = App(app_ui, server)

Which metric to choose?

Use RMSE when large errors are particularly costly (e.g., toxicity prediction).
Use MAE when you want robustness to outliers (e.g., noisy experimental data).
Use $R^2$ for a quick sense of overall fit quality, but note that it is relative to the variance of your test set: the same model scores higher on a chemically diverse set than on a narrow one, so $R^2$ values are not comparable across datasets.
Use ROC-AUC for balanced classification. For imbalanced tasks, consider Precision-Recall AUC instead.

4.4 Overfitting and Underfitting

A metric tells you how well the model scored; it does not tell you why it scored that way. The two failure modes behind a poor score are opposites, and they call for opposite remedies.

Underfitting (high bias): The model is too simple to capture the underlying pattern, and performs poorly on held-out data because it never captured the signal in the first place.

Overfitting (high variance): The model has absorbed noise specific to the training set, and performs poorly on held-out data because that noise does not recur.

Both are diagnosed from the held-out score, and distinguished by how that score responds to model complexity: it improves with added complexity if you are underfitting, and degrades if you are overfitting.

The training-error gap is not a reliable diagnostic

It is tempting to diagnose overfitting from the gap between training and validation error, but for many models a low training error is structural rather than symptomatic: a random forest of fully grown trees drives training error to nearly zero by construction, without necessarily harming the held-out score. Judge the model by the held-out score and how it responds to complexity, not by the gap.

4.4.1 Remedies

Problem	Remedies
Underfitting (held-out score poor, improves with complexity)	Increase model complexity, add more features, reduce regularization, train longer
Overfitting (held-out score degrades as complexity grows)	Add regularization (L1/L2), reduce model complexity, get more data, use dropout, early stopping

Early stopping

For iterative models (neural networks, gradient boosting), a practical remedy for overfitting is early stopping: monitor the validation loss during training and stop when it begins to increase, even if the training loss is still decreasing.

4.5 How You Split Decides What You Measure

Everything so far has assumed that the test set tells you how the model will perform in use. That is only true if the test set resembles what the model will actually see. In most machine learning applications a random split makes this approximately true, because the data are independent and identically distributed (i.i.d.). In chemistry the assumption frequently breaks down: molecules in a dataset are rarely independent draws. They share scaffolds and functional groups, arrive in congeneric series (families of close analogues built around a common core), and are the product of a project that deliberately explored one region of chemical space at a time.

A random split therefore leaves near-duplicates on both sides, and the model can score well by recognising a scaffold it has already seen. But the deeper problem is not that similar molecules appear in both train and test: it is that the molecules you will actually predict on are dissimilar to both. What the test score fails to capture is out-of-distribution (OOD) generalization.

4.5.1 Splitting strategies as probes of OOD generalization

Each split withholds a different kind of chemistry from the model, and so simulates a different distribution shift. The right one is whichever shift you expect between the data you have and the molecules you will actually predict on.

Strategy	What the model never sees in training	Choose it when
Random	Nothing in particular; the test molecules resemble the training ones	The model will be used on the same chemical space the data came from, so the task really is interpolation
Cluster	Entire regions of feature space	The model will meet broadly new chemistry; cluster the fingerprints (k-means, Butina) and hold out whole clusters
Scaffold	The molecular core	The model will see new series rather than new analogues; Murcko scaffolds are the usual heuristic for “same core”, not the only one
Temporal	Whatever changed from one period to the next: reformulated assays, a move to a new chemotype, easy targets done first	The model will run on data generated after the training data (not merely “the future”, which is true of any prediction)

Two consequences: a reported score is only interpretable together with its split, and the aggressive splits are not simply better – a scaffold split gives a pessimistic number for a task that really is interpolation.

4.5.2 Interactive example: random vs. scaffold split

This demonstration shows how random splits allow structurally similar molecules to leak across train and test sets, while scaffold-aware splits keep related molecules together. Points represent molecules; color represents their scaffold group.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 560

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({
    "font.family": "serif",
    "figure.facecolor": "#fafaf7",
    "axes.facecolor": "#fafaf7",
    "axes.grid": True,
    "grid.alpha": 0.3,
})

app_ui = ui.page_fluid(
    ui.layout_columns(
        ui.input_radio_buttons("split_type", "Split strategy",
            choices={"random": "Random", "scaffold": "Scaffold-aware"},
            selected="random"),
        ui.input_slider("n_scaffolds", "Number of scaffold groups", min=3, max=8, value=5, step=1),
        ui.input_slider("test_frac_split", "Test fraction", min=0.1, max=0.4, value=0.2, step=0.05),
    ),
    ui.output_plot("split_compare_plot"),
)

def server(input, output, session):
    @output
    @render.plot
    def split_compare_plot():
        np.random.seed(12)
        n_scaffolds = input.n_scaffolds()
        pts_per_scaffold = 15
        n = n_scaffolds * pts_per_scaffold

        # Generate clustered 2D data (simulating scaffold groups)
        centers = np.random.uniform(-3, 3, (n_scaffolds, 2))
        x_all = []
        y_all = []
        groups = []
        for i, c in enumerate(centers):
            pts = np.random.normal(c, 0.5, (pts_per_scaffold, 2))
            x_all.append(pts[:, 0])
            y_all.append(pts[:, 1])
            groups.extend([i] * pts_per_scaffold)
        x_all = np.concatenate(x_all)
        y_all = np.concatenate(y_all)
        groups = np.array(groups)

        # Property: depends on position + scaffold
        prop = 0.5 * x_all + 0.3 * y_all + groups * 0.4 + np.random.normal(0, 0.3, n)

        test_frac = input.test_frac_split()

        if input.split_type() == "random":
            idx = np.random.permutation(n)
            test_idx = idx[:int(n * test_frac)]
            train_idx = idx[int(n * test_frac):]
            title = "Random Split"
        else:
            # Scaffold split: entire scaffold groups go to test
            n_test_scaffolds = max(1, int(n_scaffolds * test_frac))
            test_scaffolds = set(range(n_test_scaffolds))
            test_idx = np.where(np.isin(groups, list(test_scaffolds)))[0]
            train_idx = np.where(~np.isin(groups, list(test_scaffolds)))[0]
            title = "Scaffold-Aware Split"

        cmap = plt.cm.Set2
        colors = [cmap(g / max(n_scaffolds - 1, 1)) for g in groups]

        fig, axes = plt.subplots(1, 2, figsize=(6, 4))

        for ax, idx_set, label in [(axes[0], train_idx, "Train"), (axes[1], test_idx, "Test")]:
            ax.set_facecolor("#fafaf7")
            for i in idx_set:
                ax.scatter(x_all[i], y_all[i], c=[colors[i]], s=35, alpha=0.8, edgecolors="white", linewidth=0.5, zorder=3)
            ax.set_title(f"{label} set ({len(idx_set)} pts)")
            ax.set_xlabel("Feature 1")
            ax.set_ylabel("Feature 2")
            ax.set_xlim(-5, 5)
            ax.set_ylim(-5, 5)

        fig.suptitle(title, fontsize=13, fontweight="bold")
        fig.tight_layout()
        return fig

app = App(app_ui, server)

With random splitting, each scaffold group has members in both train and test — the model can interpolate within known scaffolds. With scaffold splitting, entire groups are held out, forcing the model to extrapolate to genuinely new chemical series. Neither number is the “true” performance: they answer different questions, and the difference between them tells you how much of the random-split score rested on scaffolds the model had already seen.

4.6 Applicability Domains

Splitting strategies address the question at evaluation time: they let you estimate, before deployment, how the model behaves under a distribution shift you have anticipated. The applicability domain is the same question asked at prediction time, one molecule at a time — given this specific query, is the model being asked to interpolate or to extrapolate?

A model should only be trusted for predictions on molecules that are “similar enough” to the training data. The applicability domain (AD) defines the region of chemical space where the model’s predictions are reliable.

Outside the applicability domain, the model is extrapolating rather than interpolating, and predictions become unreliable. This matters most when the cost of a wrong prediction is high — a model predicting the potential toxicity of a compound must not return a confident “non-toxic” for chemistry it has never seen.

4.6.1 Common approaches to define an AD

Bounding box: The simplest approach — check whether each feature falls within the range observed in training. Crude but fast.
Distance to training data: Compute the distance (Euclidean, Tanimoto, etc.) from a query molecule to its nearest neighbor(s) in the training set. Flag if the distance exceeds a threshold.
Leverage approach: A point’s leverage measures how far it lies from the centre of the training distribution in feature space, and hence how much influence it would have on a linear fit. It is read off the diagonal of the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$, so called because it maps observed values onto fitted ones, $\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}$.
Density-based: Estimate the density of the training data in feature space. Low-density regions are outside the AD.

4.6.2 Interactive example: applicability domain

Explore how distance-based applicability domain detection works. Points far from the training data (beyond the threshold) are flagged as outside the AD, and their predictions should be treated with caution.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 560

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.patches import Circle
matplotlib.rcParams.update({
    "font.family": "serif",
    "figure.facecolor": "#fafaf7",
    "axes.facecolor": "#fafaf7",
    "axes.grid": True,
    "grid.alpha": 0.3,
})

app_ui = ui.page_fluid(
    ui.layout_columns(
        ui.input_slider("ad_threshold", "AD distance threshold", min=0.5, max=4.0, value=2.0, step=0.25),
        ui.input_slider("n_query", "Number of query points", min=5, max=30, value=15, step=1),
        ui.input_radio_buttons("ad_method", "AD method",
            choices={"knn": "k-NN distance", "centroid": "Centroid distance"},
            selected="knn"),
    ),
    ui.output_plot("ad_plot"),
)

def server(input, output, session):
    @reactive.calc
    def compute_ad():
        np.random.seed(55)
        # Training data: two clusters
        n_train = 50
        cluster1 = np.random.normal([1, 1], 0.8, (25, 2))
        cluster2 = np.random.normal([-1.5, -1], 0.8, (25, 2))
        x_train = np.vstack([cluster1, cluster2])

        # Query points: spread more widely
        n_q = input.n_query()
        x_query = np.random.uniform(-5, 5, (n_q, 2))

        threshold = input.ad_threshold()

        if input.ad_method() == "knn":
            # k-NN distance (k=3)
            k = 3
            distances = []
            for q in x_query:
                dists = np.sqrt(np.sum((x_train - q)**2, axis=1))
                distances.append(np.mean(np.sort(dists)[:k]))
            distances = np.array(distances)
        else:
            centroid = np.mean(x_train, axis=0)
            distances = np.sqrt(np.sum((x_query - centroid)**2, axis=1))

        in_ad = distances <= threshold
        return x_train, x_query, distances, in_ad, threshold

    @output
    @render.plot
    def ad_plot():
        x_train, x_query, distances, in_ad, threshold = compute_ad()

        fig, ax = plt.subplots(figsize=(6, 4.5))
        ax.scatter(x_train[:, 0], x_train[:, 1], c="#2b6cb0", s=25, alpha=0.5, label="Training data", zorder=3)

        # Query points
        for i in range(len(x_query)):
            color = "#276749" if in_ad[i] else "#c53030"
            marker = "o" if in_ad[i] else "x"
            ax.scatter(x_query[i, 0], x_query[i, 1], c=color, s=60, marker=marker,
                      linewidths=2, zorder=4)

        # Draw AD boundary approximation
        if input.ad_method() == "centroid":
            centroid = np.mean(x_train, axis=0)
            circle = Circle(centroid, threshold, fill=False, linestyle="--",
                          edgecolor="#e8833a", linewidth=1.5, alpha=0.7)
            ax.add_patch(circle)
            ax.scatter([centroid[0]], [centroid[1]], c="#e8833a", s=100, marker="+", linewidths=2, zorder=5)

        ax.set_xlabel("Feature 1")
        ax.set_ylabel("Feature 2")
        ax.set_title("Applicability Domain")
        ax.set_xlim(-5.5, 5.5)
        ax.set_ylim(-5.5, 5.5)
        ax.set_aspect("equal")

        # Custom legend
        from matplotlib.lines import Line2D
        legend_elements = [
            Line2D([0], [0], marker="o", color="w", markerfacecolor="#2b6cb0", markersize=8, label="Training"),
            Line2D([0], [0], marker="o", color="w", markerfacecolor="#276749", markersize=8, label="Query (in AD)"),
            Line2D([0], [0], marker="x", color="#c53030", markersize=8, markeredgewidth=2, linestyle="None", label="Query (outside AD)"),
        ]
        ax.legend(handles=legend_elements, fontsize=9, loc="upper right")
        fig.tight_layout()
        return fig

app = App(app_ui, server)

Every prediction should carry an uncertainty estimate, such as a confidence interval from an ensemble of models or a Gaussian process, whether or not the query falls inside the domain. The applicability domain answers a different question: not how uncertain a prediction is, but whether the model is being asked to extrapolate at all. That is what determines the action taken:

Reject: Do not provide a prediction at all.
Flag: Provide the prediction but mark it as outside the domain the model was trained for.

AD in regulatory settings

Reporting the applicability domain is increasingly expected in regulatory submissions (e.g., REACH in the EU) when QSAR models are used to fill data gaps for chemical safety assessments — which covers industrial chemicals, agrochemicals and cosmetics ingredients, not only pharmaceuticals. The OECD guidelines for QSAR validation explicitly require AD characterization.

4.7 Summary

A performance score only means something once you say what the model was kept away from. The test set stays untouched until the end, and when data is scarce, cross-validation reuses every point while also giving the uncertainty on the estimate. The split then decides what the score measures, since random, cluster, scaffold and temporal splits each withhold a different kind of chemistry, and a number is only interpretable together with the split that produced it. The metric has to match the task, and overfitting shows up in the held-out score as complexity grows rather than in the gap between training and validation. All of this describes average behaviour, which leaves the applicability domain to answer the remaining question of whether the model is interpolating or extrapolating for one particular molecule.

--- title: "Model Evaluation" format: html: toc: true theme: cosmo code-fold: true filters: - shinylive --- A model is only as useful as its ability to generalize to new, unseen data. This chapter covers the key principles of evaluating machine learning models in the context of chemistry, where naive evaluation strategies can paint a dangerously optimistic picture. ## Train / Validation / Test Philosophy The fundamental principle: **never evaluate a model on the data it was trained on.** To get an honest estimate of how well a model will perform in practice, we split our data into distinct subsets: - **Training set** (~60--70%): Used to fit the model parameters. - **Validation set** (~15--20%): Used to tune hyperparameters (e.g., regularization strength, polynomial degree) and select among competing models. - **Test set** (~15--20%): Held out until the very end. Touched exactly once to report final, unbiased performance. The validation set acts as a proxy for unseen data *during* model development. If we used the test set to choose hyperparameters, we would effectively be "training" on it and our final performance estimate would be optimistic. A common mistake is to repeatedly peek at the test set during development. Even if the test set is never directly used for gradient updates, selecting the model that performs best on it causes **information leakage** --- the test set silently guides our modeling choices. ### Interactive example: train/validation/test split Adjust the split ratios and the noise level to see how the amount of training data and the split proportions affect model performance across the three subsets. Only training and validation are set directly; the test set takes whatever remains, and the plot title reports the resulting three-way split. ```{shinylive-python} #| standalone: true #| viewerHeight: 600 from shiny import App, render, ui, reactive import numpy as np import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({ "font.family": "serif", "figure.facecolor": "#fafaf7", "axes.facecolor": "#fafaf7", "axes.grid": True, "grid.alpha": 0.3, }) app_ui = ui.page_fluid( ui.layout_columns( ui.input_slider("train_frac", "Training fraction", min=0.3, max=0.8, value=0.6, step=0.05), ui.input_slider("val_frac", "Validation fraction", min=0.05, max=0.65, value=0.2, step=0.05), ui.input_slider("noise", "Noise level", min=0.1, max=2.0, value=0.8, step=0.1), ui.input_slider("degree", "Polynomial degree", min=1, max=12, value=3, step=1), ui.input_slider("n_points", "Total data points", min=30, max=200, value=80, step=10), col_widths=(4, 4, 4, 6, 6), ), ui.output_plot("split_plot"), ) def server(input, output, session): @reactive.effect @reactive.event(input.train_frac) def _limit_val_frac(): # Cap the validation slider at what the training fraction leaves behind, # keeping at least 5% for the test set. Without this the slider can be # set to a value the plot cannot honour. cap = round(1.0 - input.train_frac() - 0.05, 2) ui.update_slider("val_frac", max=cap, value=min(input.val_frac(), cap)) @output @render.plot def split_plot(): np.random.seed(42) n = input.n_points() x = np.sort(np.random.uniform(-3, 3, n)) y_true = np.sin(x) + 0.3 * x y = y_true + np.random.normal(0, input.noise(), n) # The test set is whatever is left over, so the validation fraction is # capped to always leave at least 5% of the points for testing. Without # this the test set can be empty and its score comes out as nan. train_frac = input.train_frac() val_frac = min(input.val_frac(), 1.0 - train_frac - 0.05) test_frac = 1.0 - train_frac - val_frac train_end = int(n * train_frac) val_end = train_end + int(n * val_frac) idx = np.random.permutation(n) train_idx = idx[:train_end] val_idx = idx[train_end:val_end] test_idx = idx[val_end:] x_train, y_train = x[train_idx], y[train_idx] x_val, y_val = x[val_idx], y[val_idx] x_test, y_test = x[test_idx], y[test_idx] deg = input.degree() coeffs = np.polyfit(x_train, y_train, deg) poly = np.poly1d(coeffs) x_line = np.linspace(-3.2, 3.2, 300) y_pred_line = poly(x_line) mse_train = np.mean((y_train - poly(x_train))**2) mse_val = np.mean((y_val - poly(x_val))**2) if len(x_val) > 0 else float("nan") mse_test = np.mean((y_test - poly(x_test))**2) if len(x_test) > 0 else float("nan") fig, ax = plt.subplots(figsize=(6, 4)) ax.scatter(x_train, y_train, c="#2b6cb0", s=30, alpha=0.7, label=f"Train (MSE={mse_train:.3f})", zorder=3) ax.scatter(x_val, y_val, c="#e8833a", s=30, alpha=0.7, label=f"Val (MSE={mse_val:.3f})", zorder=3) ax.scatter(x_test, y_test, c="#c53030", s=30, alpha=0.7, label=f"Test (MSE={mse_test:.3f})", zorder=3) ax.plot(x_line, np.sin(x_line) + 0.3 * x_line, "--", color="#999", linewidth=1.2, label="True function") y_clipped = np.clip(y_pred_line, -6, 6) ax.plot(x_line, y_clipped, color="#276749", linewidth=2, label=f"Degree {deg} fit") ax.set_ylim(-5, 5) ax.set_xlabel("x") ax.set_ylabel("y") ax.set_title(f"Train {train_frac:.0%} / Val {val_frac:.0%} / Test {test_frac:.0%}") ax.legend(fontsize=9, loc="upper left") fig.tight_layout() return fig app = App(app_ui, server) ``` Notice how the test MSE gives a more honest picture of generalization than the training MSE, especially for high-degree polynomials where the model memorizes the training data. ## K-Fold Cross-Validation When data is limited (a common situation in chemistry), setting aside a fixed validation set wastes precious training examples. **K-fold cross-validation** offers a solution: partition the data into $k$ equally sized folds, then train $k$ separate models, each time holding out one fold as the validation set and training on the remaining $k - 1$ folds. The cross-validated score is the average performance over the $k$ held-out folds: $$ \text{CV}(k) = \frac{1}{k} \sum_{i=1}^{k} \mathcal{L}_i $$ where $\mathcal{L}_i$ is the loss on fold $i$. Common choices are $k = 5$ or $k = 10$. The extreme case $k = n$ is **leave-one-out cross-validation (LOOCV)**, which has low bias but high variance and is computationally expensive. (The bias--variance trade-off is introduced in the [linear models chapter](03_linear_models.qmd).) **Advantages:** - Every data point is used for both training and validation. - Yields $k$ scores rather than a single number, so the result comes with an uncertainty. We interpret the **standard error of the mean** , $\sigma/\sqrt{k}$, where $\sigma$ is the standard deviation across the $k$ folds, to give a sense of the variability in the cross-validated performance. **Disadvantages:** - $k$ times more expensive than a single train/val split. - Folds can be correlated, especially for LOOCV. ## Metrics Choosing the right metric is essential --- a model that looks excellent by one metric may be mediocre by another. The appropriate metric depends on whether the task is **regression** or **classification**. ### Regression metrics **Root Mean Squared Error (RMSE):** Penalizes large errors more heavily due to squaring. $$ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} $$ **Mean Absolute Error (MAE):** More robust to outliers; gives equal weight to all errors. $$ \text{MAE} = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i| $$ **Coefficient of Determination ($R^2$):** Measures the proportion of variance explained. $R^2 = 1$ is perfect; $R^2 = 0$ means the model is no better than predicting the mean; $R^2 < 0$ means the model is *worse* than the mean. $$ R^2 = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2} $$ ### Interactive example: regression metrics Drag the noise slider and add an outlier to see how RMSE, MAE, and $R^2$ respond differently. ```{shinylive-python} #| standalone: true #| viewerHeight: 520 from shiny import App, render, ui, reactive import numpy as np import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({ "font.family": "serif", "figure.facecolor": "#fafaf7", "axes.facecolor": "#fafaf7", "axes.grid": True, "grid.alpha": 0.3, }) app_ui = ui.page_fluid( ui.layout_columns( ui.input_slider("noise_met", "Noise level", min=0.1, max=3.0, value=0.5, step=0.1), ui.input_slider("outlier_size", "Outlier magnitude", min=0.0, max=8.0, value=0.0, step=0.5), ui.input_slider("n_pts_met", "Number of points", min=20, max=100, value=40, step=5), ), ui.output_plot("metrics_plot"), ) def server(input, output, session): @reactive.calc def compute_data(): np.random.seed(21) n = input.n_pts_met() x = np.sort(np.random.uniform(0, 10, n)) y_true = 0.8 * x + 1.0 y = y_true + np.random.normal(0, input.noise_met(), n) # Add outlier at the last point if input.outlier_size() > 0: y[-1] = y_true[-1] + input.outlier_size() # Simple linear fit coeffs = np.polyfit(x, y, 1) y_pred = np.polyval(coeffs, x) rmse = np.sqrt(np.mean((y - y_pred)**2)) mae = np.mean(np.abs(y - y_pred)) ss_res = np.sum((y - y_pred)**2) ss_tot = np.sum((y - np.mean(y))**2) r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0 return x, y, y_true, y_pred, coeffs, rmse, mae, r2 @output @render.plot def metrics_plot(): x, y, y_true, y_pred, coeffs, rmse, mae, r2 = compute_data() fig, ax = plt.subplots(figsize=(6, 4)) ax.scatter(x, y, c="#2b6cb0", s=30, alpha=0.7, zorder=3, label="Data") x_line = np.linspace(-0.5, 10.5, 100) ax.plot(x_line, np.polyval(coeffs, x_line), color="#e8833a", linewidth=2, label="Linear fit") ax.plot(x_line, 0.8 * x_line + 1.0, "--", color="#999", linewidth=1, label="True line") # Draw residuals for xi, yi, yp in zip(x, y, y_pred): ax.plot([xi, xi], [yi, yp], color="#c53030", alpha=0.3, linewidth=1) if input.outlier_size() > 0: ax.scatter([x[-1]], [y[-1]], c="#c53030", s=80, zorder=5, marker="x", linewidths=2, label="Outlier") ax.set_xlabel("x") ax.set_ylabel("y") ax.set_title("Regression Metrics Comparison") ax.legend(fontsize=9, loc="upper left") ax.text(0.98, 0.02, f"RMSE = {rmse:.3f}\nMAE = {mae:.3f}\nR² = {r2:.3f}", transform=ax.transAxes, ha="right", va="bottom", fontsize=10, family="monospace", bbox=dict(boxstyle="round,pad=0.4", facecolor="white", edgecolor="#ddd")) fig.tight_layout() return fig app = App(app_ui, server) ``` Note that RMSE is always greater than or equal to MAE, and a single large outlier inflates RMSE far more than MAE. ### Classification metrics For binary classification tasks (e.g., active/inactive, toxic/non-toxic), we use different metrics. **Confusion matrix** terminology: | | Predicted positive | Predicted negative | |---|---|---| | **Actual positive** | True Positive (TP) | False Negative (FN) | | **Actual negative** | False Positive (FP) | True Negative (TN) | **Accuracy** $= \frac{TP + TN}{TP + TN + FP + FN}$ --- misleading for imbalanced datasets. **Precision** $= \frac{TP}{TP + FP}$ --- of those predicted positive, how many are correct? **Recall (Sensitivity)** $= \frac{TP}{TP + FN}$ --- of actual positives, how many did we catch? **ROC-AUC:** The area under the Receiver Operating Characteristic curve, which plots recall vs. false positive rate at various classification thresholds. AUC = 1.0 is perfect; AUC = 0.5 is random guessing. ### Interactive example: ROC curve and threshold Adjust the classification threshold and class separability to see how the ROC curve, precision, and recall change. ```{shinylive-python} #| standalone: true #| viewerHeight: 700 from shiny import App, render, ui, reactive import numpy as np import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({ "font.family": "serif", "figure.facecolor": "#fafaf7", "axes.facecolor": "#fafaf7", "axes.grid": True, "grid.alpha": 0.3, }) app_ui = ui.page_fluid( ui.layout_columns( ui.input_slider("separability", "Class separability", min=0.5, max=4.0, value=2.0, step=0.25), ui.input_slider("threshold", "Classification threshold", min=0.0, max=1.0, value=0.5, step=0.02), ui.input_slider("imbalance", "Class imbalance (% positive)", min=10, max=90, value=50, step=5), ), ui.output_plot("roc_plot", height="560px"), ) def server(input, output, session): @reactive.calc def compute_class_data(): np.random.seed(33) n = 200 n_pos = int(n * input.imbalance() / 100) n_neg = n - n_pos sep = input.separability() scores_pos = np.random.normal(sep / 2, 1, n_pos) scores_neg = np.random.normal(-sep / 2, 1, n_neg) # Convert to probabilities via sigmoid def sigmoid(x): return 1 / (1 + np.exp(-x)) probs_pos = sigmoid(scores_pos) probs_neg = sigmoid(scores_neg) y_true = np.concatenate([np.ones(n_pos), np.zeros(n_neg)]) y_probs = np.concatenate([probs_pos, probs_neg]) return y_true, y_probs @output @render.plot def roc_plot(): from matplotlib.gridspec import GridSpec from matplotlib.colors import ListedColormap y_true, y_probs = compute_class_data() thresh = input.threshold() # ROC curve over all thresholds thresholds = np.linspace(0, 1, 200) tpr_list = [] fpr_list = [] for t in thresholds: y_pred = (y_probs >= t).astype(int) tp_t = np.sum((y_pred == 1) & (y_true == 1)) fn_t = np.sum((y_pred == 0) & (y_true == 1)) fp_t = np.sum((y_pred == 1) & (y_true == 0)) tn_t = np.sum((y_pred == 0) & (y_true == 0)) tpr_list.append(tp_t / max(tp_t + fn_t, 1)) fpr_list.append(fp_t / max(fp_t + tn_t, 1)) sorted_idx = np.argsort(fpr_list) fpr_sorted = np.array(fpr_list)[sorted_idx] tpr_sorted = np.array(tpr_list)[sorted_idx] auc = np.trapz(tpr_sorted, fpr_sorted) # Confusion matrix at the current threshold y_pred_curr = (y_probs >= thresh).astype(int) tp = int(np.sum((y_pred_curr == 1) & (y_true == 1))) fn = int(np.sum((y_pred_curr == 0) & (y_true == 1))) fp = int(np.sum((y_pred_curr == 1) & (y_true == 0))) tn = int(np.sum((y_pred_curr == 0) & (y_true == 0))) curr_tpr = tp / max(tp + fn, 1) curr_fpr = fp / max(fp + tn, 1) accuracy = (tp + tn) / len(y_true) precision = tp / max(tp + fp, 1) recall = tp / max(tp + fn, 1) fig = plt.figure(figsize=(6, 6.2)) gs = GridSpec(2, 2, figure=fig, height_ratios=[1, 1], hspace=0.55, wspace=0.4) ax = fig.add_subplot(gs[0, 0]) ax2 = fig.add_subplot(gs[0, 1]) ax_cm = fig.add_subplot(gs[1, 0]) ax_met = fig.add_subplot(gs[1, 1]) # Top-left: score distributions pos_probs = y_probs[y_true == 1] neg_probs = y_probs[y_true == 0] ax.hist(neg_probs, bins=25, alpha=0.6, color="#2b6cb0", label="Negative", density=True) ax.hist(pos_probs, bins=25, alpha=0.6, color="#c53030", label="Positive", density=True) ax.axvline(thresh, color="#e8833a", linewidth=2, linestyle="--", label=f"Threshold={thresh:.2f}") ax.set_xlabel("Predicted probability") ax.set_ylabel("Density") ax.set_title("Score distributions") ax.legend(fontsize=8) # Top-right: ROC curve ax2.plot(fpr_sorted, tpr_sorted, color="#2b6cb0", linewidth=2, label=f"ROC (AUC={auc:.3f})") ax2.plot([0, 1], [0, 1], "--", color="#999", linewidth=1, label="Random") ax2.scatter([curr_fpr], [curr_tpr], c="#e8833a", s=80, zorder=5, label="Current threshold") ax2.set_xlabel("False Positive Rate") ax2.set_ylabel("True Positive Rate") ax2.set_title("ROC curve") ax2.legend(fontsize=8, loc="lower right") ax2.set_xlim(-0.02, 1.02) ax2.set_ylim(-0.02, 1.02) # Bottom-left: confusion matrix at the current threshold cm_cmap = ListedColormap(["#f5c6c6", "#c6e6c6"]) grid = np.array([[1, 0], [0, 1]]) ax_cm.imshow(grid, cmap=cm_cmap, vmin=0, vmax=1) ax_cm.grid(False) cells = [[("TP", tp), ("FN", fn)], [("FP", fp), ("TN", tn)]] for r in range(2): for c in range(2): lbl, val = cells[r][c] ax_cm.text(c, r, f"{lbl}\n{val}", ha="center", va="center", fontsize=12, fontweight="bold", color="#222") ax_cm.set_xticks([0, 1]) ax_cm.set_xticklabels(["Pred +", "Pred -"]) ax_cm.set_yticks([0, 1]) ax_cm.set_yticklabels(["Actual +", "Actual -"]) ax_cm.set_title("Confusion matrix") # Bottom-right: derived metrics ax_met.axis("off") ax_met.text(0.0, 0.9, f"Accuracy = {accuracy:.3f}\n\n" f"Precision = {precision:.3f}\n\n" f"Recall = {recall:.3f}", transform=ax_met.transAxes, va="top", ha="left", fontsize=12, family="monospace") return fig app = App(app_ui, server) ``` ::: {.callout-tip} ## Which metric to choose? - Use **RMSE** when large errors are particularly costly (e.g., toxicity prediction). - Use **MAE** when you want robustness to outliers (e.g., noisy experimental data). - Use **$R^2$** for a quick sense of overall fit quality, but note that it is relative to the variance of *your* test set: the same model scores higher on a chemically diverse set than on a narrow one, so $R^2$ values are not comparable across datasets. - Use **ROC-AUC** for balanced classification. For imbalanced tasks, consider **Precision-Recall AUC** instead. ::: ## Overfitting and Underfitting A metric tells you how well the model scored; it does not tell you *why* it scored that way. The two failure modes behind a poor score are opposites, and they call for opposite remedies. **Underfitting** (high bias): The model is too simple to capture the underlying pattern, and performs poorly on held-out data because it never captured the signal in the first place. **Overfitting** (high variance): The model has absorbed noise specific to the training set, and performs poorly on held-out data because that noise does not recur. Both are diagnosed from the **held-out score**, and distinguished by how that score responds to model complexity: it improves with added complexity if you are underfitting, and degrades if you are overfitting. ::: {.callout-warning} ## The training-error gap is not a reliable diagnostic It is tempting to diagnose overfitting from the *gap* between training and validation error, but for many models a low training error is structural rather than symptomatic: a random forest of fully grown trees drives training error to nearly zero by construction, without necessarily harming the held-out score. Judge the model by the held-out score and how it responds to complexity, not by the gap. ::: ### Remedies | Problem | Remedies | |---|---| | **Underfitting** (held-out score poor, improves with complexity) | Increase model complexity, add more features, reduce regularization, train longer | | **Overfitting** (held-out score degrades as complexity grows) | Add regularization (L1/L2), reduce model complexity, get more data, use dropout, early stopping | ::: {.callout-note} ## Early stopping For iterative models (neural networks, gradient boosting), a practical remedy for overfitting is **early stopping**: monitor the validation loss during training and stop when it begins to increase, even if the training loss is still decreasing. ::: ## How You Split Decides What You Measure Everything so far has assumed that the test set tells you how the model will perform in use. That is only true if the test set resembles what the model will actually see. In most machine learning applications a random split makes this approximately true, because the data are independent and identically distributed (i.i.d.). In chemistry the assumption frequently breaks down: molecules in a dataset are rarely independent draws. They share scaffolds and functional groups, arrive in congeneric series (families of close analogues built around a common core), and are the product of a project that deliberately explored one region of chemical space at a time. A random split therefore leaves near-duplicates on both sides, and the model can score well by recognising a scaffold it has already seen. But the deeper problem is not that similar molecules appear in both train and test: it is that **the molecules you will actually predict on are dissimilar to both**. What the test score fails to capture is **out-of-distribution (OOD) generalization**. ### Splitting strategies as probes of OOD generalization Each split withholds a different kind of chemistry from the model, and so simulates a different distribution shift. The right one is whichever shift you expect between the data you have and the molecules you will actually predict on. | Strategy | What the model never sees in training | Choose it when | |---|---|---| | **Random** | Nothing in particular; the test molecules resemble the training ones | The model will be used on the same chemical space the data came from, so the task really is interpolation | | **Cluster** | Entire regions of feature space | The model will meet broadly new chemistry; cluster the fingerprints (k-means, Butina) and hold out whole clusters | | **Scaffold** | The molecular core | The model will see new series rather than new analogues; Murcko scaffolds are the usual heuristic for "same core", not the only one | | **Temporal** | Whatever changed from one period to the next: reformulated assays, a move to a new chemotype, easy targets done first | The model will run on data generated *after* the training data (not merely "the future", which is true of any prediction) | Two consequences: a reported score is only interpretable together with its split, and the aggressive splits are not simply better -- a scaffold split gives a pessimistic number for a task that really is interpolation. ### Interactive example: random vs. scaffold split This demonstration shows how random splits allow structurally similar molecules to leak across train and test sets, while scaffold-aware splits keep related molecules together. Points represent molecules; color represents their scaffold group. ```{shinylive-python} #| standalone: true #| viewerHeight: 560 from shiny import App, render, ui, reactive import numpy as np import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams.update({ "font.family": "serif", "figure.facecolor": "#fafaf7", "axes.facecolor": "#fafaf7", "axes.grid": True, "grid.alpha": 0.3, }) app_ui = ui.page_fluid( ui.layout_columns( ui.input_radio_buttons("split_type", "Split strategy", choices={"random": "Random", "scaffold": "Scaffold-aware"}, selected="random"), ui.input_slider("n_scaffolds", "Number of scaffold groups", min=3, max=8, value=5, step=1), ui.input_slider("test_frac_split", "Test fraction", min=0.1, max=0.4, value=0.2, step=0.05), ), ui.output_plot("split_compare_plot"), ) def server(input, output, session): @output @render.plot def split_compare_plot(): np.random.seed(12) n_scaffolds = input.n_scaffolds() pts_per_scaffold = 15 n = n_scaffolds * pts_per_scaffold # Generate clustered 2D data (simulating scaffold groups) centers = np.random.uniform(-3, 3, (n_scaffolds, 2)) x_all = [] y_all = [] groups = [] for i, c in enumerate(centers): pts = np.random.normal(c, 0.5, (pts_per_scaffold, 2)) x_all.append(pts[:, 0]) y_all.append(pts[:, 1]) groups.extend([i] * pts_per_scaffold) x_all = np.concatenate(x_all) y_all = np.concatenate(y_all) groups = np.array(groups) # Property: depends on position + scaffold prop = 0.5 * x_all + 0.3 * y_all + groups * 0.4 + np.random.normal(0, 0.3, n) test_frac = input.test_frac_split() if input.split_type() == "random": idx = np.random.permutation(n) test_idx = idx[:int(n * test_frac)] train_idx = idx[int(n * test_frac):] title = "Random Split" else: # Scaffold split: entire scaffold groups go to test n_test_scaffolds = max(1, int(n_scaffolds * test_frac)) test_scaffolds = set(range(n_test_scaffolds)) test_idx = np.where(np.isin(groups, list(test_scaffolds)))[0] train_idx = np.where(~np.isin(groups, list(test_scaffolds)))[0] title = "Scaffold-Aware Split" cmap = plt.cm.Set2 colors = [cmap(g / max(n_scaffolds - 1, 1)) for g in groups] fig, axes = plt.subplots(1, 2, figsize=(6, 4)) for ax, idx_set, label in [(axes[0], train_idx, "Train"), (axes[1], test_idx, "Test")]: ax.set_facecolor("#fafaf7") for i in idx_set: ax.scatter(x_all[i], y_all[i], c=[colors[i]], s=35, alpha=0.8, edgecolors="white", linewidth=0.5, zorder=3) ax.set_title(f"{label} set ({len(idx_set)} pts)") ax.set_xlabel("Feature 1") ax.set_ylabel("Feature 2") ax.set_xlim(-5, 5) ax.set_ylim(-5, 5) fig.suptitle(title, fontsize=13, fontweight="bold") fig.tight_layout() return fig app = App(app_ui, server) ``` With random splitting, each scaffold group has members in both train and test --- the model can interpolate within known scaffolds. With scaffold splitting, entire groups are held out, forcing the model to extrapolate to genuinely new chemical series. Neither number is the "true" performance: they answer different questions, and the difference between them tells you how much of the random-split score rested on scaffolds the model had already seen. ## Applicability Domains Splitting strategies address the question at *evaluation* time: they let you estimate, before deployment, how the model behaves under a distribution shift you have anticipated. The applicability domain is the same question asked at *prediction* time, one molecule at a time — given this specific query, is the model being asked to interpolate or to extrapolate? A model should only be trusted for predictions on molecules that are "similar enough" to the training data. The **applicability domain (AD)** defines the region of chemical space where the model's predictions are reliable. Outside the applicability domain, the model is **extrapolating** rather than interpolating, and predictions become unreliable. This matters most when the cost of a wrong prediction is high --- a model predicting the potential toxicity of a compound must not return a confident "non-toxic" for chemistry it has never seen. ### Common approaches to define an AD 1. **Bounding box:** The simplest approach --- check whether each feature falls within the range observed in training. Crude but fast. 2. **Distance to training data:** Compute the distance (Euclidean, Tanimoto, etc.) from a query molecule to its nearest neighbor(s) in the training set. Flag if the distance exceeds a threshold. 3. **Leverage approach:** A point's **leverage** measures how far it lies from the centre of the training distribution in feature space, and hence how much influence it would have on a linear fit. It is read off the diagonal of the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$, so called because it maps observed values onto fitted ones, $\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}$. 4. **Density-based:** Estimate the density of the training data in feature space. Low-density regions are outside the AD. ### Interactive example: applicability domain Explore how distance-based applicability domain detection works. Points far from the training data (beyond the threshold) are flagged as outside the AD, and their predictions should be treated with caution. ```{shinylive-python} #| standalone: true #| viewerHeight: 560 from shiny import App, render, ui, reactive import numpy as np import matplotlib.pyplot as plt import matplotlib from matplotlib.patches import Circle matplotlib.rcParams.update({ "font.family": "serif", "figure.facecolor": "#fafaf7", "axes.facecolor": "#fafaf7", "axes.grid": True, "grid.alpha": 0.3, }) app_ui = ui.page_fluid( ui.layout_columns( ui.input_slider("ad_threshold", "AD distance threshold", min=0.5, max=4.0, value=2.0, step=0.25), ui.input_slider("n_query", "Number of query points", min=5, max=30, value=15, step=1), ui.input_radio_buttons("ad_method", "AD method", choices={"knn": "k-NN distance", "centroid": "Centroid distance"}, selected="knn"), ), ui.output_plot("ad_plot"), ) def server(input, output, session): @reactive.calc def compute_ad(): np.random.seed(55) # Training data: two clusters n_train = 50 cluster1 = np.random.normal([1, 1], 0.8, (25, 2)) cluster2 = np.random.normal([-1.5, -1], 0.8, (25, 2)) x_train = np.vstack([cluster1, cluster2]) # Query points: spread more widely n_q = input.n_query() x_query = np.random.uniform(-5, 5, (n_q, 2)) threshold = input.ad_threshold() if input.ad_method() == "knn": # k-NN distance (k=3) k = 3 distances = [] for q in x_query: dists = np.sqrt(np.sum((x_train - q)**2, axis=1)) distances.append(np.mean(np.sort(dists)[:k])) distances = np.array(distances) else: centroid = np.mean(x_train, axis=0) distances = np.sqrt(np.sum((x_query - centroid)**2, axis=1)) in_ad = distances <= threshold return x_train, x_query, distances, in_ad, threshold @output @render.plot def ad_plot(): x_train, x_query, distances, in_ad, threshold = compute_ad() fig, ax = plt.subplots(figsize=(6, 4.5)) ax.scatter(x_train[:, 0], x_train[:, 1], c="#2b6cb0", s=25, alpha=0.5, label="Training data", zorder=3) # Query points for i in range(len(x_query)): color = "#276749" if in_ad[i] else "#c53030" marker = "o" if in_ad[i] else "x" ax.scatter(x_query[i, 0], x_query[i, 1], c=color, s=60, marker=marker, linewidths=2, zorder=4) # Draw AD boundary approximation if input.ad_method() == "centroid": centroid = np.mean(x_train, axis=0) circle = Circle(centroid, threshold, fill=False, linestyle="--", edgecolor="#e8833a", linewidth=1.5, alpha=0.7) ax.add_patch(circle) ax.scatter([centroid[0]], [centroid[1]], c="#e8833a", s=100, marker="+", linewidths=2, zorder=5) ax.set_xlabel("Feature 1") ax.set_ylabel("Feature 2") ax.set_title("Applicability Domain") ax.set_xlim(-5.5, 5.5) ax.set_ylim(-5.5, 5.5) ax.set_aspect("equal") # Custom legend from matplotlib.lines import Line2D legend_elements = [ Line2D([0], [0], marker="o", color="w", markerfacecolor="#2b6cb0", markersize=8, label="Training"), Line2D([0], [0], marker="o", color="w", markerfacecolor="#276749", markersize=8, label="Query (in AD)"), Line2D([0], [0], marker="x", color="#c53030", markersize=8, markeredgewidth=2, linestyle="None", label="Query (outside AD)"), ] ax.legend(handles=legend_elements, fontsize=9, loc="upper right") fig.tight_layout() return fig app = App(app_ui, server) ``` Every prediction should carry an **uncertainty estimate**, such as a confidence interval from an ensemble of models or a Gaussian process, whether or not the query falls inside the domain. The applicability domain answers a different question: not how uncertain a prediction is, but whether the model is being asked to extrapolate at all. That is what determines the action taken: - **Reject:** Do not provide a prediction at all. - **Flag:** Provide the prediction but mark it as outside the domain the model was trained for. ::: {.callout-important} ## AD in regulatory settings Reporting the applicability domain is increasingly expected in regulatory submissions (e.g., REACH in the EU) when QSAR models are used to fill data gaps for chemical safety assessments --- which covers industrial chemicals, agrochemicals and cosmetics ingredients, not only pharmaceuticals. The OECD guidelines for QSAR validation explicitly require AD characterization. ::: ## Summary A performance score only means something once you say what the model was kept away from. The test set stays untouched until the end, and when data is scarce, cross-validation reuses every point while also giving the uncertainty on the estimate. The split then decides what the score measures, since random, cluster, scaffold and temporal splits each withhold a different kind of chemistry, and a number is only interpretable together with the split that produced it. The metric has to match the task, and overfitting shows up in the held-out score as complexity grows rather than in the gap between training and validation. All of this describes average behaviour, which leaves the applicability domain to answer the remaining question of whether the model is interpolating or extrapolating for one particular molecule.