Add LOO Difference Plot #178

ParadaCarleton · 2021-05-28T17:33:51Z

Implement plot from Bayesian workflow paper

ParadaCarleton · 2021-05-28T19:52:42Z

Closes #127

codecov-commenter · 2021-05-28T20:43:36Z

Codecov Report

Merging #178 (a0417dd) into master (0117858) will decrease coverage by 2.31%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##           master     #178      +/-   ##
==========================================
- Coverage   95.15%   92.84%   -2.32%     
==========================================
  Files          28       29       +1     
  Lines        2663     2697      +34     
==========================================
- Hits         2534     2504      -30     
- Misses        129      193      +64

Impacted Files	Coverage Δ
R/loo_difference_plot.R	`0.00% <0.00%> (ø)`
R/effective_sample_sizes.R	`88.18% <0.00%> (-9.45%)`	⬇️
R/importance_sampling.R	`87.23% <0.00%> (-5.32%)`	⬇️
R/loo.R	`93.82% <0.00%> (-4.53%)`	⬇️
R/loo_moment_matching.R	`98.42% <0.00%> (-0.79%)`	⬇️

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 0117858...a0417dd. Read the comment docs.

jgabry

Hey Carlos, thanks for working on this! There are a bunch of small things I put in review comments (most of it related to things that matter for package development but not necessarily when writing code outside of package development), but this is a really really great start. Once we sort out some of smaller issues I commented on we can then get some other people, e.g. @avehtari or the other loo package authors to also take a look and see what they think about the plot itself and the different options provided.

jgabry · 2021-05-28T20:29:38Z

R/loo_difference_plot.R

+           alpha = 1,
+           jitter = 0,
+           quantiles = FALSE,
+           sortByGroup = FALSE


I would use sort_by_group instead of SortByGroup to match the style already used in the package (and also in this function, e.g. psis_object_1)

jgabry · 2021-05-28T20:30:15Z

R/loo_difference_plot.R

+#' 
+#' # Plot using groups from WHO
+#' 
+#' plot_loo_dif(factor(GM@data$super_region_name), loo3, loo2, 


It looks like here in the doc it's plot_loo_dif but in the code it's plot_loo_variation.

jgabry · 2021-05-28T20:36:25Z

R/loo_difference_plot.R

+#' 
+#' # Plot using groups from WHO
+#' 
+#' plot_loo_dif(factor(GM@data$super_region_name), loo3, loo2, 


Right now none of these examples run because the data GM and loo objects (loo2, loo3) don't exist. For the Examples section in the doc we'll need to change this to a self contained example or add all this data to the loo package itself so it can be used in the example. One possibility is to use a toy example here in the doc and then potentially add a more real example (e.g. this one from the paper) in one of the package vignettes. I'll think a bit more about this too.

I believe this is the only problem that's left unresolved; as I mentioned in the email, I'd like to see if it's possible to add this data to the LOO package, since I think the example here is really great.

jgabry · 2021-05-28T20:37:24Z

R/loo_difference_plot.R

+           psis_object_2,
+           ...,
+           group = FALSE,
+           outlier_thresh = FALSE,


I don't think this argument has documentation

jgabry · 2021-05-28T20:40:38Z

R/loo_difference_plot.R

+      values <- group_by(tibble(group, elpdDif), factor(group)) %>%
+                    arrange(.by_group = TRUE)
+
+      elpdDif <- pull(values, elpdDif)
+      group <- pull(values, group)
+
+    }


It looks like this is using tibble and dplyr code but we should avoid that if possible. I like those packages in general but for package development we don't want to add two new dependencies to the loo package for functionality that can be accomplished without those dependencies, even if the code is a bit more cumbersome.

jgabry · 2021-05-28T20:42:23Z

R/loo_difference_plot.R

+    plot <- ggplot(mapping=aes(y, elpdDif)) +
+            geom_hline(yintercept=0) + 
+            xlab(ifelse(sortByGroup, "y", "Index")) +
+            ylab(expression(ELPD[i][1] - ELPD[i][2])) + 
+            labs(color = "Groups")


In order to call functions from ggplot (or any other package) inside of the loo package source code we need to prefix them with ggplot2::, e.g. ggplot2::geom_hline().

For example, here are some calls to the parallel package from a different file in the loo package:

loo/R/loo_moment_matching.R

Lines 130 to 133 in 0117858

cl <- parallel::makePSOCKcluster(cores)

