Boundary corrected KDE values and flag: ppc_loo_pit_overlay()

[ecoronado92]

@jgabry @avehtari Included a bc_kde optional flag to ppc_loo_pit_overlay for users to decide whether to compute the boundary corrected densities or the current stat_density ones based on #171

Flag is set to FALSE by default, and appropriate documentation has been added to function params and man stating the following

For `ppc_loo_pit_overlay()`, when set to TRUE function will
compute `"boundary corrected"` density values using beta kernel estimators.
Although preferable, the current implementation is slow (default = FALSE).

[codecov-commenter]

Codecov Report

Merging #235 into master will decrease coverage by 0.20%.
The diff coverage is 93.28%.

@@            Coverage Diff             @@
##           master     #235      +/-   ##
==========================================
- Coverage   98.55%   98.34%   -0.21%     
==========================================
  Files          32       32              
  Lines        4009     4107      +98     
==========================================
+ Hits         3951     4039      +88     
- Misses         58       68      +10     

Impacted Files	Coverage Δ
R/ppc-loo.R	96.51% <93.28%> (-3.49%)	⬇️

---

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 5500d2e...6768b81. Read the comment docs.

[jgabry]

Awesome, thanks for this!! Will review as soon as possible.

[jgabry]

This is great, thanks again! I made a bunch of separate comments with small suggestions. In addition to those here are two other minor things:

1. I think the argument name should be something like boundary_correction instead of bc_kde. Even though it's not as concise, it's much clearer just from seeing the argument name what the purpose is.

2. Can you add a test? That would go in the file https://github.com/stan-dev/bayesplot/blob/master/tests/testthat/test-ppc-loo.R. You can use the quantities already calculated in that file (y, yrep, etc.) and just add a test making sure the new function runs without error.

[ecoronado92]

@jgabry Hey! Thanks for the comments. I've made some changes and added some comments to your points above.

I also tried using the https://github.com/stan-dev/bayesplot/blob/master/tests/testthat/test-ppc-loo.R and didn't get an error. Would you mind checking if the error persists with the new changes?

Main changes in past commit:

• Changed bc_kde flag to boundary_correction (made changes to documentation accordingly)

• Make sure bc_vals length consistency wouldn't be an issue by ensuring certain code is further down the line (i.e. when the xs are generated in bc_pvals)

• Changed strict logic in bc_pvals now to make sure that CDF values aren't greater than user defined xmax

[jgabry]

Thanks for the changes! They look good. Just one thing:

I also tried using the https://github.com/stan-dev/bayesplot/blob/master/tests/testthat/test-ppc-loo.R and didn't get an error. Would you mind checking if the error persists with the new changes?

Just to clarify, I wasn't getting an error running those tests, I was getting an error using the stan_glm model from those tests but setting bc_kde=TRUE. Now if I do that with boundary_correction=TRUE I don't get an error but I get a bunch of warnings of the form:

50: In if (valid_pvals > xmax) { ... :
  the condition has length > 1 and only the first element will be used

because if doesn't like that valid_pvals is a vector. I just pushed a test in ed455bc to catch this, which will currently fail because of those warnings.

[ecoronado92]

50: In if (valid_pvals > xmax) { ... :
  the condition has length > 1 and only the first element will be used

@jgabry Dumb mistake on my end, just pushed an any() there which should take care of the warning

[jgabry]

@avehtari Do you want to take a look and try this out before we merge it?

@ecoronado92 In anticipation of merging this I just pushed a commit to this branch adding an item to the NEWS file and adding you as a contributor in authors list in the DESCRIPTION file. I think this is pretty much ready to go, but let's also see what @avehtari says.

In the meantime one more question about the speed. I just tried it on another example:

fit <- stan_glm(switch ~ I(dist/100) + log(arsenic), data = wells) # wells comes with rstanarm
y <- get_y(fit)
yrep <- posterior_predict(fit)
lw <- weights(psis(-log_lik(fit), r_eff = NA, cores = 4))
system.time(
  p <- ppc_loo_pit_overlay(y, yrep, lw, boundary_correction = TRUE)
)

   user  system elapsed 
226.708   1.501 228.761 

This dataset has about 3,000 observations and it took a 228 seconds to generate the plot. Is that consistent with the timings you were seeing or is that slower than expected?

[ecoronado92]

This dataset has about 3,000 observations and it took a 228 seconds to generate the plot. Is that consistent with the timings you were seeing or is that slower than expected?

@jgabry hmm that does look slower than expected.

I ran the radon data set with 919 observations and 100 reference uniforms, it took 23 seconds (without the plotting part).

 user  system elapsed 
 22.337   0.223  23.286 

Then I just expanded it to have approx 3,000 obs using rep and that took (without the plotting part)

   user  system elapsed 
168.698   0.745 179.890 

The bottleneck was generating the reference lines but I'm wondering if plotting it might add extra time since we're not using stat_density() and instead rely on geom_line().

[avehtari]

I think this is too slow. I didn't expect that it could be this slow, but then I guess there is a reason why so many R packages include C/C++ code.

Options: 1) find another faster package that makes that makes sensible KDE in closed range, 2) write own C/C++ code, 3) drop KDE density estimates and wait for a moment for ECDF based plots. What do you think?

[ecoronado92]

@jgabry Sounds great! Thank you.

