Validation

Model validation against OECD principle 4, which asks for three distinct things — goodness-of-fit, robustness and predictivity. A single \(R^2\) addresses only the first, and it is the number most likely to be quoted.

Class

Answers

CrossValidator

Does it predict? Cross-validated \(Q^2\), with out-of-fold predictions.

YScrambling

Is the fit real? Could the model score this well on permuted labels?

BootstrapValidator

How precise is the score? An out-of-bag interval rather than a point estimate.

ExternalValidator

Does it predict on compounds it never saw? The full QSAR metric set plus the Golbraikh-Tropsha criteria.

Cross-validation

>>> from qsarkit.models import QSARRegressor
>>> from qsarkit.validation import CrossValidator
>>> X, y = demo_fingerprints(512), DEMO_Y
>>> model = QSARRegressor("rf", random_state=0)
>>> report = CrossValidator(n_splits=5, random_state=0).evaluate(model, X, y)
>>> round(report["q2"], 1), round(report["rmse_cv"], 1)
(0.7, 0.6)

The report carries the out-of-fold predictions, so any further statistic or plot needs no refitting:

>>> report["y_pred_cv"].shape
(24,)
>>> report["method"], report["n_splits"]
('kfold', 5)

n_splits comes from the splitter, not the constructor argument — leave-one-out on 24 compounds is 24 splits, and a report claiming 5 would be wrong:

>>> CrossValidator(method="loo").evaluate(QSARRegressor("ridge"), X, y)["n_splits"]
24

"repeated_kfold" averages over several partitions, which matters on small datasets where one k-fold estimate is dominated by the luck of the split. "leave_group_out" is the one to use when compounds come in groups — a scaffold series, an assay batch, a source publication.

Robustness: y-scrambling

Refit the model on randomly permuted activities. If the scrambled models score anywhere near the real one, the apparent performance came from the model’s flexibility relative to the dataset size, not from a structure-activity relationship.

>>> from qsarkit.validation import YScrambling
>>> scramble = YScrambling(n_iterations=50, random_state=0).run(model, X, y)
>>> round(scramble["real_score"], 3)
0.953
>>> round(scramble["mean_scrambled_score"], 3)
0.839
>>> scramble["p_value"] < 0.05
True

Danger

Read those numbers again. The model fits randomly permuted activities to \(R^2 = 0.84\) on average, against 0.95 on the real ones. The p-value clears 0.05, but the honest reading is that most of this model’s apparent fit is capacity: 24 compounds described by 512 features will fit almost anything.

This is exactly the failure y-scrambling exists to expose, and it is invisible in the training \(R^2\) that would otherwise be reported.

The p-value can never be exactly zero — a permutation test cannot distinguish “very unlikely” from “impossible”, so reporting 0 would claim more than was measured:

>>> small = YScrambling(n_iterations=20, random_state=0).run(model, X, y)
>>> round(small["p_value"], 4) == round(1 / 21, 4)
True

YScrambling.plot() shows the scrambled distribution with the real score marked, which is more informative than either number alone:

>>> scrambler = YScrambling(n_iterations=20, random_state=0)
>>> _ = scrambler.run(model, X, y)
>>> type(scrambler.plot()).__name__
'Figure'

Precision: the bootstrap

A single cross-validated \(Q^2\) is one number with no error bar. Resampling the training set and scoring out-of-bag gives the spread:

>>> from qsarkit.validation import BootstrapValidator
>>> boot = BootstrapValidator(n_iterations=30, random_state=0).run(model, X, y)
>>> round(boot["mean_score"], 2)
0.34
>>> round(boot["ci_upper"] - boot["ci_lower"], 1)
1.9

An interval nearly two \(R^2\) units wide. Any comparison between two models on this dataset that turns on less than that is noise — and the interval is the only thing that says so.

Scoring is out-of-bag, not in-bag: about 36.8% of the data is left out of each resample, and scoring there rather than on the fitted rows is what makes this an estimate of generalization instead of of fit.