on.exit(parallel::stopCluster(cl))

mm_list <- parallel::parLapply(cl = cl, X = I,

fun = function(i) loo_moment_match_i_fun(i))

We need to do the same thing with ggplot2.

jgabry · 2021-05-28T20:47:07Z

R/loo_difference_plot.R

+           psis_object_1,
+           psis_object_2,
+           ...,
+           group = FALSE,


I would use NULL instead of FALSE. To me a default of FALSE gives the impression that the possible values are TRUE and FALSE, but group=TRUE isn't a possible value. So I think NULL makes more sense as the default. Same with outlier_thresh. On the other hand, FALSE makes sense for quantiles since it seems to actually be an argument that needs a logical/boolean value.

avehtari · 2021-05-30T18:01:00Z

R/loo_difference_plot.R

+#' points. `jitter` can be either a number or a vector of numbers.
+#' Passing a single number will jitter variables along the x axis only, while 
+#' passing a vector will jitter along both axes.
+#' @param quantiles Boolean that determines whether to plot the quantiles of


For me quantiles sounds like it would accept vector of quantiles (real valued) as, e.g., dnorm and pnorm. I suggest that quantiles parameter would actually allow to define which quantiles are plotted and the type of the plot would be selected with argument called, e.g., type. With options type='y' (default) and type='quantiles', which would allow then also extending to other possible types- Can you also include example of how these plots look like in the discussion thread of this PR?

For me quantiles sounds like it would accept vector of quantiles (real valued) as, e.g., dnorm and pnorm.

This does seem like it might be confusing; I think it makes sense to remove this, since users can provide transformations of the values themselves by replacing y.

I suggest that quantiles parameter would actually allow to define which quantiles are plotted and the type of the plot would be selected with argument called, e.g., type. With options type='y' (default) and type='quantiles', which would allow then also extending to other possible types

I'm not sure what you mean here; could you elaborate?

Can you also include example of how these plots look like in the discussion thread of this PR?

Will do.

If users can provide transformations of the values themselves by replacing y, and this option is not needed then I think removing it as you did is the correct action and you don't need to care about the rest I said.

ParadaCarleton · 2021-05-31T06:44:24Z

@avehtari I've added an example of a plot using the IQ dataset from rstanarm:

The first model only has a constant intercept; the second model also includes a term for the mother's IQ. At the edges, the constant-based model is more accurate than the full model, while the IQ-based model is stronger in the center. This seems counterintuitive, until you realize that this reveals there is significant measurement error in IQ tests: extreme values are probably caused by mismeasurement and need to be adjusted for regression to the mean.

Other examples:

avehtari · 2021-05-31T08:03:30Z

The plots looks great! Thanks for making this PR!

ParadaCarleton · 2021-05-31T16:57:06Z

The plots looks great! Thanks for making this PR!

Thanks! Do you happen to know if the spatial data from the Bayesian visualization paper can be added to the loo package? I think it's a good example, but if it can't, I can remove it.

jgabry · 2021-06-07T18:09:18Z

Sorry for the delay in responding to this. Regarding the issue of including the data from the visualization paper, I think we should avoid adding it to the package but I think there's still a way we can still use it. Basically, we should use a different example in the Examples section in the documentation, but then we can add an example using the visualization paper data in one of the vignettes. This works because in the vignette we're allowed to download data that isn't included in the package. For example, we do this already in one of the vignettes with some other data:

loo/vignettes/loo2-with-rstan.Rmd

Lines 95 to 96 in 0117858

    
           url <- "http://stat.columbia.edu/~gelman/arm/examples/arsenic/wells.dat" 
        
           wells <- read.table(url)

So we could do something similar and download the visualization paper data from https://github.com/jgabry/bayes-vis-paper. @ParadaCarleton What do you think about that option?

ParadaCarleton · 2021-06-07T18:29:47Z

Sorry for the delay in responding to this. Regarding the issue of including the data from the visualization paper, I think we should avoid adding it to the package but I think there's still a way we can still use it. Basically, we should use a different example in the Examples section in the documentation, but then we can add an example using the visualization paper data in one of the vignettes. This works because in the vignette we're allowed to download data that isn't included in the package. For example, we do this already in one of the vignettes with some other data:

loo/vignettes/loo2-with-rstan.Rmd

Lines 95 to 96 in 0117858

url <- "http://stat.columbia.edu/~gelman/arm/examples/arsenic/wells.dat"

