{ "cells": [ { "cell_type": "markdown", "id": "f07f14d1b152062e", "metadata": {}, "source": [ "# Results comparison framework for benchmarking using ``Cinnabar``\n", "\n", "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.\n", "\n", "This tutorial walks through the methodology and demonstrates how to use the ``compare_and_rank_results`` function with example RBFE data.\n", "\n", "
\n", "Note: The comparison methodology implemented here is inspired by the approach described in Valsson et al.\n", "
" ] }, { "cell_type": "markdown", "id": "f9e3ff514773b116", "metadata": {}, "source": [ "## Understanding the comparison method\n", "\n", "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:\n", "\n", "1. **Define** the evaluation metric (e.g. RMSE, MUE) and the sets of results to compare — either ``edgewise`` (ΔΔG) or ``nodewise`` (ΔG).\n", "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.\n", "3. **Compute two-sided p-values** from the fraction of bootstrap differences that cross zero.\n", "4. **Apply the Holm multiple testing correction** when comparing more than two sets of results.\n", "5. **Report confidence intervals** for the metric differences directly from the bootstrap distribution.\n", "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." ] }, { "cell_type": "markdown", "id": "88669d7f9af4abf8", "metadata": {}, "source": [ "## Limitations\n", "\n", "
\n", "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.\n", "
\n", "\n", "
\n", "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.\n", "
\n", "\n", "### Requirements\n", "\n", "- Experimental measurements must be available for every ligand in the ``FEMap``.\n", "- All sources must provide predictions for exactly the **same** set of ligands or ligand pairs." ] }, { "cell_type": "markdown", "id": "35186bc4db8fa6ca", "metadata": {}, "source": [ "## Loading example RBFE data\n", "\n", "For this tutorial we use RBFE data generated with [OpenFE](https://docs.openfree.energy/en/latest/) which is included with ``cinnabar``'s test suite. You can replace these files with your own outputs. See the {doc}`cinnabar API tutorial ` for details on building an ``FEMap`` from different data formats.\n", "\n", "The two data files have a simple tabular format:" ] }, { "cell_type": "code", "execution_count": 1, "id": "1d4764ff8e34299c", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:46.574932Z", "start_time": "2026-05-19T09:21:46.423124Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ligand,expt_DG,expt_dDG\n", "CAT-13a,-8.83,0.10\n", "CAT-13b,-9.11,0.10\n", "CAT-13c,-9.31,0.10\n", "CAT-13d,-10.46,0.10\n", "CAT-13e,-9.95,0.10\n", "CAT-13f,-9.08,0.10\n", "CAT-13g,-9.08,0.10\n", "CAT-13h,-9.62,0.10\n", "CAT-13i,-9.26,0.10\n" ] } ], "source": [ "! head ../cinnabar/data/experimental_data.csv" ] }, { "cell_type": "code", "execution_count": 2, "id": "68317fab7ba5d0c8", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:46.713070Z", "start_time": "2026-05-19T09:21:46.575967Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ligand1,Ligand2,calc_DDG,calc_dDDG(MBAR),calc_dDDG(additional)\n", "CAT-13b,CAT-17g,0.36,0.11,0.0\n", "CAT-13a,CAT-17g,-0.02,0.1,0.0\n", "CAT-13e,CAT-17g,1.5,0.11,0.0\n", "CAT-4m,CAT-4c,0.78,0.1,0.0\n", "CAT-13k,CAT-4d,-0.59,0.11,0.0\n", "CAT-24,CAT-17e,1.98,0.08,0.0\n", "CAT-13g,CAT-17g,0.86,0.15,0.0\n", "CAT-13d,CAT-13h,1.46,0.1,0.0\n", "CAT-13a,CAT-17i,-0.76,0.11,0.0\n" ] } ], "source": [ "! head ../cinnabar/data/computational_data.csv" ] }, { "cell_type": "code", "execution_count": 3, "id": "cd9c5ddc07ceeb9", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:47.429987Z", "start_time": "2026-05-19T09:21:46.714383Z" } }, "outputs": [], "source": [ "# fmt: off\n", "from cinnabar.compare import compare_and_rank_results\n", "from cinnabar.femap import FEMap\n", "\n", "%matplotlib inline\n", "import numpy as np\n", "import pandas as pd\n", "from openff.units import unit\n", "\n", "femap = FEMap()\n", "\n", "# load the computational results (source 1)\n", "rbfe_results = pd.read_csv(\"../cinnabar/data/computational_data.csv\")\n", "\n", "for _, result in rbfe_results.iterrows():\n", " femap.add_relative_calculation(\n", " labelA=result[\"Ligand1\"],\n", " labelB=result[\"Ligand2\"],\n", " value=result[\"calc_DDG\"] * unit.kilocalorie_per_mole,\n", " uncertainty=result[\"calc_dDDG(MBAR)\"] * unit.kilocalorie_per_mole,\n", " source=\"OpenFE\",\n", " )\n", "\n", "# load the experimental values\n", "experimental_results = pd.read_csv(\"../cinnabar/data/experimental_data.csv\")\n", "\n", "for _, exp_row in experimental_results.iterrows():\n", " femap.add_experimental_measurement(\n", " label=exp_row[\"Ligand\"],\n", " value=exp_row[\"expt_DG\"] * unit.kilocalorie_per_mole,\n", " uncertainty=exp_row[\"expt_dDG\"] * unit.kilocalorie_per_mole,\n", " source=\"Experimental\",\n", " )\n", "# fmt: on" ] }, { "cell_type": "markdown", "id": "7dc78019e981e2f9", "metadata": {}, "source": [ "## Creating additional result sets for comparison\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 4, "id": "71a8a6b537f85455", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:47.502421Z", "start_time": "2026-05-19T09:21:47.468442Z" } }, "outputs": [], "source": [ "np.random.seed(42) # for reproducibility\n", "\n", "# Source 2: slight perturbation (similar accuracy to source 1)\n", "for _, result in rbfe_results.iterrows():\n", " femap.add_relative_calculation(\n", " labelA=result[\"Ligand1\"],\n", " labelB=result[\"Ligand2\"],\n", " value=np.random.normal(loc=result[\"calc_DDG\"], scale=result[\"calc_dDDG(MBAR)\"]) * unit.kilocalorie_per_mole,\n", " uncertainty=result[\"calc_dDDG(MBAR)\"] * unit.kilocalorie_per_mole,\n", " source=\"OpenFE_perturbed\",\n", " )\n", "\n", "# Source 3: large perturbation (noticeably worse accuracy)\n", "for _, result in rbfe_results.iterrows():\n", " femap.add_relative_calculation(\n", " labelA=result[\"Ligand1\"],\n", " labelB=result[\"Ligand2\"],\n", " value=np.random.normal(loc=result[\"calc_DDG\"], scale=12.0 * result[\"calc_dDDG(MBAR)\"])\n", " * unit.kilocalorie_per_mole,\n", " uncertainty=result[\"calc_dDDG(MBAR)\"] * unit.kilocalorie_per_mole,\n", " source=\"OpenFE_noisy\",\n", " )" ] }, { "cell_type": "markdown", "id": "4df80ede3e08ddf0", "metadata": {}, "source": [ "## Running the comparison\n", "\n", "We now call ``compare_and_rank_results`` with the ``FEMap`` containing all three sources. Key parameters:\n", "\n", "| Parameter | Description |\n", "| --- |--------------------------------------------------------------------------------------|\n", "| ``prediction_type`` | ``\"edgewise\"`` (ΔΔG) or ``\"nodewise\"`` (ΔG) |\n", "| ``rank_metric`` | The metric used to rank and compare methods. |\n", "| ``metrics_to_compute`` | Metrics reported in the summary table. Defaults to ``[\"MUE\", \"RMSE\"]`` for edgewise. |\n", "| ``num_bootstraps`` | Number of joint bootstrap resamples (1 000 is a sensible default). |\n", "| ``confidence_level`` | Width of the reported confidence intervals, default ``0.95`` (95 %). |\n", "| ``alpha`` | The significance level for determining statistical significance, default ``0.05``. |\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "3d987cc6b1ddf37e", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:49.156781Z", "start_time": "2026-05-19T09:21:47.504921Z" } }, "outputs": [], "source": [ "summary_df, comparison_df = compare_and_rank_results(\n", " femap=femap,\n", " prediction_type=\"edgewise\",\n", " rank_metric=\"MUE\", # the metric used for ranking the results\n", " metrics_to_compute=[\"MUE\", \"RMSE\"], # metrics to report in the summary table\n", " num_bootstraps=1_000,\n", " confidence_level=0.95, # report 95 % confidence intervals\n", " alpha=0.05, # for statistical significance testing (e.g. p-value threshold)\n", ")" ] }, { "cell_type": "markdown", "id": "48fb1d4e6b173827", "metadata": {}, "source": [ "### The summary table\n", "\n", "The summary table has one row per source. For each requested metric you get:\n", "\n", "- the **sample value** (metric computed on all data points),\n", "- ``*_CI_Lower`` / ``*_CI_Upper``: the **bootstrap confidence interval** bounds, and\n", "- ``CLD``: the compact letter display assignment.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "id": "cb590e7befe1d980", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:49.183272Z", "start_time": "2026-05-19T09:21:49.157518Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ModelMUEMUE_CI_LowerMUE_CI_UpperRMSERMSE_CI_LowerRMSE_CI_UpperCLD
0OpenFE0.8675860.7105131.0214141.0530020.8697011.221763a
1OpenFE_noisy1.1669300.9714671.3617601.3865751.1752431.577053b
2OpenFE_perturbed0.8599560.6985991.0234841.0565100.8676841.228928a
\n", "
" ], "text/plain": [ " Model MUE MUE_CI_Lower MUE_CI_Upper RMSE \\\n", "0 OpenFE 0.867586 0.710513 1.021414 1.053002 \n", "1 OpenFE_noisy 1.166930 0.971467 1.361760 1.386575 \n", "2 OpenFE_perturbed 0.859956 0.698599 1.023484 1.056510 \n", "\n", " RMSE_CI_Lower RMSE_CI_Upper CLD \n", "0 0.869701 1.221763 a \n", "1 1.175243 1.577053 b \n", "2 0.867684 1.228928 a " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "summary_df" ] }, { "cell_type": "markdown", "id": "c9ec8f1dd33a19bb", "metadata": {}, "source": [ "### The comparison table\n", "\n", "The comparison table has one row per unique source pair and records:\n", "\n", "- ``Diff in ``: observed difference in the ranking metric (source 1 minus source 2, on all data),\n", "- ``CI Lower`` / ``CI Upper``: bootstrap confidence interval around that difference,\n", "- ``p-value``: two-sided bootstrap p-value,\n", "- ``p-value corrected``: Holm-corrected p-value *(appears automatically when >2 sources are present)*,\n", "- ``significant``: ``True`` if the (corrected) p-value is below ``alpha``." ] }, { "cell_type": "code", "execution_count": 7, "id": "330b38a6f59bf249", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:49.211626Z", "start_time": "2026-05-19T09:21:49.190516Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Model 1Model 2Diff in MUECI LowerCI Upperp-valuesignificantp-value corrected
0OpenFEOpenFE_noisy-0.299344-0.511742-0.0886890.008True0.024
1OpenFEOpenFE_perturbed0.007630-0.0159440.0300250.562False0.562
2OpenFE_noisyOpenFE_perturbed0.3069740.0962800.5186300.008True0.024
\n", "
" ], "text/plain": [ " Model 1 Model 2 Diff in MUE CI Lower CI Upper p-value \\\n", "0 OpenFE OpenFE_noisy -0.299344 -0.511742 -0.088689 0.008 \n", "1 OpenFE OpenFE_perturbed 0.007630 -0.015944 0.030025 0.562 \n", "2 OpenFE_noisy OpenFE_perturbed 0.306974 0.096280 0.518630 0.008 \n", "\n", " significant p-value corrected \n", "0 True 0.024 \n", "1 False 0.562 \n", "2 True 0.024 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "comparison_df" ] }, { "cell_type": "markdown", "id": "46a1d694c36bee59", "metadata": {}, "source": [ "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.\n", "\n", "
\n", "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.\n", "
\n", "\n", "
\n", "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.\n", "
" ] }, { "cell_type": "markdown", "id": "4d0c3678b72cf1b7", "metadata": {}, "source": [ "## Additional options\n", "\n", "### Choosing the right metrics\n", "\n", "The available metrics depend on ``prediction_type``:\n", "\n", "| Metric | Edgewise | Nodewise | Notes |\n", "| --- | :---: | :---: | --- |\n", "| ``MUE`` | ✓ | ✓ | Mean unsigned error |\n", "| ``RMSE`` | ✓ | ✓ | Root mean squared error |\n", "| ``RAE`` | ✓ | ✓ | Relative absolute error vs. a naïve mean predictor |\n", "| ``R2`` | — | ✓ | R² (Pearson r²); not meaningful for relative ΔΔG data |\n", "| ``rho`` | — | ✓ | Pearson r |\n", "| ``KTAU`` | — | ✓ | Kendall's τ |\n", "| ``PI`` | — | ✓ | Predictive index (Pearlman *et al.*) |\n", "\n", "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.\n", "\n", "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.\n", "\n", "### Adjusting the confidence level\n", "\n", "The ``confidence_level`` parameter (default ``0.95``) controls the width of the reported CIs. Pass a different value to tighten or loosen the intervals:" ] }, { "cell_type": "code", "execution_count": 8, "id": "e5f167ac67af2e25", "metadata": { "ExecuteTime": { "end_time": "2026-05-19T09:21:50.857180Z", "start_time": "2026-05-19T09:21:49.212873Z" } }, "outputs": [], "source": [ "summary_df, comparison_df = compare_and_rank_results(\n", " femap=femap,\n", " rank_metric=\"MUE\",\n", " metrics_to_compute=[\"MUE\", \"RMSE\"],\n", " confidence_level=0.90, # 90 % CIs — narrower, easier to detect differences\n", ")" ] }, { "cell_type": "markdown", "id": "3265c28aba157f15", "metadata": {}, "source": [ "## Recap\n", "\n", "- ``compare_and_rank_results`` performs **joint bootstrapping** across all sources, ensuring paired, fair comparisons on the same resampled data.\n", "- 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).\n", "- The **CLD** encodes statistical groupings: sources sharing a letter are *not* significantly different; sources with no shared letter *are* significantly different.\n", "- With **more than two sources**, Holm-corrected p-values appear automatically in ``p-value corrected`` and drive the ``significant`` flags.\n", "- Use ``prediction_type=\"nodewise\"`` only for independent per-ligand ΔG estimates (e.g. ABFE), *not* for MLE-derived values from a relative network.\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "openfe_dev", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }