Repository: IAMconsortium/pyam
Pull request #523 (2021-05-04 11:12:42)

Add feature to compute weights for de-biasing using count
# Please confirm that this PR has done the following:

- [x] Tests Added
- [x] Documentation Added
- ~Name of contributors Added to AUTHORS.rst~
- [x] Description in RELEASE_NOTES.md Added

# Description of PR

Based on a rekindling of the discussion on de-biasing of scenario ensembles by @Rlamboll, this is an initial implementation using the most basic approach to compute de-biasing weights: counting the number of model or scenario names. My hope is that having this initial PR can help to structure the discussion on provide a template to think about more advanced methods for debiasing (e.g., based on generic meta-indicators like projects, differences between certain variables, etc.)

Intended results
----------------

In the test, there is a simple IamDataFrame with three scenarios (i.e., scenario runs) from two models, where `model_a` implements scenarios (in the sense of a scenario protocol) `scen_a` and `scen_b` and `model_b` implements only `scen_a`. The method counts the number of occurrences by name on the `axis` (in the test case: "scenario"). The expected output is

|  model    | scenario | bias |
|----------|--------|-------|
|  model_a | scen_a |  0.5   |
|  model_b | scen_a |  0.5  |
|  model_a | scen_b |  1     |

Note: setting `axis=["model", "scenario"]` returns a bias-vector of 1 in a standard IAMC-template IamDataFrame.

Open questions
----------------

- ~Should the computed bias weights be inverted? i.e., `0.5, 0.5, 1`~
  -> Per discussion below: yes, changed from initial implementation
- Does the function signature make sense?
- The function directly adds the bias-weights to `df.meta` so that they can be retrieved with `df[<name>]`.
  An alternative implementation could be:
  The function `compute_bias()` returns the weights as pd.Series, so that a used would have to add it to `meta` explicitly
  e.g., `df.set_meta(meta=df.compute_bias(...), name="bias")`