Predictivity: external validation

>>> from qsarkit.model_selection import RandomSplitter
>>> from qsarkit.validation import ExternalValidator
>>> train, test = next(RandomSplitter(test_size=0.25, random_state=0).split(X, y))
>>> cv = CrossValidator(n_splits=5, random_state=0).evaluate(model, X[train], y[train])
>>> fitted = QSARRegressor("rf", random_state=0).fit(X[train], y[train])
>>> result = ExternalValidator(q2=cv["q2"]).validate(
...     fitted, X[test], y[test], y[train])
>>> round(result["r2"], 2), round(result["q2_f1"], 2)
(0.82, 0.82)

Supply y_train so Q²F1 is scaled by the training set variance, which is what makes it comparable across differently-centred test sets. Supply q2 so Golbraikh-Tropsha criterion 1 can be evaluated rather than reporting None:

>>> result["golbraikh_tropsha"]["passed"]
False
>>> gt = result["golbraikh_tropsha"]
>>> [k for k in sorted(gt) if k.startswith("criterion") and gt[k] is False]
['criterion_1_q2', 'criterion_3_r0']

An \(R^2\) of 0.82 on the test set looks respectable and would have been reported as a success. Criterion 1 fails because the cross-validated \(Q^2\) of 0.24 is below the 0.5 threshold; criterion 3 concerns regression through the origin — the predictions correlate with the truth but are systematically offset. Running the full check is what turns a respectable-looking number into an accurate picture.

Note

None of these validators mutates the estimator you hand them: each clones it before fitting, so the same configured model can be passed to all four.

>>> template = QSARRegressor("ridge")
>>> _ = YScrambling(n_iterations=5, random_state=0).run(template, X, y)
>>> hasattr(template, "estimator_")
False

Choosing the metric

All four validators take scoring. It defaults to \(R^2\), which is right for a regression QSAR and wrong for everything else: a toxicity classifier has to be argued in ROC-AUC or average precision, and a regulator asking for RMSE is not asking for \(R^2\) reported next to it.

>>> from qsarkit.validation import available_metrics
>>> len(available_metrics())
18
>>> [m for m in available_metrics() if "auc" in m]
['pr_auc', 'roc_auc']

Pass several and every score becomes an array in the order given, so one pass reports them all:

>>> from qsarkit.models import QSARRegressor
>>> from qsarkit.validation import CrossValidator
>>> cv = CrossValidator(n_splits=5, random_state=0, scoring=["r2", "rmse", "mae"])
>>> result = cv.evaluate(QSARRegressor("rf", random_state=0), X, y)
>>> result["metric"]
('r2', 'rmse', 'mae')
>>> result["score"].round(2)
array([0.69, 0.6 , 0.44])

One metric returns a float; an iterable returns an array even when it holds a single entry, so adding a second metric never changes the shape of your code.

Classification metrics see probabilities

roc_auc, pr_auc and brier rank or calibrate, so they are given predict_proba’s positive-class column, never a thresholded label. Thresholding first throws away the ranking that ROC-AUC exists to measure.

>>> import numpy as np
>>> from qsarkit.models import QSARClassifier
>>> labels = (DEMO_Y > np.median(DEMO_Y)).astype(int)
>>> model = QSARClassifier("rf", random_state=0)
>>> cv = CrossValidator(n_splits=5, random_state=0,
...                     scoring=["roc_auc", "pr_auc", "mcc"])
>>> report = cv.evaluate(model, X, labels)
>>> report["score"].round(2)
array([0.93, 0.94, 0.75])
>>> report["y_score_cv"].shape          # out-of-fold probabilities
(24,)

For a metric that is not on the list, or one that needs probabilities:

>>> from sklearn.metrics import average_precision_score
>>> from qsarkit.validation import make_scorer
>>> scorer = make_scorer(average_precision_score, needs_proba=True, name="ap")
>>> CrossValidator(n_splits=5, random_state=0, scoring=scorer).evaluate(
...     model, X, labels)["metric"]
'ap'

