Results comparison framework for benchmarking using Cinnabar#

Comparing two sets of free energy predictions by eye is deceptively difficult. Two methods can show different RMSE or MUE values on a plot yet be statistically indistinguishable given the size of the dataset. Cinnabar addresses this with a rigorous comparison framework built on a joint bootstrapping method to build a distribution of metric differences and identify statistically significant differences between methods.

This tutorial walks through the methodology and demonstrates how to use the compare_and_rank_results function with example RBFE data.

Note: The comparison methodology implemented here is inspired by the approach described in Valsson et al.

Understanding the comparison method#

The comparison is built around a joint bootstrapping procedure that generates a distribution of differences in the chosen evaluation metric. The full workflow can be broken down into the following steps:

  1. Define the evaluation metric (e.g. RMSE, MUE) and the sets of results to compare — either edgewise (ΔΔG) or nodewise (ΔG).

  2. Bootstrap jointly draw paired resamples of the data and recompute the evaluation metric for every source on the same resample. This preserves the correlation structure between methods and avoids artificially inflating differences.

  3. Compute two-sided p-values from the fraction of bootstrap differences that cross zero.

  4. Apply the Holm multiple testing correction when comparing more than two sets of results.

  5. Report confidence intervals for the metric differences directly from the bootstrap distribution.

  6. Rank and group results using the compact letter display (CLD) via the insert-absorb algorithm. Methods that share a letter are not significantly different; methods that share no letters are significantly different.

Limitations#

Important: RBFE calculations are inherently correlated because a weakly connected network is required to estimate ΔG for each ligand. Bootstrapping edges therefore tends to underestimate confidence intervals as the effective sample size is closer to N_ligands than to N_edges as assumed by the method.
Important: ΔG estimates back-calculated from ΔΔGs via the MLE estimator are not independent across ligands as the estimator uses information from all edges simultaneously. Nodewise comparisons are most appropriate for absolute binding free energy (ABFE) predictions, where each ligand's estimate is obtained independently.

Requirements#

  • Experimental measurements must be available for every ligand in the FEMap.

  • All sources must provide predictions for exactly the same set of ligands or ligand pairs.

Loading example RBFE data#

For this tutorial we use RBFE data generated with OpenFE which is included with cinnabar’s test suite. You can replace these files with your own outputs. See the cinnabar API tutorial for details on building an FEMap from different data formats.

The two data files have a simple tabular format:

! head ../cinnabar/data/experimental_data.csv

Ligand,expt_DG,expt_dDG
CAT-13a,-8.83,0.10
CAT-13b,-9.11,0.10
CAT-13c,-9.31,0.10
CAT-13d,-10.46,0.10
CAT-13e,-9.95,0.10
CAT-13f,-9.08,0.10
CAT-13g,-9.08,0.10
CAT-13h,-9.62,0.10
CAT-13i,-9.26,0.10

! head ../cinnabar/data/computational_data.csv

Ligand1,Ligand2,calc_DDG,calc_dDDG(MBAR),calc_dDDG(additional)
CAT-13b,CAT-17g,0.36,0.11,0.0
CAT-13a,CAT-17g,-0.02,0.1,0.0
CAT-13e,CAT-17g,1.5,0.11,0.0
CAT-4m,CAT-4c,0.78,0.1,0.0
CAT-13k,CAT-4d,-0.59,0.11,0.0
CAT-24,CAT-17e,1.98,0.08,0.0
CAT-13g,CAT-17g,0.86,0.15,0.0
CAT-13d,CAT-13h,1.46,0.1,0.0
CAT-13a,CAT-17i,-0.76,0.11,0.0

# fmt: off
from cinnabar.compare import compare_and_rank_results
from cinnabar.femap import FEMap

%matplotlib inline
import numpy as np
import pandas as pd
from openff.units import unit

femap = FEMap()

# load the computational results (source 1)
rbfe_results = pd.read_csv("../cinnabar/data/computational_data.csv")