wells <- read.table(url)

So we could do something similar and download the visualization paper data from https://github.com/jgabry/bayes-vis-paper. @ParadaCarleton What do you think about that option?

Sounds good! I've removed the example. Where do you think I should put this -- should I add it to an existing vignette or create a new one?

jgabry · 2021-06-07T19:02:52Z

Where do you think I should put this -- should I add it to an existing vignette or create a new one?

Hmm, on the one hand I think it could be nice to have it in the main introduction vignette https://mc-stan.org/loo/articles/loo2-example.html. On the other hand, this would be using a totally different data set than the example in that vignette, so that's not ideal and maybe a new vignette is preferable.

I think if we go the route of a new vignette then it shouldn't just be about this one plot but rather about visualizing loo output in general (e.g., this plot, the Pareto k diagnostic plot, perhaps some of the loo related plots in bayesplot, etc.). But that would be more work than just adding this to an existing vignette so it depends how much you feel like working on this (no pressure!).

ParadaCarleton · 2021-08-14T01:14:19Z

@jgabry btw, I think we should merge this, unless there's some change you'd like to see. I might be able to get to building a vignette before classes start again, but given that I'm focused on adding new features to ParetoSmooth.jl I'm not sure I'll be able to. We can probably add this to some other vignette later.

avehtari · 2023-03-24T11:56:50Z

@ParadaCarleton and @jgabry , what is the status of this PR?

ParadaCarleton added 4 commits May 28, 2021 10:31

Add LOO variation plot

dca229b

typo

8786be3

Typo

8e329d0

Rename file

90ecb40

jgabry requested changes May 28, 2021

View reviewed changes

Made recommended changes

48b97b6

ParadaCarleton requested a review from jgabry May 28, 2021 22:08

avehtari reviewed May 30, 2021

View reviewed changes

Added continuous example, removed quantile option

5399cb8

ParadaCarleton added 2 commits May 31, 2021 09:54

Fixed error in example

4d085d6

Typo

d4b8641

Typo/formatting

a0417dd

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add LOO Difference Plot #178

Add LOO Difference Plot #178

ParadaCarleton commented May 28, 2021

ParadaCarleton commented May 28, 2021

codecov-commenter commented May 28, 2021 •

edited

Loading

jgabry left a comment

jgabry May 28, 2021

jgabry May 28, 2021

jgabry May 28, 2021

ParadaCarleton May 29, 2021

jgabry May 28, 2021

jgabry May 28, 2021 •

edited

Loading

jgabry May 28, 2021

jgabry May 28, 2021

jgabry May 28, 2021

avehtari May 30, 2021

ParadaCarleton May 30, 2021 •

edited

Loading

avehtari May 31, 2021

ParadaCarleton commented May 31, 2021

avehtari commented May 31, 2021

ParadaCarleton commented May 31, 2021

jgabry commented Jun 7, 2021

ParadaCarleton commented Jun 7, 2021 •

edited

Loading

jgabry commented Jun 7, 2021 •

edited

Loading

ParadaCarleton commented Aug 14, 2021

avehtari commented Mar 24, 2023

	cl <- parallel::makePSOCKcluster(cores)
	on.exit(parallel::stopCluster(cl))
	mm_list <- parallel::parLapply(cl = cl, X = I,
	fun = function(i) loo_moment_match_i_fun(i))

Add LOO Difference Plot #178

Are you sure you want to change the base?

Add LOO Difference Plot #178

Conversation

ParadaCarleton commented May 28, 2021

ParadaCarleton commented May 28, 2021

codecov-commenter commented May 28, 2021 • edited Loading

Codecov Report

jgabry 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

jgabry May 28, 2021 • edited Loading

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

ParadaCarleton May 30, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ParadaCarleton commented May 31, 2021

avehtari commented May 31, 2021

ParadaCarleton commented May 31, 2021

jgabry commented Jun 7, 2021

ParadaCarleton commented Jun 7, 2021 • edited Loading

jgabry commented Jun 7, 2021 • edited Loading

ParadaCarleton commented Aug 14, 2021

avehtari commented Mar 24, 2023

codecov-commenter commented May 28, 2021 •

edited

Loading

jgabry May 28, 2021 •

edited

Loading

ParadaCarleton May 30, 2021 •

edited

Loading

ParadaCarleton commented Jun 7, 2021 •

edited

Loading

jgabry commented Jun 7, 2021 •

edited

Loading