I don't think I can approve the changes, but feel free to merge them.

[ecoronado92]

Options: 1) find another faster package that makes that makes sensible KDE in closed range, 2) write own C/C++ code, 3) drop KDE density estimates and wait for a moment for ECDF based plots. What do you think?

@avehtari @jgabry I can take a stab at option 2) this weekend using Rcpp. I'd probably use the Eigen library and would focus on speeding up the bc_dunif function (i.e. iterating through each element of the CDF value vectors).

Then depending on the speed gain we can decide to continue via this path or wait until option 3) is available. Thoughts?

[avehtari]

Great if you have time to test Rcpp speed. No pressure, so if it seems it's taking too much of your time, just let us know. If it's fast, it would stay even after we would have 3

[codecov-io]

Codecov Report

Merging #235 into master will decrease coverage by 0.93%.
The diff coverage is 73.15%.

@@            Coverage Diff             @@
##           master     #235      +/-   ##
==========================================
- Coverage   98.55%   97.61%   -0.94%     
==========================================
  Files          32       32              
  Lines        4009     4107      +98     
==========================================
+ Hits         3951     4009      +58     
- Misses         58       98      +40     

Impacted Files	Coverage Δ
R/ppc-loo.R	86.06% <73.15%> (-13.94%)	⬇️

---

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 5500d2e...6a60492. Read the comment docs.

[ecoronado92]

@jgabry I just changed boundary correction to TRUE in functions and all documentation. I also made sure to pull from master before adding this commit to make sure it's up to date with any changes there.

[jgabry]

@jgabry I just changed boundary correction to TRUE in functions and all documentation. I also made sure to pull from master before adding this commit to make sure it's up to date with any changes there.

Thanks! I think this is ready to go, just one last thing that I hadn't noticed before (sorry!). Right now there's a comment in the code

  # TODO currently boundary correction generates NAs for values near ends of vector
  # after convolution . Important to note that # of NAs "varies" so na.omit()
  # won't work and will cause size inconsistencies
  # Current fix, use last non-null value at head and tail of vector as padding

Is there anything that still needs to be done here? If so, can you move it from a code comment to a new issue? We can merge this with the current fix, but we'd prefer to avoid using code comments to keep track of TODO items.

Otherwise I think that's it. Thanks again for working on this, it's a big improvement!!

[ecoronado92]

Is there anything that still needs to be done here? If so, can you move it from a code comment to a new issue? We can merge this with the current fix, but we'd prefer to avoid using code comments to keep track of TODO items.

@jgabry Yes, I wasn’t able to figure why this was happening so I’m happy to move this to a separate issue. I’ll make changes later today and will ping you

[jgabry]

Awesome, thanks!

[ecoronado92]

@jgabry Removed the TODO from the code and created a new issue pointing at where the solution should be addressed in the code.

Let me know if we need to do anything else :)

[avehtari]

@ecoronado92 It would be useful to mention in the help text that the boundary correction is using the "reflection method" (method is better word than trick here). Even better if there would also be an article reference. Did ArviZ have a reference?

@jgabry would you have time to re-create the bayes-vis paper LOO-PIT figures with the new version? If you don't have time I'll do it.

[ecoronado92]

@avehtari Just added reflection method to the documentation and code comments.

Did ArviZ have a reference?

Here's the only references I found related to their LOO PIT function and/or the reflection method.

References
    ----------
    * Gabry et al. (2017) see https://arxiv.org/abs/1709.01449
    * https://mc-stan.org/bayesplot/reference/PPC-loo.html
    * Gelman et al. BDA (2014) Section 6.3

Happy to add as well to the help text

[jgabry]

Even better if there would also be an article reference. Did ArviZ have a reference?

Yeah I agree it would be good to cite the reflection method but I don't see a reference for that in arviz either. Just the reference to the preprint of our paper, like @ecoronado92 found.

@jgabry would you have time to re-create the bayes-vis paper LOO-PIT figures with the new version? If you don't have time I'll do it.

Yeah I have a few min right now so I will try to do that quickly.

[jgabry]

@avehtari Ok here are PIT plots from our paper using the new version:

One thing I notice is that we still have the axis limits cut off before 0 and 1, but now that we have the boundary correction we could let the axis limits include 0 and 1. What do you think? This would need to get changed in the ppc_loo_pit_overlay() source code but it's a one line change.

Also, just as a reference, here are the same plots with boundary_correction=FALSE:

[avehtari]

I asked in ArviZ slack for a reference

Thanks for the plots. They look better, and I think the interval should from 0 to 1.

[avehtari]

@tomicapretto and @aloctavodia helped to find a useful reference:
Boneva, L. I., Kendall, D., & Stefanov, I. (1971). Spline transformations: Three new diagnostic aids for the statistical data-analyst. Journal of the Royal Statistical Society. Series B (Methodological), 33(1), 1-71.
https://www.jstor.org/stable/2986005
See figs 5 and 6

[jgabry]

Ok great, I just pushed a commit citing that Boneva et al paper and I made the x-axis limits go out to 0 and 1 if boundary_correction=TRUE. Here's what the plots look like now:

@avehtari Anything else before we merge this?

[avehtari]

Looks good.

[jgabry]

Ok I think we're ready to merge this. Thanks again for implementing this!