Skip to content

Take3 model averaging #414

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
Sep 21, 2022
Merged

Take3 model averaging #414

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
c75abb3
notebook: model_averaging (header,footer updates)
reshamas Aug 9, 2022
8451466
model_averaging: adding myst file
reshamas Aug 9, 2022
f4762e7
add watermark header; add footer
reshamas Aug 15, 2022
8c621a6
rerunning myst file
reshamas Aug 15, 2022
ff16430
update notebook tags
reshamas Aug 15, 2022
83b6f25
recreate myst file after updating notebook tags
reshamas Aug 15, 2022
71621b3
minor text fixes
reshamas Aug 15, 2022
34191d8
adding myst file
reshamas Aug 15, 2022
9dad6d6
add References section at the end
reshamas Aug 15, 2022
f191cd7
add myst file
reshamas Aug 15, 2022
dcf922b
add 2 more urls to list in pre-commit yaml file
reshamas Aug 17, 2022
942af75
update yaml, bibtext, citation references
reshamas Aug 21, 2022
5b538b2
adding myst file
reshamas Aug 21, 2022
e6cc033
remove end of |
reshamas Aug 21, 2022
7cd0054
revert unwanted updates to moderation_analysis + add excluded nb link
reshamas Sep 9, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
363 changes: 243 additions & 120 deletions examples/diagnostics_and_criticism/model_averaging.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

89 changes: 68 additions & 21 deletions myst_nbs/diagnostics_and_criticism/model_averaging.myst.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,20 @@ jupytext:
format_version: 0.13
jupytext_version: 1.13.7
kernelspec:
display_name: Python PyMC3 (Dev)
display_name: Python 3 (ipykernel)
language: python
name: pymc3-dev-py38
name: python3
---

(model_averaging)=
# Model Averaging

:::{post} Aug 2022
:tags: model comparison, model averaging
:category: intermediate
:author: Osvaldo Martin
:::

```{code-cell} ipython3
---
papermill:
Expand All @@ -27,7 +36,7 @@ import numpy as np
import pandas as pd
import pymc3 as pm

print(f"Running on PyMC3 v{pm.__version__}")
print(f"Running on PyMC3 v{pm.__version__}")
```