A bare callable is assumed to take (y_true, y_pred) and to improve as it grows; make_scorer() is how the other cases are declared. Declaring a loss matters more than it looks: greater_is_better=False is what keeps “did the scrambled model do at least as well” comparing in the right direction, so a y-randomization p-value computed on RMSE is not reported backwards.

>>> from qsarkit.validation import YScrambling
>>> scramble = YScrambling(n_iterations=20, random_state=0,
...                        scoring=["r2", "rmse"]).run(
...     QSARRegressor("rf", random_state=0), X, y)
>>> scramble["p_value"].round(3)        # same verdict from a gain and a loss
array([0.048, 0.048])

Warning

Score out of fold before reading anything into a ranking metric. YScrambling scores the apparent, in-sample fit by default, which is what earlier releases did and what makes the classic over-fitting demonstration work for \(R^2\). It cannot work for ROC-AUC: a random forest separates permuted labels in-sample as perfectly as real ones, so both sides read near 1.0 and the test reports nothing.

>>> in_sample = YScrambling(n_iterations=20, random_state=0,
...                         scoring="roc_auc").run(model, X, labels)
>>> round(in_sample["real_score"], 2), round(in_sample["mean_scrambled_score"], 2)
(1.0, 1.0)

Pass cv — and stratify=True on an imbalanced endpoint — and the same test becomes informative:

>>> honest = YScrambling(n_iterations=20, random_state=0, scoring="roc_auc",
...                      cv=5, stratify=True).run(model, X, labels)
>>> round(honest["real_score"], 2), round(honest["mean_scrambled_score"], 2)
(0.93, 0.43)
>>> honest["scored_out_of_fold"]
True

See OECD validation for how these fit together into a reportable validation, and Metrics for the statistics they compute.

API

Model validation against the OECD principles.

Every validator here takes a scoring argument. It defaults to \(R^2\), which is the right default for a regression QSAR and wrong for everything else – a toxicity classifier has to be argued in ROC-AUC or average precision. Name a metric from available_metrics(), pass a (y_true, y_pred) callable, or wrap one with make_scorer() when it needs probabilities or is a loss. Pass several metrics and every score in the result becomes an array in the order given, so one pass reports them all:

CrossValidator(scoring=["roc_auc", "pr_auc", "mcc"]).evaluate(model, X, y)

References

class qsarkit.validation.CrossValidator(method='kfold', n_splits=5, n_repeats=10, random_state=None, scoring=None)[source]

Bases: object

Cross-validated Q^2, RMSE_CV and MAE_CV for a scikit-learn estimator.

Supports the four fold-splitting schemes routinely used in QSAR validation studies: k-fold, leave-one-out, repeated k-fold (to average away the fold-assignment randomness of a single k-fold run) and leave-one-group-out (for scaffold- or assay-based splits passed via groups). Every fold clones the estimator with sklearn.base.clone before fitting, so the caller’s own estimator object is never mutated.

Parameters:

  • method (Literal['kfold', 'loo', 'repeated_kfold', 'leave_group_out']) – Cross-validation scheme.

  • n_splits (int) – Number of folds for "kfold" and "repeated_kfold". Ignored by "loo" and "leave_group_out" (both use one fold per sample / per group).

  • n_repeats (int) – Number of repeats for "repeated_kfold". Ignored otherwise.

  • random_state (Optional[int]) – Seed controlling the fold shuffling of "kfold" and "repeated_kfold".

Examples

>>> import numpy as np
>>> from sklearn.datasets import make_regression
>>> from sklearn.linear_model import Ridge
>>> from qsarkit.validation import CrossValidator
>>> X, y = make_regression(n_samples=60, n_features=5, noise=1.0, random_state=0)
>>> result = CrossValidator(method="kfold", n_splits=5, random_state=0).evaluate(
...     Ridge(), X, y
... )
>>> result["q2"] > 0.5
True