for _, result in rbfe_results.iterrows():
    femap.add_relative_calculation(
        labelA=result["Ligand1"],
        labelB=result["Ligand2"],
        value=result["calc_DDG"] * unit.kilocalorie_per_mole,
        uncertainty=result["calc_dDDG(MBAR)"] * unit.kilocalorie_per_mole,
        source="OpenFE",
    )

# load the experimental values
experimental_results = pd.read_csv("../cinnabar/data/experimental_data.csv")

for _, exp_row in experimental_results.iterrows():
    femap.add_experimental_measurement(
        label=exp_row["Ligand"],
        value=exp_row["expt_DG"] * unit.kilocalorie_per_mole,
        uncertainty=exp_row["expt_dDG"] * unit.kilocalorie_per_mole,
        source="Experimental",
    )
# fmt: on

Creating additional result sets for comparison#

In a real benchmarking study you would load a second (and third) set of calculated values here from a different force field, a different simulation protocol, or a different software package. You just need to add them to the same FEMap with a distinct source string.

For this tutorial we generate two synthetic alternatives by perturbing the original values with random noise, giving us full control over how similar or different the methods appear.

np.random.seed(42)  # for reproducibility

# Source 2: slight perturbation (similar accuracy to source 1)
for _, result in rbfe_results.iterrows():
    femap.add_relative_calculation(
        labelA=result["Ligand1"],
        labelB=result["Ligand2"],
        value=np.random.normal(loc=result["calc_DDG"], scale=result["calc_dDDG(MBAR)"]) * unit.kilocalorie_per_mole,
        uncertainty=result["calc_dDDG(MBAR)"] * unit.kilocalorie_per_mole,
        source="OpenFE_perturbed",
    )

# Source 3: large perturbation (noticeably worse accuracy)
for _, result in rbfe_results.iterrows():
    femap.add_relative_calculation(
        labelA=result["Ligand1"],
        labelB=result["Ligand2"],
        value=np.random.normal(loc=result["calc_DDG"], scale=12.0 * result["calc_dDDG(MBAR)"])
        * unit.kilocalorie_per_mole,
        uncertainty=result["calc_dDDG(MBAR)"] * unit.kilocalorie_per_mole,
        source="OpenFE_noisy",
    )

Running the comparison#

We now call compare_and_rank_results with the FEMap containing all three sources. Key parameters:

Parameter

Description

prediction_type

"edgewise" (ΔΔG) or "nodewise" (ΔG)

rank_metric

The metric used to rank and compare methods.

metrics_to_compute

Metrics reported in the summary table. Defaults to ["MUE", "RMSE"] for edgewise.

num_bootstraps

Number of joint bootstrap resamples (1 000 is a sensible default).

confidence_level

Width of the reported confidence intervals, default 0.95 (95 %).

alpha

The significance level for determining statistical significance, default 0.05.

 
summary_df, comparison_df = compare_and_rank_results(
    femap=femap,
    prediction_type="edgewise",
    rank_metric="MUE",  # the metric used for ranking the results
    metrics_to_compute=["MUE", "RMSE"],  # metrics to report in the summary table
    num_bootstraps=1_000,
    confidence_level=0.95,  # report 95 % confidence intervals
    alpha=0.05,  # for statistical significance testing (e.g. p-value threshold)
)

The summary table#

The summary table has one row per source. For each requested metric you get:

  • the sample value (metric computed on all data points),

  • *_CI_Lower / *_CI_Upper: the bootstrap confidence interval bounds, and

  • CLD: the compact letter display assignment.

Sources that share a CLD letter are not significantly different from each other for the chosen rank_metric (after multiple testing correction). Sources with no letters in common are significantly different. For example, three methods labelled "a", "b", and "b" tell you: the first and second/third are significantly different from each other, but the second is not significantly different from the third method.