```{code-cell} ipython3
Expand All @@ -47,17 +56,16 @@ az.style.use("arviz-darkgrid")

+++ {"papermill": {"duration": 0.068882, "end_time": "2020-11-29T12:13:08.020372", "exception": false, "start_time": "2020-11-29T12:13:07.951490", "status": "completed"}, "tags": []}

# Model averaging

When confronted with more than one model we have several options. One of them is to perform model selection, using for example a given Information Criterion as exemplified [in this notebook](model_comparison.ipynb) and this other [example](GLM-model-selection.ipynb). Model selection is appealing for its simplicity, but we are discarding information about the uncertainty in our models. This is somehow similar to computing the full posterior and then just keep a point-estimate like the posterior mean; we may become overconfident of what we really know.
When confronted with more than one model we have several options. One of them is to perform model selection, using for example a given Information Criterion as exemplified the [model comparison notebook](https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/model_comparison.html) and the [GLM model selection notebook](https://www.pymc.io/projects/examples/en/latest/generalized_linear_models/GLM-model-selection.html). Model selection is appealing for its simplicity, but we are discarding information about the uncertainty in our models. This is somehow similar to computing the full posterior and then just keep a point-estimate like the posterior mean; we may become overconfident of what we really know.

One alternative is to perform model selection but discuss all the different models together with the computed values of a given Information Criterion. It is important to put all these numbers and tests in the context of our problem so that we and our audience can have a better feeling of the possible limitations and shortcomings of our methods. If you are in the academic world you can use this approach to add elements to the discussion section of a paper, presentation, thesis, and so on.

Yet another approach is to perform model averaging. The idea now is to generate a meta-model (and meta-predictions) using a weighted average of the models. There are several ways to do this and PyMC3 includes 3 of them that we are going to briefly discuss, you will find a more thorough explanation in the work by [Yuling Yao et. al.](https://arxiv.org/abs/1704.02030)
Yet another approach is to perform model averaging. The idea now is to generate a meta-model (and meta-predictions) using a weighted average of the models. There are several ways to do this and PyMC includes 3 of them that we are going to briefly discuss, you will find a more thorough explanation in the work by [Yuling Yao et. al](https://arxiv.org/abs/1704.02030).


## Pseudo Bayesian model averaging

Bayesian models can be weighted by their marginal likelihood, this is known as Bayesian Model Averaging. While this is theoretically appealing, is problematic in practice: on the one hand the marginal likelihood is highly sensible to the specification of the prior, in a way that parameter estimation is not, and on the other computing the marginal likelihood is usually a challenging task. An alternative route is to use the values of WAIC (Widely Applicable Information Criterion) or LOO (pareto-smoothed importance sampling Leave-One-Out cross-validation), which we will call generically IC, to estimate weights. We can do this by using the following formula:
Bayesian models can be weighted by their marginal likelihood, this is known as Bayesian Model Averaging. While this is theoretically appealing, it is problematic in practice: on the one hand the marginal likelihood is highly sensible to the specification of the prior, in a way that parameter estimation is not, and on the other, computing the marginal likelihood is usually a challenging task. An alternative route is to use the values of WAIC (Widely Applicable Information Criterion) or LOO (pareto-smoothed importance sampling Leave-One-Out cross-validation), which we will call generically IC, to estimate weights. We can do this by using the following formula:

$$w_i = \frac {e^{ - \frac{1}{2} dIC_i }} {\sum_j^M e^{ - \frac{1}{2} dIC_j }}$$

Expand All @@ -71,7 +79,7 @@ The above formula for computing weights is a very nice and simple approach, but

## Stacking

The third approach implemented in PyMC3 is know as _stacking of predictive distributions_ and it has been recently [proposed](https://arxiv.org/abs/1704.02030). We want to combine several models in a metamodel in order to minimize the diverge between the meta-model and the _true_ generating model, when using a logarithmic scoring rule this is equivalently to:
The third approach implemented in PyMC is known as [_stacking of predictive distributions_](https://arxiv.org/abs/1704.02030). We want to combine several models in a metamodel in order to minimize the divergence between the meta-model and the _true_ generating model, when using a logarithmic scoring rule this is equivalent to:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here (it's actually the same paper as above so so far 1 reference to add to bibtex)


$$\max_{w} \frac{1}{n} \sum_{i=1}^{n}log\sum_{k=1}^{K} w_k p(y_i|y_{-i}, M_k)$$

Expand All @@ -81,11 +89,11 @@ The quantity $p(y_i|y_{-i}, M_k)$ is the leave-one-out predictive distribution f

## Weighted posterior predictive samples

Once we have computed the weights, using any of the above 3 methods, we can use them to get a weighted posterior predictive samples. PyMC3 offers functions to perform these steps in a simple way, so let see them in action using an example.
Once we have computed the weights, using any of the above 3 methods, we can use them to get a weighted posterior predictive samples. PyMC offers functions to perform these steps in a simple way, so let see them in action using an example.

The following example is taken from the superb book [Statistical Rethinking](http://xcelab.net/rm/statistical-rethinking/) by Richard McElreath. You will find more PyMC3 examples from this book in this [repository](https://github.com/aloctavodia/Statistical-Rethinking-with-Python-and-PyMC3). We are going to explore a simplified version of it. Check the book for the whole example and a more thorough discussion of both, the biological motivation for this problem and a theoretical/practical discussion of using Information Criteria to compare, select and average models.
The following example is taken from the superb book [Statistical Rethinking](http://xcelab.net/rm/statistical-rethinking/) by Richard McElreath. You will find more PyMC examples from this book in this [repository](https://github.com/aloctavodia/Statistical-Rethinking-with-Python-and-PyMC3). We are going to explore a simplified version of it. Check the book for the whole example and a more thorough discussion of both, the biological motivation for this problem and a theoretical/practical discussion of using Information Criteria to compare, select and average models.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

statistical rethinking should be a citation (it already is in the bibtex file but it would be good to add the url to the bibtex citation). The link to the pymc port of the book code should point to https://github.com/pymc-devs/pymc-resources now instead.


Briefly, our problem is as follows: We want to explore the composition of milk across several primate species, it is hypothesized that females from species of primates with larger brains produce more _nutritious_ milk (loosely speaking this is done _in order to_ support the development of such big brains). This is an important question for evolutionary biologists and try to give and answer we will use 3 variables, two predictor variables: the proportion of neocortex compare to the total mass of the brain and the logarithm of the body mass of the mothers. And for predicted variable, the kilocalories per gram of milk. With these variables we are going to build 3 different linear models:
Briefly, our problem is as follows: We want to explore the composition of milk across several primate species, it is hypothesized that females from species of primates with larger brains produce more _nutritious_ milk (loosely speaking this is done _in order to_ support the development of such big brains). This is an important question for evolutionary biologists and try to give an answer we will use 3 variables, two predictor variables: the proportion of neocortex compare to the total mass of the brain and the logarithm of the body mass of the mothers. And for predicted variable, the kilocalories per gram of milk. With these variables we are going to build 3 different linear models:

1. A model using only the neocortex variable
2. A model using only the logarithm of the mass variable
Expand Down Expand Up @@ -211,7 +219,7 @@ az.plot_forest(traces, figsize=(10, 5));

+++ {"papermill": {"duration": 0.052958, "end_time": "2020-11-29T12:14:55.196722", "exception": false, "start_time": "2020-11-29T12:14:55.143764", "status": "completed"}, "tags": []}

Another option is to plot several traces in a same plot is to use `densityplot`. This plot is somehow similar to a forestplot, but we get truncated KDE plots (by default 95% credible intervals) grouped by variable names together with a point estimate (by default the mean).
Another option is to plot several traces in a same plot is to use `plot_density`. This plot is somehow similar to a forestplot, but we get truncated KDE (kernel density estimation) plots (by default 95% credible intervals) grouped by variable names together with a point estimate (by default the mean).

```{code-cell} ipython3
---
Expand All @@ -223,12 +231,25 @@ papermill:
status: completed
tags: []
---
az.plot_density(traces, var_names=["alpha", "sigma"]);
ax = az.plot_density(
traces,
var_names=["alpha", "sigma"],
shade=0.1,
data_labels=["Model 0 (neocortex)", "Model 1 (log_mass)", "Model 2 (neocortex+log_mass)"],
)