References

method: Literal['kfold', 'loo', 'repeated_kfold', 'leave_group_out']
n_splits: int
n_repeats: int
random_state: int | None
evaluate(estimator, X, y, groups=None)[source]

Cross-validate estimator and compute Q^2, RMSE_CV and MAE_CV.

Parameters:
Returns:

"q2" (float), "rmse_cv" (float), "mae_cv" (float), "y_pred_cv" (ndarray of out-of-fold predictions aligned to the original row order – for "repeated_kfold" this is the average out-of-fold prediction over all repeats), "method" (str) and "n_splits" (int – the number of folds actually run, which for "loo" is the sample count and for "leave_group_out" the number of groups).

Return type:

Dict[str, Any]

Raises:

ValueError – If method="leave_group_out" and groups is None, or if method is not one of the four supported schemes.

Examples

>>> import numpy as np
>>> from sklearn.datasets import make_regression
>>> from sklearn.linear_model import Ridge
>>> from qsarkit.validation import CrossValidator
>>> X, y = make_regression(n_samples=50, n_features=4, noise=1.0, random_state=0)
>>> result = CrossValidator(method="loo").evaluate(Ridge(), X, y)
>>> sorted(result)
['mae_cv', 'method', 'n_splits', 'q2', 'rmse_cv', 'y_pred_cv']
class qsarkit.validation.YScrambling(n_iterations=100, random_state=None, scoring=None, cv=None, stratify=False)[source]

Bases: object

Test whether a model can fit randomly permuted labels as well as real ones.

Also called y-randomization. Refit the model many times on shuffled activities: if the scrambled models score anywhere near the real one, the apparent performance came from the model’s flexibility relative to the dataset size, not from a structure-activity relationship.

This is the check that catches the classic QSAR failure – a few dozen compounds described by thousands of descriptors, where something will always correlate. It is required evidence under OECD principle 4, and it is cheap, so there is no excuse for omitting it.

Parameters:

  • n_iterations (int) – Number of permutations. The smallest p-value obtainable is 1 / (n_iterations + 1), so 100 iterations cannot report anything below 0.0099.

  • random_state (Optional[int]) – Seed for the permutations.

  • scoring (Union[None, str, Scorer, Callable[..., float], Iterable[Union[str, Scorer, Callable[..., float]]]]) – The metric to argue in. Defaults to \(R^2\). Name one of available_metrics(), pass a (y_true, y_pred) callable, or use make_scorer() for a metric that needs probabilities or is a loss. Pass several and every score in the result becomes an array in the order given.

  • cv (Optional[int]) – Score out of fold over this many folds instead of on the training data. Strongly recommended for any flexible model, and required for a ranking metric to mean anything: a random forest reaches an in-sample ROC-AUC near 1.0 on permuted labels just as it does on real ones, so the in-sample comparison shows no gap and the test reports nothing. The default is None – the apparent, in-sample fit – because that is what earlier releases computed.

  • stratify (bool) – Use stratified folds when cv is set. Needed on an imbalanced classification endpoint, where an unstratified fold can contain no positives at all.

Variables:

  • real_score (float or ndarray) – The model’s score on the true labels: a float for one metric, an array in the given order for several.

  • scrambled_scores (ndarray) – Shape (n_iterations,) for one metric, (n_iterations, n_metrics) for several.

Examples

A real relationship survives the test:

>>> import numpy as np
>>> from sklearn.linear_model import Ridge
>>> from qsarkit.validation import YScrambling
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(60, 4))
>>> y = X[:, 0] * 3.0 + rng.normal(scale=0.2, size=60)
>>> result = YScrambling(n_iterations=50, random_state=0).run(Ridge(), X, y)
>>> result["p_value"] < 0.05
True

Pure noise does not:

>>> y_noise = rng.normal(size=60)
>>> noise = YScrambling(n_iterations=50, random_state=0).run(Ridge(), X, y_noise)
>>> noise["p_value"] > 0.05
True

