Metrics

QSAR-specific performance measures. R² alone does not establish that a model predicts — these are the statistics the QSAR literature and the OECD guidance actually ask for.

Regression

>>> from qsarkit.metrics import ccc, q2_f1, qsar_regression_report, rmse
>>> from qsarkit.models import QSARRegressor
>>> X, y = demo_fingerprints(256), DEMO_Y
>>> y_pred = QSARRegressor("rf", random_state=0).fit(X, y).predict(X)
>>> round(rmse(y, y_pred), 2)
0.24
>>> round(ccc(y, y_pred), 3)
0.972

Q²F1–F3 differ in what they compare the model against. The choice matters: an external Q² computed against the training mean flatters a model whose test set happens to be centred differently.

>>> round(q2_f1(y, y_pred, y), 3)
0.952

The Golbraikh–Tropsha criteria

Five conditions a predictive QSAR model must satisfy simultaneously. They exist because a high R² is achievable by a model with a systematic slope error, which the regression-through-origin terms catch.

>>> from qsarkit.metrics import golbraikh_tropsha_criteria
>>> result = golbraikh_tropsha_criteria(y, y_pred)
>>> result["passed"]
True
>>> [k for k in sorted(result) if k.startswith("criterion")]
['criterion_1_q2', 'criterion_2_r2', 'criterion_3_r0', 'criterion_4_slope', 'criterion_5_delta_r0']

Criterion 1 needs a cross-validated Q², which this call was not given, so it reports None rather than quietly passing:

>>> result["criterion_1_q2"] is None, result["q2_available"]
(True, False)

The full report

>>> report = qsar_regression_report(y, y_pred)
>>> sorted(report)
['average_r2m', 'ccc', 'delta_r2m', 'golbraikh_tropsha', 'mae', 'q2_f2', 'r2', 'rmse']

Classification and virtual screening

>>> import numpy as np
>>> from qsarkit.metrics import bedroc, enrichment_factor, roc_auc
>>> scores = np.linspace(1.0, 0.0, 100)
>>> labels = np.zeros(100); labels[:10] = 1        # actives ranked first
>>> round(roc_auc(labels, scores), 3)
1.0
>>> round(enrichment_factor(labels, scores, fraction=0.1), 2)
10.0

ROC AUC weights every rank equally, which is the wrong emphasis for screening: only the top of the list will ever be tested. BEDROC applies an exponential weight so early recognition dominates.

>>> round(bedroc(labels, scores, alpha=20.0), 3)
1.0

API

QSAR-specific regression and classification metrics.

Covers the classical statistics (RMSE, MAE, R^2, MCC, …) plus the QSAR-specific external-predictivity coefficients (Q^2_F1/F2/F3, CCC, r_m^2, the Golbraikh-Tropsha criteria) and the early-recognition metrics used in virtual screening (enrichment factor, RIE, BEDROC).

qsarkit.metrics.calibration_curve(y_true, y_prob, n_bins=10, strategy='uniform')[source]

Observed frequency against predicted probability, per bin.

The data behind a reliability diagram. A perfectly calibrated model lies on the diagonal: among the compounds it scored 0.7, 70% are active.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Binary labels, 0 or 1.

  • y_prob (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Predicted probability of the positive class.

  • n_bins (int) – Number of bins.

  • strategy (Literal['uniform', 'quantile']) – "uniform" splits [0, 1] into equal-width bins, which shows where on the probability scale the model is wrong. "quantile" puts an equal number of samples in each bin, which gives every point a comparable error bar – the better choice when predictions cluster near 0, as they do in virtual screening.

Returns:

mean_predicted and observed_frequency (one entry per non-empty bin), counts, and bin_edges.

Return type:

Dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If the inputs are not a matching pair of binary labels and probabilities, or n_bins is below 2.

Examples

A perfectly calibrated set of predictions lies on the diagonal:

>>> import numpy as np
>>> from qsarkit.metrics import calibration_curve
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(size=4000)
>>> y = (rng.uniform(size=4000) < p).astype(int)
>>> curve = calibration_curve(y, p, n_bins=5)
>>> bool(np.allclose(curve["mean_predicted"], curve["observed_frequency"],
...                  atol=0.05))
True

An over-confident model bends away from it:

>>> squashed = np.clip(p * 1.6 - 0.3, 0, 1)
>>> curve = calibration_curve(y, squashed, n_bins=5)
>>> bool((curve["observed_frequency"][0] > curve["mean_predicted"][0]))
True

References

qsarkit.metrics.expected_calibration_error(y_true, y_prob, n_bins=10, strategy='uniform')[source]

Sample-weighted mean gap between confidence and accuracy.

\[\mathrm{ECE} = \sum_{b=1}^{B} \frac{n_b}{N} \bigl| \bar{p}_b - \bar{y}_b \bigr|\]

where \(\bar{p}_b\) is the mean predicted probability in bin \(b\), \(\bar{y}_b\) the observed frequency, and \(n_b\) the bin’s size. 0 is perfect.

Parameters:
Returns:

In [0, 1].

Return type:

float

Notes

ECE depends on the binning, and a model can lower it by concentrating its predictions rather than by improving. Read it beside the curve from calibration_curve(), not on its own.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import expected_calibration_error
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(size=4000)
>>> y = (rng.uniform(size=4000) < p).astype(int)
>>> round(expected_calibration_error(y, p, n_bins=10), 2) < 0.05
True

A model whose probabilities are all shifted upward scores worse:

>>> shifted = np.clip(p + 0.25, 0, 1)
>>> expected_calibration_error(y, shifted) > expected_calibration_error(y, p)
True

References

qsarkit.metrics.maximum_calibration_error(y_true, y_prob, n_bins=10, strategy='uniform', min_count=1)[source]

Largest single-bin gap between confidence and accuracy.

The worst case rather than the average. Useful when a decision will be made at one particular probability: an ECE of 0.02 is no comfort if the bin you actually threshold on is off by 0.3.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Binary labels.

  • y_prob (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Predicted probabilities.

  • n_bins (int)

  • strategy (Literal['uniform', 'quantile'])

  • min_count (int) – Ignore bins with fewer samples than this. A bin holding two compounds can only report frequencies of 0, 0.5 or 1, so it produces a large gap by arithmetic rather than by miscalibration; raising this suppresses that artefact.

Returns:

In [0, 1]. 0.0 when no bin meets min_count.

Return type:

float

Examples

>>> import numpy as np
>>> from qsarkit.metrics import (
...     expected_calibration_error, maximum_calibration_error)
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(size=4000)
>>> y = (rng.uniform(size=4000) < p).astype(int)
>>> mce = maximum_calibration_error(y, p, n_bins=10, min_count=20)
>>> ece = expected_calibration_error(y, p, n_bins=10)
>>> mce >= ece        # the worst bin is at least as bad as the average
True

References

qsarkit.metrics.calibration_report(y_true, y_prob, n_bins=10, strategy='uniform')[source]

Everything needed to judge whether probabilities can be believed.

Parameters:
Returns:

ece, mce, brier, brier_skill_score, mean_predicted, observed_frequency, base_rate, n_samples and n_bins_used.

Return type:

Dict[str, Any]

Notes

brier_skill_score compares the Brier score against always predicting the base rate: positive means the model beats that baseline, 0 or below means it does not. On an imbalanced screening set a raw Brier score near 0.05 looks excellent and is often worse than the constant prediction – the skill score is what exposes that.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import calibration_report
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(size=2000)
>>> y = (rng.uniform(size=2000) < p).astype(int)
>>> report = calibration_report(y, p, n_bins=10)
>>> report["ece"] < 0.05
True
>>> report["brier_skill_score"] > 0
True

A constant prediction at the base rate has no skill at all:

>>> flat = np.full(2000, y.mean())
>>> round(calibration_report(y, flat)["brier_skill_score"], 6)
0.0

References

qsarkit.metrics.qq_data(residuals, standardize=True)[source]

Theoretical against observed quantiles, for a normal Q-Q plot.

Every regression metric here – RMSE, \(R^2\), the Golbraikh-Tropsha criteria – assumes roughly normal, homoscedastic errors. A Q-Q plot is the quickest check. Points on the diagonal mean normal residuals; an S-shape means heavy tails; a curve at one end means skew, usually from a handful of badly mispredicted compounds that RMSE alone will not name.

Parameters:

  • residuals (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Observed minus predicted. Non-finite entries are dropped.

  • standardize (bool) – Divide by the standard deviation, so the reference line is \(y = x\) regardless of the residuals’ scale.

Returns:

theoretical_quantiles, sample_quantiles (both sorted ascending) and reference_line as (slope, intercept).

Return type:

Dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If fewer than three finite residuals remain.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import qq_data
>>> rng = np.random.default_rng(0)
>>> data = qq_data(rng.normal(size=500))
>>> corr = np.corrcoef(data["theoretical_quantiles"],
...                    data["sample_quantiles"])[0, 1]
>>> bool(corr > 0.99)            # normal residuals track the diagonal
True

Heavy-tailed residuals do not:

>>> heavy = qq_data(rng.standard_t(df=2, size=500))
>>> bool(np.corrcoef(heavy["theoretical_quantiles"],
...                  heavy["sample_quantiles"])[0, 1] < corr)
True

References

  • Wilk, M. B. & Gnanadesikan, R. (1968). “Probability Plotting Methods for the Analysis of Data.” Biometrika, 55(1), 1-17. https://doi.org/10.1093/biomet/55.1.1

  • Blom, G. (1958). “Statistical Estimates and Transformed Beta-Variables.” Wiley.

qsarkit.metrics.residual_normality(y_true, y_pred)[source]

Test whether regression residuals are normal and homoscedastic.

Parameters:
Returns:

shapiro_statistic and shapiro_p (normality; p below 0.05 argues against it), skew, excess_kurtosis, heteroscedasticity_r (Spearman correlation between the fitted value and the absolute residual – non-zero means error size depends on the prediction), and n.

Return type:

Dict[str, float]

Notes

Shapiro-Wilk cuts both ways at QSAR sample sizes, and neither direction should be read as a verdict:

  • On a few dozen compounds it has little power, so shapiro_p > 0.05 is weak evidence of normality rather than a clearance.

  • On a few hundred it starts rejecting samples that are normal, for departures far too small to affect an RMSE.

The skew, kurtosis and heteroscedasticity terms describe the shape rather than testing a hypothesis, which is more useful here – and the Q-Q plot from qq_data() more useful still, because it shows where the departure is.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import residual_normality
>>> rng = np.random.default_rng(0)
>>> truth = rng.normal(size=500)
>>> report = residual_normality(truth, truth + rng.normal(scale=0.1, size=500))
>>> report["shapiro_p"] > 0.05         # residuals are normal
True
>>> abs(report["skew"]) < 0.2
True

A badly skewed residual distribution is caught by every term at once:

>>> skewed = residual_normality(np.zeros(300), -rng.exponential(size=300))
>>> skewed["shapiro_p"] < 0.001, abs(skewed["skew"]) > 1
(True, True)

Errors that grow with the prediction show up as heteroscedasticity:

>>> fitted = np.linspace(1, 10, 300)
>>> noisy = fitted + rng.normal(scale=fitted * 0.3)
>>> abs(residual_normality(noisy, fitted)["heteroscedasticity_r"]) > 0.2
True

References

qsarkit.metrics.threshold_sweep(y_true, y_score, n_thresholds=None, pos_label=None)[source]

Every operating point of a binary classifier, as arrays.

Evaluates the confusion matrix at each candidate threshold, so a criterion can be maximized or a curve plotted without refitting.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Binary labels. Numeric or boolean labels need no further specification – the larger value is the positive class.

  • y_score (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Scores or probabilities. Only their order matters.

  • n_thresholds (Optional[int]) – Evaluate this many evenly-spaced quantiles of the score instead of every distinct value. Use it on large sets, where the default gives one column per unique score.

  • pos_label (Optional[Any]) – Which label is the positive class. Required for string labels.

Returns:

thresholds and, aligned with it, tp, fp, tn, fn, sensitivity (recall, TPR), specificity, precision, f1, mcc, balanced_accuracy, youden_j and accuracy.

Return type:

Dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If the inputs mismatch, or do not contain exactly two classes.

Notes

A prediction is positive when score >= threshold, so the sweep includes a threshold above every score (predict nothing) to make the degenerate end of the curve explicit rather than absent.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import threshold_sweep
>>> y = np.array([0, 0, 1, 1])
>>> scores = np.array([0.1, 0.4, 0.6, 0.9])
>>> sweep = threshold_sweep(y, scores)
>>> best = int(np.argmax(sweep["youden_j"]))
>>> float(sweep["thresholds"][best]), float(sweep["youden_j"][best])
(0.6, 1.0)

Sensitivity falls and specificity rises as the threshold increases:

>>> bool(np.all(np.diff(sweep["sensitivity"]) <= 0))
True
>>> bool(np.all(np.diff(sweep["specificity"]) >= 0))
True

References

qsarkit.metrics.optimal_threshold(y_true, y_score, criterion='youden', cost_fn=1.0, cost_fp=1.0, min_precision=None, min_recall=None, n_thresholds=None, pos_label=None)[source]

Choose the decision threshold that best serves a stated objective.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Binary labels.

  • y_score (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Scores or probabilities.

  • criterion (Literal['youden', 'f1', 'mcc', 'balanced_accuracy', 'cost', 'precision', 'recall']) –

    What “best” means:

    "youden"

    Maximize Youden’s \(J = \mathrm{sensitivity} + \mathrm{specificity} - 1\). Treats both error types as equally costly and both classes as equally important; the usual default when you have no cost information.

    "mcc"

    Maximize the Matthews correlation coefficient. The most informative single number on imbalanced data, because it uses all four cells of the confusion matrix.

    "f1"

    Maximize F1. Ignores true negatives, so it suits screening where the inactive majority is uninteresting.

    "balanced_accuracy"

    Maximize the mean of sensitivity and specificity.

    "cost"

    Minimize cost_fn * FN + cost_fp * FP. The honest choice when the two errors have different consequences.

    "precision" / "recall"

    Maximize the other of the pair subject to min_precision or min_recall. For “find me 200 compounds to test, as pure as possible” and its mirror image.

  • cost_fn (float) – Relative cost of a false negative and a false positive. Used by criterion="cost".

  • cost_fp (float) – Relative cost of a false negative and a false positive. Used by criterion="cost".

  • min_precision (Optional[float]) – Required with criterion="precision": the floor precision must clear, among which recall is maximized.

  • min_recall (Optional[float]) – Required with criterion="recall".

  • n_thresholds (Optional[int]) – Passed to threshold_sweep().

  • pos_label (Optional[Any]) – Which label is the positive class. Required for string labels.

Returns:

threshold, the criterion used, its score, and the confusion matrix and derived rates at that threshold.

Return type:

Dict[str, Any]

Raises:

ValueError – If criterion is unknown, a required constraint is missing, or no threshold satisfies the constraint.

Notes

Select the threshold on validation data, never on the test set. A threshold tuned on the same data it is scored on is a fitted parameter, and the resulting performance is optimistic in exactly the way an untuned 0.5 cut is not.

Examples

>>> import numpy as np
>>> from qsarkit.metrics import optimal_threshold
>>> y = np.array([0, 0, 0, 0, 1, 1])
>>> scores = np.array([0.1, 0.2, 0.3, 0.55, 0.6, 0.8])
>>> best = optimal_threshold(y, scores, criterion="youden")
>>> float(best["threshold"]), round(best["score"], 3)
(0.6, 1.0)

On an imbalanced set the best threshold is nowhere near 0.5:

>>> rng = np.random.default_rng(0)
>>> y = np.zeros(1000, dtype=int); y[:30] = 1
>>> scores = rng.beta(2, 8, size=1000) + y * 0.25
>>> chosen = optimal_threshold(y, scores, criterion="mcc")
>>> bool(chosen["threshold"] < 0.5)
True

Asymmetric costs move it. Making a missed active ten times as expensive as a false positive lowers the bar:

>>> cheap = optimal_threshold(y, scores, criterion="cost",
...                           cost_fn=1.0, cost_fp=1.0)
>>> costly = optimal_threshold(y, scores, criterion="cost",
...                            cost_fn=10.0, cost_fp=1.0)
>>> bool(costly["threshold"] <= cheap["threshold"])
True

And a hard constraint is respected rather than traded away:

>>> pure = optimal_threshold(y, scores, criterion="precision",
...                          min_precision=0.3)
>>> bool(pure["precision"] >= 0.3)
True

An unreachable constraint is reported, with the best actually available, rather than quietly relaxed:

>>> optimal_threshold(y, scores, criterion="precision", min_precision=0.9)
Traceback (most recent call last):
    ...
ValueError: No threshold reaches precision 0.9. The best available is ...

References

qsarkit.metrics.threshold_report(y_true, y_score, n_thresholds=None, pos_label=None)[source]

Compare what each criterion would choose, against the 0.5 default.

The most useful output when you do not yet know which criterion you want: it shows how much the choice actually matters on your data, and what the untuned 0.5 cut is costing.

Parameters:
Returns:

One entry per criterion (youden, mcc, f1, balanced_accuracy), plus default_0.5 evaluated at the conventional cut, base_rate, and roc_auc/pr_auc for the threshold-free picture.

Return type:

Dict[str, Any]

Examples

>>> import numpy as np
>>> from qsarkit.metrics import threshold_report
>>> rng = np.random.default_rng(0)
>>> y = np.zeros(600, dtype=int); y[:30] = 1
>>> scores = rng.beta(2, 8, size=600) + y * 0.3
>>> report = threshold_report(y, scores)
>>> sorted(k for k in report if isinstance(report[k], dict))
['balanced_accuracy', 'default_0.5', 'f1', 'mcc', 'youden']
>>> report["youden"]["mcc"] > report["default_0.5"]["mcc"]
True

The default cut can miss nearly every active on imbalanced data:

>>> report["default_0.5"]["sensitivity"] < report["youden"]["sensitivity"]
True

References

qsarkit.metrics.mse(y_true, y_pred)[source]

Mean squared error.

Parameters:
Returns:

mean((y_true - y_pred) ** 2).

Return type:

float

Examples

>>> from qsarkit.metrics import mse
>>> round(mse([1.0, 2.0, 3.0], [1.0, 2.0, 5.0]), 4)
1.3333

References

qsarkit.metrics.rmse(y_true, y_pred)[source]

Root mean squared error.

Parameters:
Returns:

sqrt(mean((y_true - y_pred) ** 2)), in the units of y.

Return type:

float

Examples

>>> from qsarkit.metrics import rmse
>>> rmse([1.0, 2.0, 3.0], [2.0, 3.0, 4.0])
1.0

References

qsarkit.metrics.rmsep(y_true, y_pred)[source]

Root mean squared error of prediction (RMSE on an external set).

Numerically identical to rmse(); the separate name is retained because the QSAR literature consistently distinguishes RMSEC (calibration/training), RMSECV (cross-validation) and RMSEP (external prediction), and reporting code reads better when the intent is explicit.

Parameters:
Returns:

The root mean squared error of prediction.

Return type:

float

Examples

>>> from qsarkit.metrics import rmsep
>>> rmsep([1.0, 2.0], [1.5, 2.5])
0.5

References

qsarkit.metrics.mae(y_true, y_pred)[source]

Mean absolute error.

Parameters:
Returns:

mean(|y_true - y_pred|).

Return type:

float

Examples

>>> from qsarkit.metrics import mae
>>> mae([1.0, 2.0, 3.0], [1.0, 4.0, 3.0])
0.6666666666666666

References

qsarkit.metrics.median_ae(y_true, y_pred)[source]

Median absolute error (outlier-robust location of the error).

Parameters:
Returns:

median(|y_true - y_pred|).

Return type:

float

Examples

>>> from qsarkit.metrics import median_ae
>>> median_ae([1.0, 2.0, 3.0], [1.0, 2.0, 30.0])
0.0

References

qsarkit.metrics.bias(y_true, y_pred)[source]

Systematic error (mean signed residual mean(y_pred - y_true)).

A non-zero bias indicates the model systematically over- (positive) or under-predicts (negative) the endpoint, which RMSE alone hides.

Parameters:
Returns:

The mean signed residual.

Return type:

float

Examples

>>> from qsarkit.metrics import bias
>>> bias([1.0, 2.0, 3.0], [2.0, 3.0, 4.0])
1.0

References

qsarkit.metrics.press(y_true, y_pred)[source]

Predictive residual sum of squares, sum((y_true - y_pred) ** 2).

Parameters:
Returns:

The PRESS statistic.

Return type:

float

Examples

>>> from qsarkit.metrics import press
>>> press([1.0, 2.0, 3.0], [1.0, 2.0, 5.0])
4.0

References

qsarkit.metrics.see(y_true, y_pred, n_parameters=0)[source]

Standard error of estimate, sqrt(RSS / (n - p - 1)).

Parameters:
Returns:

The standard error of estimate.

Return type:

float

Raises:

ValueError – If n_samples - n_parameters - 1 <= 0.

Examples

>>> from qsarkit.metrics import see
>>> round(see([1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 5.0], n_parameters=1), 4)
0.7071

References

qsarkit.metrics.r2_score(y_true, y_pred)[source]

Coefficient of determination R^2 (fraction of variance explained).

R^2 = 1 - RSS / TSS where TSS uses the mean of y_true. Note this is not the squared Pearson correlation unless the predictions are unbiased and unit-slope; QSAR papers that report “R^2” for a test set usually mean this quantity (equivalently q2_f2()).

Parameters:
Returns:

The coefficient of determination. Can be negative.

Return type:

float

Raises:

ValueError – If y_true has zero variance.

Examples

>>> from qsarkit.metrics import r2_score
>>> r2_score([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.adjusted_r2_score(y_true, y_pred, n_features)[source]

R^2 penalized for the number of descriptors in the model.

R^2_adj = 1 - (1 - R^2) * (n - 1) / (n - p - 1).

Parameters:
Returns:

The adjusted coefficient of determination.

Return type:

float

Raises:

ValueError – If n_samples - n_features - 1 <= 0.

Examples

>>> from qsarkit.metrics import adjusted_r2_score
>>> round(adjusted_r2_score([1.0, 2.0, 3.0, 4.5], [1.1, 2.0, 2.9, 4.4], 1), 3)
0.993

References

qsarkit.metrics.ccc(y_true, y_pred)[source]

Lin’s concordance correlation coefficient.

CCC = 2 * cov(y, yhat) / (var(y) + var(yhat) + (mean(y) - mean(yhat))^2)

CCC simultaneously penalizes loss of precision (correlation) and loss of accuracy (deviation from the 45-degree line), which makes it stricter than the Pearson correlation and a recommended external-validation statistic for QSAR. Population (biased, ddof=0) moments are used, as in the original paper.

Parameters:
Returns:

The concordance correlation coefficient, in [-1, 1].

Return type:

float

Examples

>>> from qsarkit.metrics import ccc
>>> ccc([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

  • Lin, L. I.-K. (1989). “A Concordance Correlation Coefficient to Evaluate Reproducibility.” Biometrics, 45(1), 255-268. https://doi.org/10.2307/2532051

qsarkit.metrics.q2_f1(y_true, y_pred, y_train)[source]

External predictivity Q^2_F1 (Shi/Schuurmann formulation).

Q^2_F1 = 1 - sum((y_ext - yhat_ext)^2) / sum((y_ext - mean(y_train))^2)

The residual sum of squares of the external set is normalized by its variation around the training set mean. This is the original “R^2_pred” of Shi et al. and is the most commonly reported form.

Parameters:
Returns:

The Q^2_F1 statistic.

Return type:

float

Raises:

ValueError – If the normalizing sum of squares is zero.

Examples

>>> from qsarkit.metrics import q2_f1
>>> q2_f1([2.0, 4.0], [2.0, 3.0], y_train=[0.0, 2.0, 4.0])
0.75

References

qsarkit.metrics.q2_f2(y_true, y_pred)[source]

External predictivity Q^2_F2 (Schuurmann formulation).

Q^2_F2 = 1 - sum((y_ext - yhat_ext)^2) / sum((y_ext - mean(y_ext))^2)

Uses the external set mean, making it identical to the ordinary coefficient of determination computed on the test set. It is the most conservative of the three and is invariant to the training/test split ratio.

Parameters:
Returns:

The Q^2_F2 statistic.

Return type:

float

Examples

>>> from qsarkit.metrics import q2_f2
>>> q2_f2([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.q2_f3(y_true, y_pred, y_train)[source]

External predictivity Q^2_F3 (Consonni/Todeschini formulation).

Q^2_F3 = 1 - [sum((y_ext - yhat_ext)^2) / n_ext]
           / [sum((y_train - mean(y_train))^2) / n_train]

Both sums of squares are divided by their own sample size, which makes the statistic independent of the size of the external set and (per Consonni et al.) monotonically related to RMSEP for a fixed training set.

Parameters:
Returns:

The Q^2_F3 statistic.

Return type:

float

Raises:

ValueError – If y_train has zero variance.

Examples

>>> from qsarkit.metrics import q2_f3
>>> q2_f3([2.0, 4.0], [2.0, 3.0], y_train=[0.0, 2.0, 4.0])
0.8125

References

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

qsarkit.metrics.k_slope(y_true, y_pred)[source]

Slope k of the regression of predicted on observed through the origin.

k = sum(y * yhat) / sum(y^2)

Parameters:
Returns:

The through-origin slope k. Golbraikh & Tropsha require 0.85 <= k <= 1.15.

Return type:

float

Raises:

ValueError – If sum(y_true ** 2) is zero.

Examples

>>> from qsarkit.metrics import k_slope
>>> k_slope([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.k_prime_slope(y_true, y_pred)[source]

Slope k' of the regression of observed on predicted through the origin.

k' = sum(y * yhat) / sum(yhat^2)

Parameters:
Returns:

The through-origin slope k'.

Return type:

float

Raises:

ValueError – If sum(y_pred ** 2) is zero.

Examples

>>> from qsarkit.metrics import k_prime_slope
>>> k_prime_slope([1.0, 2.0, 3.0], [2.0, 4.0, 6.0])
0.5

References

qsarkit.metrics.r0_squared(y_true, y_pred)[source]

Determination coefficient of the through-origin fit yhat = k * y.

R0^2 = 1 - sum((yhat - k*y)^2) / sum((yhat - mean(yhat))^2)

Parameters:
Returns:

The through-origin R0^2.

Return type:

float

Examples

>>> from qsarkit.metrics import r0_squared
>>> r0_squared([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.r0_prime_squared(y_true, y_pred)[source]

Determination coefficient of the through-origin fit y = k' * yhat.

R0'^2 = 1 - sum((y - k'*yhat)^2) / sum((y - mean(y))^2)

This is the axis-swapped counterpart of r0_squared().

Parameters:
Returns:

The through-origin R0'^2.

Return type:

float

Examples

>>> from qsarkit.metrics import r0_prime_squared
>>> r0_prime_squared([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.r2m(y_true, y_pred)[source]

Roy’s r_m^2 metric, r^2 * (1 - sqrt(|r^2 - r0^2|)).

r^2 is the squared Pearson correlation and r0^2 the through-origin determination coefficient (r0_squared()). A model is considered acceptable when the average of r_m^2 and r_m'^2 exceeds 0.5.

Parameters:
Returns:

The r_m^2 value.

Return type:

float

Examples

>>> from qsarkit.metrics import r2m
>>> r2m([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.r2m_prime(y_true, y_pred)[source]

Roy’s axis-swapped r_m'^2, r^2 * (1 - sqrt(|r^2 - r0'^2|)).

Parameters:
Returns:

The r_m'^2 value.

Return type:

float

Examples

>>> from qsarkit.metrics import r2m_prime
>>> r2m_prime([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.average_r2m(y_true, y_pred)[source]

Mean of r_m^2 and r_m'^2; should exceed 0.5.

Parameters:
Returns:

(r_m^2 + r_m'^2) / 2.

Return type:

float

Examples

>>> from qsarkit.metrics import average_r2m
>>> average_r2m([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
1.0

References

qsarkit.metrics.delta_r2m(y_true, y_pred)[source]

Absolute difference |r_m^2 - r_m'^2|; should stay below 0.2.

A large delta signals that the observed-vs-predicted relationship is strongly asymmetric, i.e. the model is systematically compressing or expanding the response range.

Parameters:
Returns:

|r_m^2 - r_m'^2|.

Return type:

float

Examples

>>> from qsarkit.metrics import delta_r2m
>>> delta_r2m([1.0, 2.0, 3.0], [1.0, 2.0, 3.0])
0.0

References

qsarkit.metrics.golbraikh_tropsha_criteria(y_true, y_pred, q2=None)[source]

Evaluate the five Golbraikh-Tropsha acceptability criteria.

A QSAR model is deemed externally predictive when all of

  1. q^2 > 0.5 (leave-one-out / cross-validated Q^2 of the training set);

  2. r^2 > 0.6 (squared correlation between observed and predicted on the external set);

  3. (r^2 - R0^2) / r^2 < 0.1 or (r^2 - R0'^2) / r^2 < 0.1;

  4. 0.85 <= k <= 1.15 or 0.85 <= k' <= 1.15;

  5. |R0^2 - R0'^2| < 0.3

hold simultaneously.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Observed and predicted responses of the external validation set.

  • y_pred (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Observed and predicted responses of the external validation set.

  • q2 (Optional[float]) – Cross-validated Q^2 of the training set. If None, criterion 1 is reported as None and excluded from passed (the returned dict then also carries "q2_available": False).

Returns:

Keys r2, r0_squared, r0_prime_squared, k, k_prime, q2, the booleans criterion_1_q2, criterion_2_r2, criterion_3_r0, criterion_4_slope, criterion_5_delta_r0, q2_available and the overall passed.

Return type:

Dict[str, Any]

Examples

>>> from qsarkit.metrics import golbraikh_tropsha_criteria
>>> res = golbraikh_tropsha_criteria([1.0, 2.0, 3.0, 4.0],
...                                  [1.1, 2.0, 2.9, 4.05], q2=0.9)
>>> res["passed"]
True

References

qsarkit.metrics.confusion_counts(y_true, y_pred)[source]

Return the binary confusion-matrix counts as a dictionary.

Parameters:
Returns:

Keys "tp", "tn", "fp", "fn".

Return type:

Dict[str, int]

Examples

>>> from qsarkit.metrics import confusion_counts
>>> confusion_counts([1, 1, 0, 0], [1, 0, 0, 0]) == {
...     "tp": 1, "tn": 2, "fp": 0, "fn": 1}
True

References

qsarkit.metrics.accuracy(y_true, y_pred)[source]

Fraction of correctly classified samples.

Parameters:
Returns:

(TP + TN) / n.

Return type:

float

Examples

>>> from qsarkit.metrics import accuracy
>>> accuracy([1, 1, 0, 0], [1, 0, 0, 0])
0.75

References

qsarkit.metrics.balanced_accuracy(y_true, y_pred)[source]

Mean of sensitivity and specificity.

Preferred over plain accuracy for the strongly imbalanced datasets that are the norm in QSAR classification (e.g. toxicity endpoints with 5% actives), where a trivial majority classifier already scores high accuracy but only 0.5 balanced accuracy.

Parameters:
Returns:

(sensitivity + specificity) / 2.

Return type:

float

Examples

>>> from qsarkit.metrics import balanced_accuracy
>>> balanced_accuracy([1, 1, 0, 0], [1, 0, 0, 0])
0.75

References

  • Brodersen, K. H., Ong, C. S., Stephan, K. E. & Buhmann, J. M. (2010). “The Balanced Accuracy and Its Posterior Distribution.” ICPR 2010, 3121-3124. https://doi.org/10.1109/ICPR.2010.764

qsarkit.metrics.sensitivity(y_true, y_pred)[source]

True-positive rate TP / (TP + FN) (recall of the active class).

Parameters:
Returns:

The sensitivity; 0.0 when there are no positives.

Return type:

float

Examples

>>> from qsarkit.metrics import sensitivity
>>> sensitivity([1, 1, 0, 0], [1, 0, 0, 0])
0.5

References

qsarkit.metrics.specificity(y_true, y_pred)[source]

True-negative rate TN / (TN + FP).

Parameters:
Returns:

The specificity; 0.0 when there are no negatives.

Return type:

float

Examples

>>> from qsarkit.metrics import specificity
>>> specificity([1, 1, 0, 0], [1, 0, 0, 0])
1.0

References

qsarkit.metrics.precision(y_true, y_pred)[source]

Positive predictive value TP / (TP + FP).

Parameters:
Returns:

The precision; 0.0 when nothing is predicted positive.

Return type:

float

Examples

>>> from qsarkit.metrics import precision
>>> precision([1, 1, 0, 0], [1, 0, 0, 0])
1.0

References

qsarkit.metrics.recall(y_true, y_pred)[source]

Alias of sensitivity().

Parameters:
Returns:

TP / (TP + FN).

Return type:

float

Examples

>>> from qsarkit.metrics import recall
>>> recall([1, 1, 0, 0], [1, 0, 0, 0])
0.5

References

qsarkit.metrics.f1_score(y_true, y_pred)[source]

Harmonic mean of precision and recall.

Parameters:
Returns:

2 * P * R / (P + R); 0.0 when both are zero.

Return type:

float

Examples

>>> from qsarkit.metrics import f1_score
>>> round(f1_score([1, 1, 0, 0], [1, 0, 0, 0]), 4)
0.6667

References

qsarkit.metrics.matthews_corrcoef(y_true, y_pred)[source]

Matthews correlation coefficient (MCC).

MCC = (TP*TN - FP*FN) / sqrt((TP+FP)(TP+FN)(TN+FP)(TN+FN))

MCC is a correlation coefficient between observed and predicted binary classifications; it is high only when all four confusion-matrix quadrants are good, which is why it is the recommended single-number summary for imbalanced QSAR classification.

Parameters:
Returns:

The MCC in [-1, 1]; 0.0 if any marginal is degenerate.

Return type:

float

Examples

>>> from qsarkit.metrics import matthews_corrcoef
>>> matthews_corrcoef([1, 1, 0, 0], [1, 1, 0, 0])
1.0

References

qsarkit.metrics.cohen_kappa(y_true, y_pred)[source]

Cohen’s kappa: agreement corrected for chance.

kappa = (p_o - p_e) / (1 - p_e) where p_o is the observed agreement and p_e the agreement expected from the marginal label frequencies.

Parameters:
Returns:

Cohen’s kappa; 1.0 when observed agreement is perfect and chance agreement is degenerate (p_e == 1).

Return type:

float

Examples

>>> from qsarkit.metrics import cohen_kappa
>>> cohen_kappa([1, 1, 0, 0], [1, 1, 0, 0])
1.0

References

qsarkit.metrics.roc_auc(y_true, y_score)[source]

Area under the receiver-operating-characteristic curve.

Parameters:
Returns:

The ROC-AUC.

Return type:

float

Raises:

ValueError – If only one class is present in y_true.

Examples

>>> from qsarkit.metrics import roc_auc
>>> roc_auc([0, 0, 1, 1], [0.1, 0.2, 0.8, 0.9])
1.0

References

qsarkit.metrics.pr_auc(y_true, y_score)[source]

Area under the precision-recall curve (average precision).

More informative than ROC-AUC when actives are rare, because the precision axis is sensitive to the large number of true negatives that ROC-AUC dilutes away.

Parameters:
Returns:

The average precision.

Return type:

float

Raises:

ValueError – If only one class is present in y_true.

Examples

>>> from qsarkit.metrics import pr_auc
>>> pr_auc([0, 0, 1, 1], [0.1, 0.2, 0.8, 0.9])
1.0

References

qsarkit.metrics.brier_score(y_true, y_prob)[source]

Brier score: mean squared error of the predicted probabilities.

Measures calibration as well as discrimination; lower is better.

Parameters:
Returns:

mean((y_prob - y_true) ** 2).

Return type:

float

Examples

>>> from qsarkit.metrics import brier_score
>>> brier_score([0, 1], [0.0, 1.0])
0.0

References

qsarkit.metrics.enrichment_factor(y_true, y_score, fraction=0.01)[source]

Enrichment factor at a given fraction of the ranked list.

EF(chi) = (n_actives_in_top / n_top) / (n_actives_total / N)

An EF of 10 at 1% means the screen finds ten times as many actives in the top 1% as random selection would.

Parameters:

  • y_true (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Binary ground-truth labels in {0, 1}.

  • y_score (Union[Buffer, _SupportsArray[dtype[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], complex, bytes, str, _NestedSequence[complex | bytes | str]]) – Continuous scores; the list is ranked in decreasing score order.

  • fraction (float) – Fraction chi of the ranked list to inspect, in (0, 1]. The top-k cut-off is max(1, round(fraction * N)).

Returns:

The enrichment factor. Its maximum attainable value is min(1 / fraction, N / n_actives).

Return type:

float

Raises:

ValueError – If fraction is outside (0, 1] or there are no actives.

Examples

>>> from qsarkit.metrics import enrichment_factor
>>> enrichment_factor([1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
...                   [0.9, 0.8, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1],
...                   fraction=0.2)
5.0

References

qsarkit.metrics.robust_initial_enhancement(y_true, y_score, alpha=20.0)[source]

Robust initial enhancement (RIE) of Sheridan et al.

RIE = sum_i exp(-alpha * r_i / N)
      / [ (n/N) * (1 - exp(-alpha)) / (exp(alpha/N) - 1) ]

where r_i are the 1-based ranks of the n actives among N compounds. The exponential weight makes RIE a continuous, threshold-free generalization of the enrichment factor: alpha sets how sharply early ranks are rewarded (the top 1/alpha of the list carries most of the weight). RIE is 1 for a random ranking.

Parameters:
Returns:

The RIE value.

Return type:

float

Raises:

ValueError – If alpha <= 0 or there are no actives.

Examples

>>> from qsarkit.metrics import robust_initial_enhancement
>>> scores = [1.0, 0.9, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
>>> labels = [1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
>>> robust_initial_enhancement(labels, scores, alpha=20.0) > 1.0
True

References

  • Sheridan, R. P., Singh, S. B., Fluder, E. M. & Kearsley, S. K. (2001). “Protocols for Bridging the Peptide to Nonpeptide Gap in Topological Similarity Searches.” J. Chem. Inf. Comput. Sci., 41(5), 1395-1406. https://doi.org/10.1021/ci0100144

  • Truchon, J.-F. & Bayly, C. I. (2007). J. Chem. Inf. Model., 47(2), 488-508. https://doi.org/10.1021/ci600426e

qsarkit.metrics.bedroc(y_true, y_score, alpha=20.0)[source]

Boltzmann-enhanced discrimination of ROC (BEDROC).

BEDROC rescales robust_initial_enhancement() onto [0, 1], removing RIE’s dependence on the fraction of actives:

BEDROC = RIE * Ra * sinh(alpha/2)
         / (cosh(alpha/2) - cosh(alpha/2 - alpha*Ra))
         + 1 / (1 - exp(alpha * (1 - Ra)))

with Ra = n_actives / N. A perfect early-recognition ranking gives ~1, a random ranking gives ~``Ra`` and the worst ranking ~0. The default alpha=20 concentrates 80% of the weight in the top 8% of the list, the usual choice in virtual-screening benchmarks.

Parameters:
Returns:

The BEDROC score.

Return type:

float

Raises:

ValueError – If alpha <= 0, there are no actives, or every compound is active (Ra == 1, for which the metric is undefined).

Examples

>>> from qsarkit.metrics import bedroc
>>> labels = [1] * 5 + [0] * 95
>>> scores = list(range(100, 0, -1))
>>> bedroc(labels, scores, alpha=20.0) > 0.99
True

References

  • Truchon, J.-F. & Bayly, C. I. (2007). “Evaluating Virtual Screening Methods: Good and Bad Metrics for the ‘Early Recognition’ Problem.” J. Chem. Inf. Model., 47(2), 488-508. https://doi.org/10.1021/ci600426e

qsarkit.metrics.qsar_regression_report(y_true, y_pred, y_train=None, q2=None)[source]

Compute the standard bundle of QSAR regression validation statistics.

Parameters:
Returns:

Keys "r2", "rmse", "mae", "ccc", "q2_f2", "average_r2m", "delta_r2m", "golbraikh_tropsha" and, when y_train is given, "q2_f1" and "q2_f3".

Return type:

Dict[str, Any]

Examples

>>> from qsarkit.metrics import qsar_regression_report
>>> report = qsar_regression_report(
...     [1.0, 2.0, 3.0, 4.0], [1.1, 2.0, 2.9, 4.05],
...     y_train=[0.0, 1.0, 2.0, 3.0, 4.0], q2=0.9,
... )
>>> report["golbraikh_tropsha"]["passed"]
True

References

qsarkit.metrics.qsar_classification_report(y_true, y_pred, y_score=None)[source]

Compute the standard bundle of QSAR classification validation statistics.

Parameters:
Returns:

Keys "confusion", "accuracy", "balanced_accuracy", "sensitivity", "specificity", "precision", "recall", "f1", "mcc", "cohen_kappa" and, when y_score is given, "roc_auc" and "pr_auc".

Return type:

Dict[str, Any]

Examples

>>> from qsarkit.metrics import qsar_classification_report
>>> report = qsar_classification_report(
...     [1, 1, 0, 0], [1, 0, 0, 0], y_score=[0.9, 0.4, 0.2, 0.1],
... )
>>> report["mcc"] > 0
True

References

References

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

  • 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

  • Roy, K. et al. (2012). “Comparative Studies on Some Metrics for External Validation of QSPR Models.” J. Chem. Inf. Model., 52(2), 396-408. doi:10.1021/ci200520g

  • Truchon, J.-F. & Bayly, C. I. (2007). “Evaluating Virtual Screening Methods: Good and Bad Metrics for the Early Recognition Problem.” J. Chem. Inf. Model., 47(2), 488-508. doi:10.1021/ci600426e

  • Lin, L. I. (1989). “A Concordance Correlation Coefficient to Evaluate Reproducibility.” Biometrics, 45(1), 255-268. doi:10.2307/2532051