Draft plot_bpv (bayesian p-value) #1222

aloctavodia · 2020-06-04T19:33:28Z

This is a new kind of posterior predictive check plot, focused toward p_values. This a super early draft (the code is still a mess). I have doubts about the arguments names, suggestions are more than welcome. What do you think of these examples. We could also add other like ecdf, but I prefer to concentrate first on getting these right.

BTW, the first too are the same as loo-pit, but using the data instead of the IS weights.

A few examples
kind="u_value", reference="analytical". We want a uniform distribution over the [0, 1] interval

kind="u_value", reference="samples" We want a uniform distribution over the [0, 1] interval

kind="p_value", reference="analytical", t_stat=None. We want a symmetric distribution centered at 0.5

kind="p_value", reference="samples", t_stat=None We want a symmetric distribution centered at 0.5

kind="t_stat", reference=None, t_stat="median", other values are "mean", "std", quantiles can be passed as string of number between 0 and 1. Finally the user can pass an arbitrary function. The density is the distribution of sampled T_stat. The dot represent the mean of the observed T_stat. The legend shows the bayesian p_value (bpv)

Follows official PR format
Includes a sample plot to visually illustrate the changes (only for plot-related functions)
New features are properly documented (with an example if appropriate)?
Includes new or updated tests to cover the new feature
Code style correct (follows pylint and black guidelines)
Changes are listed in changelog

aloctavodia · 2020-06-22T20:00:32Z

The basic elements are in place. So this is ready for review.

ColCarroll

Looks good to me! I didn't check the bokeh code at all, and most of the comments were on naming.

ColCarroll · 2020-06-25T13:06:33Z

CHANGELOG.md