The diagnostic value is in the gap between the two scores:

>>> round(result["real_score"] - result["mean_scrambled_score"], 2) > 0.5
True

References

real_score_: Any
scrambled_scores_: ndarray[tuple[Any, ...], dtype[float64]]
run(estimator, X, y)[source]

Fit on the real labels and on n_iterations permutations of them.

Parameters:
Returns:

real_score, mean_scrambled_score, std_scrambled_score, max_scrambled_score, best_scrambled_score (the largest for a metric where more is better, the smallest for a loss), p_value (the fraction of permutations scoring at least as well as the real fit, with the conventional +1 correction), n_iterations, metric (the name, or the tuple of names) and scored_out_of_fold.

Every score is a float when one metric was requested and an ndarray in the requested order when several were.

Return type:

Dict[str, Any]

Raises:

ValueError – If n_iterations is not positive.

plot(title='y-scrambling')[source]

Histogram of scrambled scores with the real score marked.

Parameters:

title (str)

Return type:

Any

Raises:

ModelNotFittedError – If run() has not been called.

class qsarkit.validation.ExternalValidator(q2=None, scoring=None)[source]

Bases: object

Score a fitted model on a held-out set with QSAR-appropriate metrics.

OECD principle 4’s predictivity requirement. Wraps qsar_regression_report() and the Golbraikh-Tropsha criteria so an external evaluation reports the same statistics every time, rather than whichever ones happened to look best.

Parameters:

  • q2 (Optional[float]) – A cross-validated Q² from the training set. Supply it so that Golbraikh-Tropsha criterion 1 can be evaluated; without it that criterion reports None rather than silently passing.

  • scoring (Union[None, str, Scorer, Callable[..., float], Iterable[Union[str, Scorer, Callable[..., float]]]]) – Report these metrics instead of the regression report. The default (None) keeps the QSAR regression report and the Golbraikh-Tropsha criteria, which is what a regression submission needs; naming metrics is how a classification endpoint is validated, and then score and metric replace the report.

Examples

>>> import numpy as np
>>> from sklearn.linear_model import Ridge
>>> from qsarkit.validation import ExternalValidator
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(60, 4))
>>> y = X[:, 0] * 3.0 + rng.normal(scale=0.2, size=60)
>>> model = Ridge().fit(X[:45], y[:45])
>>> report = ExternalValidator(q2=0.9).validate(model, X[45:], y[45:], y[:45])
>>> report["r2"] > 0.9
True
>>> report["golbraikh_tropsha"]["passed"]
True

References

validate(estimator, X_test, y_test, y_train=None)[source]

Evaluate a fitted model on the test set.

Parameters:
Returns:

The regression report, plus q2_f1 when y_train is given and golbraikh_tropsha.

Return type:

Dict[str, Any]

class qsarkit.validation.BootstrapValidator(n_iterations=100, random_state=None, scoring=None)[source]

Bases: object

Bootstrap the training set to estimate how stable a score is.

A single cross-validated Q² is one number with no error bar. Resampling the training set with replacement and refitting gives the spread, which is what tells you whether a 0.02 difference between two models means anything on this much data – usually it does not.

Parameters:

  • n_iterations (int) – Number of bootstrap resamples.

  • random_state (Optional[int]) – Seed.

Variables:

scores (ndarray of shape (n_iterations,)) – Out-of-bag score from each resample.

Examples

>>> import numpy as np
>>> from sklearn.linear_model import Ridge
>>> from qsarkit.validation import BootstrapValidator
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(60, 4))
>>> y = X[:, 0] * 3.0 + rng.normal(scale=0.2, size=60)
>>> result = BootstrapValidator(n_iterations=25, random_state=0).run(Ridge(), X, y)
>>> result["mean_score"] > 0.9
True
>>> result["ci_lower"] <= result["mean_score"] <= result["ci_upper"]
True

References