danielhuppmann: cc @brutschki, @gidden
codecov[bot]: # [Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) Report
> Merging [#523](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) (71145ef) into [main](https://codecov.io/gh/IAMconsortium/pyam/commit/cf350e736a6d60001ea734c5215934a4eb3e388f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) (cf350e7) will **increase** coverage by `0.0%`.
> The diff coverage is `93.7%`.

[![Impacted file tree graph](https://codecov.io/gh/IAMconsortium/pyam/pull/523/graphs/tree.svg?width=650&height=150&src=pr&token=ThxWlNkonJ&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)

```diff
@@          Coverage Diff          @@
##            main    #523   +/-   ##
=====================================
  Coverage   93.1%   93.1%           
=====================================
  Files         43      45    +2     
  Lines       4945    4961   +16     
=====================================
+ Hits        4607    4622   +15     
- Misses       338     339    +1     
```


| [Impacted Files](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) | Coverage Δ | |
|---|---|---|
| [pyam/\_debiasing.py](https://codecov.io/gh/IAMconsortium/pyam/pull/523/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-cHlhbS9fZGViaWFzaW5nLnB5) | `83.3% <83.3%> (ø)` | |
| [pyam/core.py](https://codecov.io/gh/IAMconsortium/pyam/pull/523/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-cHlhbS9jb3JlLnB5) | `92.4% <100.0%> (+<0.1%)` | :arrow_up: |
| [tests/test\_feature\_debiasing.py](https://codecov.io/gh/IAMconsortium/pyam/pull/523/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-dGVzdHMvdGVzdF9mZWF0dXJlX2RlYmlhc2luZy5weQ==) | `100.0% <100.0%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium). Last update [cf350e7...71145ef](https://codecov.io/gh/IAMconsortium/pyam/pull/523?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium).

Rlamboll: In conversation with @gidden we have constructed something that would take a similar input here: [https://github.com/GranthamImperial/silicone/pull/127](https://github.com/GranthamImperial/silicone/pull/127) which used a dictionary and has the weighting factor rather than the number of scenarios - I feel the latter point is always better because it will be less confusing with non-integer biases, should one have reasons for weighting scenarios based on, say, historic performance accuracy or a continuous degree-of-overlap metric. 

I vaguely plan to allow DataFrame-style inputs at some point but the dictionary {(model, scenario): weight} makes more sense when only a few model/scenarios have weights other than 1 and people are hard-coding them, if there's a pyam tool to generate this table then maybe the increased readability wins out. 
danielhuppmann: Thanks for pointing to the PR, @Rlamboll - if I understand your statement and the code in silicone correctly, you think this function should return the inverse of the count? (i.e. the weighting factor).
Rlamboll: Yes, in this case you'd return weight=1/count. 
Rlamboll: Looks good to me! 

Patch:
diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md
index 166ff07..32c7d2b 100644
--- a/RELEASE_NOTES.md
+++ b/RELEASE_NOTES.md
@@ -1,6 +1,7 @@
 # Next Release
 
 - [#527](https://github.com/IAMconsortium/pyam/pull/527) Add an in-dataframe basic mathematical operations `subtract`, `add`, `multiply`, `divide`
+- [#523](https://github.com/IAMconsortium/pyam/pull/523) Add feature to compute weights for de-biasing using count
 - [#519](https://github.com/IAMconsortium/pyam/pull/519) Enable explicit `label` and fix for non-string items in plot legend 
 
 # Release v0.11.0
diff --git a/pyam/_debiasing.py b/pyam/_debiasing.py
new file mode 100644
index 0000000..7378c8e
--- /dev/null
+++ b/pyam/_debiasing.py
@@ -0,0 +1,9 @@
+def _compute_bias(df, name, method, axis):
+    """Internal implementation for computing bias weights"""
+    if method == "count":
+        # invert from the count to obtain the weighting factor
+        count = 1 / df.meta.groupby(axis).count().exclude
+        count.name = name
+        df.meta = df.meta.join(count, on=axis, how="outer")
+    else:
+        raise ValueError(f"Unknown method {method} for computing bias weights!")
diff --git a/pyam/core.py b/pyam/core.py
index 655da42..ccd05df 100755
--- a/pyam/core.py
+++ b/pyam/core.py
@@ -11,6 +11,8 @@ import pandas as pd
 from pathlib import Path
 from tempfile import TemporaryDirectory
 
+from pyam._debiasing import _compute_bias
+
 try:
     from datapackage import Package
 
@@ -958,6 +960,52 @@ class IamDataFrame(object):
                 self._exclude_on_fail(df)
             return df.reset_index()
 
+    def compute_bias(self, name, method, axis):
+        """Compute the bias weights and add to 'meta'
+
+        Parameters
+        ----------
+        name : str
+           Column name in the 'meta' dataframe
+        method : str
+            Method to compute the bias weights, see the notes
+        axis : str
+            Index dimensions on which to apply the `method`
+
+        Notes
+        -----
+
+        The following methods are implemented:
+
+        - "count": use the inverse of the number of scenarios grouped by `axis` names.
+
+          Using the following method on an IamDataFrame with three scenarios
+
+          .. code-block:: python
+
+              df.compute_bias(name="bias-weight", method="count", axis="scenario")
+
+          results in the following column to be added to *df.meta*:
+
+          .. list-table::
+             :header-rows: 1
+
+             * - model
+               - scenario
+               - bias-weight
+             * - model_a
+               - scen_a
+               - 0.5
+             * - model_a
+               - scen_b
+               - 1
+             * - model_b
+               - scen_a
+               - 0.5
+
+        """
+        _compute_bias(self, name, method, axis)
+
     def rename(
         self, mapping=None, inplace=False, append=False, check_duplicates=True, **kwargs
     ):
@@ -1064,7 +1112,7 @@ class IamDataFrame(object):
     def convert_unit(
         self, current, to, factor=None, registry=None, context=None, inplace=False
     ):
-        r"""Convert all data having *current* units to new units.
+        """Convert all timeseries data having *current* units to new units.
 
         If *factor* is given, existing values are multiplied by it, and the
         *to* units are assigned to the 'unit' column.


Test patch:
diff --git a/tests/test_feature_debiasing.py b/tests/test_feature_debiasing.py
new file mode 100644
index 0000000..590eac9
--- /dev/null
+++ b/tests/test_feature_debiasing.py
@@ -0,0 +1,25 @@
+from pyam import IamDataFrame
+import pytest
+from numpy.testing import assert_array_equal
+
+
+@pytest.mark.parametrize(
+    "axis, exp",
+    (["scenario", [0.5, 0.5, 1]], [["model", "scenario"], [1, 1, 1]]),
+)
+def test_debiasing_count(test_pd_df, axis, exp):
+    """Check computing bias weights counting the number of scenarios by scenario name"""
+
+    # modify the default test data to have three distinct scenarios
+    test_pd_df.loc[1, "model"] = "model_b"
+    df = IamDataFrame(test_pd_df)
+    df.compute_bias(method="count", name="bias", axis=axis)
+
+    assert_array_equal(df["bias"].values, exp)
+
+
+def test_debiasing_unknown_method(test_df_year):
+    """Check computing bias weights counting the number of scenarios by scenario name"""
+    msg = "Unknown method foo for computing bias weights!"
+    with pytest.raises(ValueError, match=msg):
+        test_df_year.compute_bias(method="foo", name="bias", axis="scenario")