ax[0, 0].set_xlabel("Density")
ax[0, 0].set_ylabel("")
ax[0, 0].set_title("95% Credible Intervals: alpha")

ax[0, 1].set_xlabel("Density")
ax[0, 1].set_ylabel("")
ax[0, 1].set_title("95% Credible Intervals: sigma")
```

+++ {"papermill": {"duration": 0.055089, "end_time": "2020-11-29T12:14:57.977616", "exception": false, "start_time": "2020-11-29T12:14:57.922527", "status": "completed"}, "tags": []}

Now that we have sampled the posterior for the 3 models, we are going to use WAIC (Widely applicable information criterion) to compare the 3 models. We can do this using the `compare` function included with PyMC3.
Now that we have sampled the posterior for the 3 models, we are going to use WAIC (Widely applicable information criterion) to compare the 3 models. We can do this using the `compare` function included with PyMC.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

similar comment to the one in the introduction, compare is an arviz function now.

We should probably also add a note or comment on the code. I think the code does use waic, but az.compare now defaults to using loo instead, so running the same code will not use waic anymore.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

or maybe leave a note on the issue for whoever updates the code and reruns it on pymc v4? (which can't really be done for now as sample_posterior_predictive_w doesn't work on v4 yet)


```{code-cell} ipython3
---
Expand All @@ -247,11 +268,11 @@ comp

+++ {"papermill": {"duration": 0.056609, "end_time": "2020-11-29T12:14:58.387481", "exception": false, "start_time": "2020-11-29T12:14:58.330872", "status": "completed"}, "tags": []}

We can see that the best model is `model_2`, the one with both predictor variables. Notice the DataFrame is ordered from lowest to highest WAIC (_i.e_ from _better_ to _worst_ model). Check [this notebook](model_comparison.ipynb) for a more detailed discussing on model comparison.
We can see that the best model is `model_2`, the one with both predictor variables. Notice the DataFrame is ordered from lowest to highest WAIC (_i.e_ from _better_ to _worst_ model). Check [model_comparison notebook](model_comparison.ipynb) for a more detailed discussion on model comparison.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should also be a sphinx cross reference, not a markdown local link


We can also see that we get a column with the relative `weight` for each model (according to the first equation at the beginning of this notebook). This weights can be _vaguely_ interpreted as the probability that each model will make the correct predictions on future data. Of course this interpretation is conditional on the models used to compute the weights, if we add or remove models the weights will change. And also is dependent on the assumptions behind WAIC (or any other Information Criterion used). So try to do not overinterpret these `weights`.
We can also see that we get a column with the relative `weight` for each model (according to the first equation at the beginning of this notebook). This weights can be _vaguely_ interpreted as the probability that each model will make the correct predictions on future data. Of course this interpretation is conditional on the models used to compute the weights, if we add or remove models the weights will change. And also is dependent on the assumptions behind WAIC (or any other Information Criterion used). So try to not overinterpret these `weights`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This also needs a note, not sure how to rewrite it though, I'll try to come back later. The weight-probability interpretation is only valid for bma, not for stacking. The notebook should be clear on this because is is a common source of confusion, see arviz-devs/arviz#2077 or https://discourse.pymc.io/t/bayesian-model-averaging-ranking-of-model-weights-and-loo-dont-match/4658