scores_: ndarray[tuple[Any, ...], dtype[float64]]
run(estimator, X, y, confidence=0.95)[source]

Refit on bootstrap resamples and score on the out-of-bag remainder.

Parameters:
Returns:

mean_score, std_score, ci_lower, ci_upper, confidence, n_iterations, n_effective (resamples that produced a usable out-of-bag set) and metric.

Every score is a float when one metric was requested and an ndarray in the requested order when several were.

Return type:

Dict[str, Any]

Raises:

ValueError – If n_iterations is not positive, confidence is not in (0, 1), or no resample left any out-of-bag samples.

class qsarkit.validation.Scorer(name, func, needs_proba=False, greater_is_better=True)[source]

Bases: object

One metric, with what a validation method needs to know about it.

Variables:

  • name (str) – How the metric is reported.

  • func (callable) – func(y_true, y_pred_or_score) -> float.

  • needs_proba (bool) – Whether func wants predict_proba(X)[:, 1] rather than predict(X).

  • greater_is_better (bool) – Whether a larger value is a better model. False for losses such as RMSE, which is what keeps a permutation p-value the right way round.

name: str
func: Callable[[...], float]
needs_proba: bool
greater_is_better: bool
is_at_least_as_good_as(candidate, reference)[source]

Whether candidate matches or beats reference for this metric.

Return type:

bool

qsarkit.validation.make_scorer(func, *, needs_proba=False, greater_is_better=True, name=None)[source]

Wrap a metric for use as scoring=.

Parameters:

  • func (Callable[..., float]) – func(y_true, y_pred) -> float, or func(y_true, y_score) when needs_proba is set.

  • needs_proba (bool) – Pass the positive-class probability instead of the predicted label.

  • greater_is_better (bool) – Set False for a loss, so that “at least as good” compares the right way round.

  • name (Optional[str]) – Defaults to the function’s __name__.

Return type:

Scorer

Examples

A metric a validation method can use directly:

>>> from sklearn.metrics import average_precision_score
>>> from qsarkit.validation import make_scorer
>>> scorer = make_scorer(average_precision_score, needs_proba=True)
>>> scorer.name
'average_precision_score'

A loss, declared as one:

>>> from qsarkit.metrics import rmse
>>> make_scorer(rmse, greater_is_better=False).greater_is_better
False
qsarkit.validation.available_metrics()[source]

The metric names scoring= accepts.

Return type:

List[str]

Examples

>>> from qsarkit.validation import available_metrics
>>> names = available_metrics()
>>> "r2" in names, "roc_auc" in names, "rmse" in names
(True, True, True)

References

  • OECD (2007). Guidance Document on the Validation of (Quantitative) Structure-Activity Relationship [(Q)SAR] Models, ENV/JM/MONO(2007)2. doi:10.1787/9789264085442-en

  • Golbraikh, A. & Tropsha, A. (2002). “Beware of q2!” J. Mol. Graph. Model., 20(4), 269-276. doi:10.1016/S1093-3263(01)00123-1

  • Rücker, C., Rücker, G. & Meringer, M. (2007). “y-Randomization and Its Variants in QSPR/QSAR.” J. Chem. Inf. Model., 47(6), 2345-2357. doi:10.1021/ci700157b

  • Tropsha, A., Gramatica, P. & Gombar, V. K. (2003). “The Importance of Being Earnest.” QSAR Comb. Sci., 22(1), 69-77. doi:10.1002/qsar.200390007

  • Efron, B. & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall. doi:10.1201/9780429246593

  • Consonni, V., Ballabio, D. & Todeschini, R. (2009). “Comments on the Definition of the Q2 Parameter for QSAR Validation.” J. Chem. Inf. Model., 49(7), 1669-1678. doi:10.1021/ci900115y

  • Varma, S. & Simon, R. (2006). “Bias in Error Estimation When Using Cross-Validation for Model Selection.” BMC Bioinformatics, 7, 91. doi:10.1186/1471-2105-7-91