@@ -10,6 +10,7 @@
 * `plot_trace` now supports multiple aesthetics to identify chain and variable
  shape and support matplotlib aliases (#1253)
 * `plot_hdi` can now take already computed HDI values (#1241)
+* `plot_bpv`. A new plot intented to plot Bayesian p-values (#1222)


Suggested change

* `plot_bpv`. A new plot intented to plot Bayesian p-values (#1222)

* `plot_bpv`. A new plot for Bayesian p-values (#1222)

ColCarroll · 2020-06-25T13:08:03Z

arviz/plots/bpvplot.py

+    data : az.InferenceData object
+        InferenceData object containing the observed and posterior/prior predictive data.
+    kind : str
+        Type of plot to display (u_value, p_value, t_stat). Defaults to u_value.


Can you add what u_value, p_value, t_stat are (or references)

ColCarroll · 2020-06-25T13:08:47Z

arviz/plots/bpvplot.py

+        acepted, see examples section for details.
+    bpv : bool
+        If True (default) add the bayesian p_value to the legend when kind = t_stat.
+    mean : bool


Not a strong feeling, but maybe change to plot_mean?

ColCarroll · 2020-06-25T13:09:24Z

arviz/plots/bpvplot.py

+        Whether or not to plot the mean T statistic. Defaults to True.
+    reference : str
+        How to compute the distributions used as reference for u_values or p_values. Allowed values
+        are "analytical" (default) and "samples". Use `None` to do not plot any reference.


Suggested change

are "analytical" (default) and "samples". Use `None` to do not plot any reference.

are "analytical" (default) and "samples". Use `None` to do not plot any reference. Defaults to "samples".

ColCarroll · 2020-06-25T13:09:33Z

arviz/plots/bpvplot.py

+        How to compute the distributions used as reference for u_values or p_values. Allowed values
+        are "analytical" (default) and "samples". Use `None` to do not plot any reference.
+    n_ref : int, optional
+        Number of reference distributions to sample when `reference=samples`


Suggested change

Number of reference distributions to sample when `reference=samples`

Number of reference distributions to sample when `reference=samples`. Defaults to 100.

ColCarroll · 2020-06-25T13:11:57Z

arviz/plots/plot_utils.py

+    """Check if value is a number between 0 and 1."""
+    try:
+        value = float(value)
+        return 0 < value < 1


I think you want this return after the except, since you are watching for a ValueError while casting to float, right?

I think this works ok.

This will work as written

ColCarroll · 2020-06-25T13:14:12Z

examples/bokeh/bokeh_plot_bpv_tstat.py

+import arviz as az
+
+data = az.load_arviz_data("regression1d")
+az.plot_bpv(data, kind="t_stat", t_stat="0.5", backend="bokeh")


Suggested change

az.plot_bpv(data, kind="t_stat", t_stat="0.5", backend="bokeh")

az.plot_bpv(data, kind="t_stat", t_stat=0.5, backend="bokeh")

Passing a string is actually valid

ColCarroll · 2020-06-25T13:14:23Z

examples/matplotlib/mpl_plot_bpv_tstat.py

+az.style.use("arviz-darkgrid")
+
+data = az.load_arviz_data("regression1d")
+az.plot_bpv(data, kind="t_stat", t_stat="0.5")


Suggested change

az.plot_bpv(data, kind="t_stat", t_stat="0.5")

az.plot_bpv(data, kind="t_stat", t_stat=0.5)

ColCarroll · 2020-06-25T13:14:50Z

examples/bokeh/bokeh_plot_bpv_tstat.py

@@ -0,0 +1,10 @@
+"""
+Bayesian p-value Posterior plot


Should this have a different name from the other, to indicate it is using a T stat?

ColCarroll · 2020-06-25T13:15:10Z

examples/matplotlib/mpl_plot_bpv_tstat.py

@@ -0,0 +1,15 @@
+"""
+Bayesian p-value Posterior plot


Same question on naming here

…ment to plot_mean, fix English

aloctavodia · 2020-06-25T19:26:51Z

maybe @ahartikainen wants to improve the bokeh code in a future PR ;-)

codecov · 2020-06-25T19:43:31Z

Codecov Report

Merging #1222 into master will decrease coverage by 0.13%.
The diff coverage is 81.99%.

@@            Coverage Diff             @@
##           master    #1222      +/-   ##
==========================================
- Coverage   93.21%   93.08%   -0.14%     
==========================================
  Files          98      101       +3     
  Lines        9734     9995     +261     
==========================================
+ Hits         9074     9304     +230     
- Misses        660      691      +31

Impacted Files	Coverage Δ
arviz/plots/bpvplot.py	`77.14% <77.14%> (ø)`
arviz/plots/backends/matplotlib/bpvplot.py	`82.14% <82.14%> (ø)`
arviz/plots/backends/bokeh/bpvplot.py	`84.26% <84.26%> (ø)`
arviz/plots/plot_utils.py	`94.20% <87.50%> (-0.42%)`	⬇️
arviz/plots/__init__.py	`100.00% <100.00%> (ø)`
arviz/plots/backends/matplotlib/__init__.py	`97.05% <100.00%> (+0.08%)`	⬆️
arviz/plots/backends/matplotlib/posteriorplot.py	`98.09% <0.00%> (ø)`
... and 1 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e5e4eab...73c89b6. Read the comment docs.

OriolAbril

Looks great and very useful, looking forward to the section on EABM too :)

There are many documentation nits that do not need to be addressed straight away, they could well be recommendations for future functions to slowly move towards numpydoc best practices and recommendations. On noting the default for example, I have noticed pandas uses both:

engine : {‘c’, ‘python’}, optional
Parser engine to use. The C engine is faster while the python engine is currently more feature-complete.

where the user has to understand that c is the default. And

compression : {‘infer’, ‘gzip’, ‘bz2’, ‘zip’, ‘xz’, None}, default ‘infer’

where the user has to understand that having a default means its optional.

OriolAbril · 2020-06-29T14:12:17Z

arviz/plots/bpvplot.py

+    ----------
+    data : az.InferenceData object
+        InferenceData object containing the observed and posterior/prior predictive data.
+    kind : str


Nit: I would modify this following numpydoc advise:

When a parameter can only assume one of a fixed set of values, those values can be listed in braces, with the default appearing first:

order : {'C', 'F', 'A'} Description of `order`.

OriolAbril · 2020-06-29T14:12:45Z

arviz/plots/bpvplot.py

+        If True (default) add the bayesian p_value to the legend when kind = t_stat.
+    plot_mean : bool
+        Whether or not to plot the mean T statistic. Defaults to True.
+    reference : str


same as above

OriolAbril · 2020-06-29T14:13:18Z

arviz/plots/bpvplot.py

+        computing u_values. Should be in the interval (0, 1]. Defaults to
+        0.94.


defaults to stats.hdi_prob rcParam

OriolAbril · 2020-06-29T14:14:12Z

arviz/plots/bpvplot.py

+        For "u_value" we compute pi := p(yi* ≤ yi | y). i.e. like a p_value but per observation yi.
+        This is also known as marginal p_value. The ideal distribution is uniform. This is similar
+        to the LOO-pit calculation/plot, the difference is than in LOO-pit plot we compute
+        pi = p(yi* r ≤ yi | y-i ), where y-i, is all other data except yi.


I think there is one extra r before the lower or equal than sign

OriolAbril · 2020-06-29T14:15:37Z

arviz/plots/bpvplot.py

+    backend_kwargs : bool, optional
+        These are kwargs specific to the backend being used. For additional documentation
+        check the plotting method of the backend.
+    group : {"prior", "posterior"}, optional


along the idea of comments above, this should be inverted to show the default first

OriolAbril · 2020-06-29T14:15:54Z

arviz/plots/bpvplot.py

+    ax : numpy array-like of matplotlib axes or bokeh figures, optional
+        A 2D array of locations into which to plot the densities. If not supplied, Arviz will create
+        its own array of plot areas (and return it).
+    backend : str, optional


same as above

OriolAbril · 2020-06-29T14:25:13Z

arviz/plots/bpvplot.py

@@ -0,0 +1,294 @@
+"""Bayesian p-value Posterior/Prior predictive plot."""
+import numpy as np
+from matplotlib.colors import to_hex


we can use az.plots.plot_utils.vectorized_to_hex. It basically calls matplotlib.colors.to_hex but works with lists of colors too (not relevant here) and we could at some point (far future) make a version of the function not requiring matplotlib, i.e. try importing from matplotlib and if not possible use bokeh version (as far as I know bokeh has no function that does this yet)

OriolAbril · 2020-06-29T14:31:36Z

arviz/plots/plot_utils.py

+    for idx in range(shape[1]):
+        density, xmin, xmax = _fast_kde(dist_rvs[:, idx])
+        x_s = np.linspace(xmin, xmax, len(density))
+        x_ss.append(x_s)
+        densities.append(density)


we should check memory performance and see if its possible to preallocate. Not sure how it will perform, but specially the append looks suspicious.

my memory is foggy but I think that was my first attempt, but then I changed because the size of the grid used by _fast_kde is not fixed and this make cause trouble. That should be fixed with the new upcoming kde. Nevertheless I agree this is something we need to revisit.

OriolAbril · 2020-06-29T14:37:52Z

arviz/tests/base_tests/test_plots_bokeh.py

@@ -966,3 +967,18 @@ def test_plot_rank(models, kwargs):
 def test_plot_dist_comparison_warn(models):
    with pytest.raises(NotImplementedError, match="The bokeh backend.+Use matplotlib bakend."):
        plot_dist_comparison(models.model_1, backend="bokeh")
+
+


I think everything will work with multidimensional data (that is chain, draw, *shape with shape being multidimensional) thanks to the skip_dims and careful reshaping, however, it would be great to test on multidimensional models too. Should be the same as test below but using multidim_models fixture instead, probably a couple cases are enough for multidim.

OriolAbril · 2020-06-29T14:44:36Z

arviz/plots/bpvplot.py

+        pi = p(yi* r ≤ yi | y-i ), where y-i, is all other data except yi.
+        For "t_stat" we compute := p(T(y)* ≤ T(y) | y) where T is any T statistic. See t_stat
+        argument below for details of available options.
+    t_stat : str, float, or callable


How about adding a "identity" option too to encourage users to calculate T statistics manually and store them in idata object and then use these new variables as variable names with t_stat="identity"?

I know it is already possible to do using t_stat=lambda x: x, hence the encourage above. This would be helpful if T statistic function was expensive to calculate or for convenience to use T statistic with both plot_bpv and loo_pit

aloctavodia · 2020-06-29T15:22:03Z

Thanks for the comments @OriolAbril!

aloctavodia added 2 commits June 4, 2020 16:11

draft new plot

f6b112d

add matplotlib backend

0283beb

aloctavodia changed the title ~~[WIP] draft new plot (t_stat_ppc, p_values)~~ [WIP] draft plot_bpv (bayesian p-value) Jun 8, 2020

aloctavodia added 5 commits June 8, 2020 17:21

fix doc style

40f8931

clean arguments

48066d3

add bokeh plot

5a27a6f

fix conflict plot_utils

b021105

add examples

8089468

aloctavodia force-pushed the plot_bpv branch from 0765d73 to 8089468 Compare June 22, 2020 19:33

aloctavodia added 2 commits June 22, 2020 16:35

Merge branch 'master' into plot_bpv

d3a55b4

add tests

1066491

aloctavodia changed the title ~~[WIP] draft plot_bpv (bayesian p-value)~~ Draft plot_bpv (bayesian p-value) Jun 22, 2020

aloctavodia added 2 commits June 23, 2020 09:54

remove logging fix doc style

59c5f6b

blackify

d06cd5e

ColCarroll approved these changes Jun 25, 2020

View reviewed changes

improve docstring, explain option and add reference, rename mean argu…

73c89b6

…ment to plot_mean, fix English

canyon289 approved these changes Jun 26, 2020

View reviewed changes

aloctavodia merged commit 8099e93 into master Jun 29, 2020

aloctavodia deleted the plot_bpv branch June 29, 2020 11:33

OriolAbril reviewed Jun 29, 2020

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Draft plot_bpv (bayesian p-value) #1222

Draft plot_bpv (bayesian p-value) #1222

aloctavodia commented Jun 4, 2020 •

edited

aloctavodia commented Jun 22, 2020

ColCarroll left a comment

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

aloctavodia Jun 25, 2020

canyon289 Jun 26, 2020

ColCarroll Jun 25, 2020

aloctavodia Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

ColCarroll Jun 25, 2020

aloctavodia commented Jun 25, 2020

codecov bot commented Jun 25, 2020

OriolAbril left a comment

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

aloctavodia Jun 29, 2020

OriolAbril Jun 29, 2020

OriolAbril Jun 29, 2020

aloctavodia commented Jun 29, 2020

	* `plot_bpv`. A new plot intented to plot Bayesian p-values (#1222)
	* `plot_bpv`. A new plot for Bayesian p-values (#1222)

	are "analytical" (default) and "samples". Use `None` to do not plot any reference.
	are "analytical" (default) and "samples". Use `None` to do not plot any reference. Defaults to "samples".

	Number of reference distributions to sample when `reference=samples`
	Number of reference distributions to sample when `reference=samples`. Defaults to 100.

	az.plot_bpv(data, kind="t_stat", t_stat="0.5", backend="bokeh")
	az.plot_bpv(data, kind="t_stat", t_stat=0.5, backend="bokeh")

		computing u_values. Should be in the interval (0, 1]. Defaults to
		0.94.

Draft plot_bpv (bayesian p-value) #1222

Draft plot_bpv (bayesian p-value) #1222

Conversation

aloctavodia commented Jun 4, 2020 • edited

aloctavodia commented Jun 22, 2020

ColCarroll left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aloctavodia commented Jun 25, 2020

codecov bot commented Jun 25, 2020

Codecov Report

OriolAbril left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

aloctavodia commented Jun 29, 2020

aloctavodia commented Jun 4, 2020 •

edited