Now we are going to use copmuted `weights` to generate predictions based not on a single model but on the weighted set of models. This is one way to perform model averaging. Using PyMC3 we can call the `sample_posterior_predictive_w` function as follows:
Now we are going to use computed `weights` to generate predictions based not on a single model, but on the weighted set of models. This is one way to perform model averaging. Using PyMC we can call the `sample_posterior_predictive_w` function as follows:

```{code-cell} ipython3
---
Expand All @@ -275,7 +296,7 @@ ppc_w = pm.sample_posterior_predictive_w(

Notice that we are passing the weights ordered by their index. We are doing this because we pass `traces` and `models` ordered from model 0 to 2, but the computed weights are ordered from lowest to highest WAIC (or equivalently from larger to lowest weight). In summary, we must be sure that we are correctly pairing the weights and models.

We are also going to compute PPCs for the lowest-WAIC model
We are also going to compute PPCs for the lowest-WAIC model.

```{code-cell} ipython3
---
Expand All @@ -292,7 +313,7 @@ ppc_2 = pm.sample_posterior_predictive(trace=trace_2, model=model_2, progressbar

+++ {"papermill": {"duration": 0.058214, "end_time": "2020-11-29T12:15:55.404271", "exception": false, "start_time": "2020-11-29T12:15:55.346057", "status": "completed"}, "tags": []}

A simple way to compare both kind of predictions is to plot their mean and hpd interval
A simple way to compare both kind of predictions is to plot their mean and hpd interval.

```{code-cell} ipython3
---
Expand Down Expand Up @@ -329,7 +350,30 @@ As we can see the mean value is almost the same for both predictions but the unc

There are other ways to average models such as, for example, explicitly building a meta-model that includes all the models we have. We then perform parameter inference while jumping between the models. One problem with this approach is that jumping between models could hamper the proper sampling of the posterior.

Besides averaging discrete models we can sometimes think of continuous versions of them. A toy example is to imagine that we have a coin and we want to estimated it's degree of bias, a number between 0 and 1 being 0.5 equal chance of head and tails. We could think of two separated models one with a prior biased towards heads and one towards tails. We could fit both separate models and then average them using, for example, IC-derived weights. An alternative, is to build a hierarchical model to estimate the prior distribution, instead of contemplating two discrete models we will be computing a continuous model that includes these the discrete ones as particular cases. Which approach is better? That depends on our concrete problem. Do we have good reasons to think about two discrete models, or is our problem better represented with a continuous bigger model?
Besides averaging discrete models we can sometimes think of continuous versions of them. A toy example is to imagine that we have a coin and we want to estimated its degree of bias, a number between 0 and 1 having a 0.5 equal chance of head and tails (fair coin). We could think of two separate models one with a prior biased towards heads and one towards tails. We could fit both separate models and then average them using, for example, IC-derived weights. An alternative, is to build a hierarchical model to estimate the prior distribution, instead of contemplating two discrete models we will be computing a continuous model that includes these the discrete ones as particular cases. Which approach is better? That depends on our concrete problem. Do we have good reasons to think about two discrete models, or is our problem better represented with a continuous bigger model?

+++

## Authors

* Authored by Osvaldo Martin in June 2017 ([pymc#2273](https://github.com/pymc-devs/pymc/pull/2273))
* Updated by Osvaldo Martin in December 2017 ([pymc#2741](https://github.com/pymc-devs/pymc/pull/2741))
* Updated by Marco Gorelli in November 2020 ([pymc#4271](https://github.com/pymc-devs/pymc/pull/4271))
* Moved from pymc to pymc-examples repo in December 2020 ([pymc-examples#8](https://github.com/pymc-devs/pymc-examples/pull/8))
* Updated by Raul Maldonado in February 2021 ([pymc#25](https://github.com/pymc-devs/pymc-examples/pull/25))
* Updated Markdown and styling by @reshamas in August 2022, ([pymc-examples#411](https://github.com/pymc-devs/pymc-examples/pull/411))

+++

## References

:::{bibliography}
:filter: docname in docnames
:::

+++

## Watermark

```{code-cell} ipython3
---
Expand All @@ -344,3 +388,6 @@ tags: []
%load_ext watermark
%watermark -n -u -v -iv -w
```

:::{include} ../page_footer.md
:::