summary_df
Model MUE MUE_CI_Lower MUE_CI_Upper RMSE RMSE_CI_Lower RMSE_CI_Upper CLD
0 OpenFE 0.867586 0.710513 1.021414 1.053002 0.869701 1.221763 a
1 OpenFE_noisy 1.166930 0.971467 1.361760 1.386575 1.175243 1.577053 b
2 OpenFE_perturbed 0.859956 0.698599 1.023484 1.056510 0.867684 1.228928 a

The comparison table#

The comparison table has one row per unique source pair and records:

  • Diff in <rank_metric>: observed difference in the ranking metric (source 1 minus source 2, on all data),

  • CI Lower / CI Upper: bootstrap confidence interval around that difference,

  • p-value: two-sided bootstrap p-value,

  • p-value corrected: Holm-corrected p-value (appears automatically when >2 sources are present),

  • significant: True if the (corrected) p-value is below alpha.

comparison_df
Model 1 Model 2 Diff in MUE CI Lower CI Upper p-value significant p-value corrected
0 OpenFE OpenFE_noisy -0.299344 -0.511742 -0.088689 0.008 True 0.024
1 OpenFE OpenFE_perturbed 0.007630 -0.015944 0.030025 0.562 False 0.562
2 OpenFE_noisy OpenFE_perturbed 0.306974 0.096280 0.518630 0.008 True 0.024

Because we have three sources here the p-value corrected column is present along with the p-value from the bootstrap test. OpenFE_noisy — with 12× the per-edge noise is in a distinct CLD group from the other two, while OpenFE and OpenFE_perturbed share a letter as their difference is not large enough to be resolved at this sample size.

Tip: When comparing only two sources the p-value corrected column is absent (no correction needed) and the raw p-value drives the significant flag directly.
Important: In cases where one method clearly outperforms another the bootstrap p-value may be reported as 0.0 because none of the bootstrap differences cross zero. This is a limitation of the bootstrap testing method. In such cases the confidence interval around the difference metric should be used to interpret the significance of the difference instead and a p-value of 0.0 should not be reported as the level of significance.

Additional options#

Choosing the right metrics#

The available metrics depend on prediction_type:

Metric

Edgewise

Nodewise

Notes

MUE

✓

✓

Mean unsigned error

RMSE

✓

✓

Root mean squared error

RAE

✓

✓

Relative absolute error vs. a naïve mean predictor

R2

—

✓

R² (Pearson r²); not meaningful for relative ΔΔG data

rho

—

✓

Pearson r

KTAU

—

✓

Kendall’s τ

PI

—

✓

Predictive index (Pearlman et al.)

Correlation metrics (R², ρ, KTAU, PI) are only meaningful for nodewise ΔG comparisons where the data have an absolute scale. For edgewise comparisons the sign of a ΔΔG is arbitrary, so error metrics (MUE, RMSE) are the default.

If metrics_to_compute=None (the default), cinnabar automatically selects ["MUE", "RMSE"] for edgewise and ["MUE", "RMSE", "RAE", "R2", "rho", "KTAU", "PI"] for nodewise comparisons.

Adjusting the confidence level#

The confidence_level parameter (default 0.95) controls the width of the reported CIs. Pass a different value to tighten or loosen the intervals:

summary_df, comparison_df = compare_and_rank_results(
    femap=femap,
    rank_metric="MUE",
    metrics_to_compute=["MUE", "RMSE"],
    confidence_level=0.90,  # 90 % CIs — narrower, easier to detect differences
)

Recap#

  • compare_and_rank_results performs joint bootstrapping across all sources, ensuring paired, fair comparisons on the same resampled data.

  • Two tables are returned: a summary table (per-source metric values, CIs, and CLD letters) and a comparison table (pairwise differences, CIs, and p-values).

  • The CLD encodes statistical groupings: sources sharing a letter are not significantly different; sources with no shared letter are significantly different.

  • With more than two sources, Holm-corrected p-values appear automatically in p-value corrected and drive the significant flags.

  • Use prediction_type="nodewise" only for independent per-ligand ΔG estimates (e.g. ABFE), not for MLE-derived values from a relative network.