adding ppc_rootogram_grouped function

[behramulukir]

This PR adds ppc_rootogram_grouped function requested at #377. If you have any feedback, I am happy to hear @mhollanders.

ppc_rootogram_grouped simply creates multiple facets (plots) of the rootogram based on the grouping user passes to the function as an argument.

Something to note here is that in ppc_bars_grouped, force_axes_in_facets function is used to make sure all facets have axes. This was resulting in errors when used with ppc_rootogram_grouped(style="discrete") since -Inf results in NaN when used with scale_y_sqrt(). It was still used in other rootogram styles, but so far, I hadn't managed to get it working for discrete style. I'll mark this PR as ready once I find a solution to that.

Tasks:

• implement functions
• implement tests
• update documentation
• fixing axes in discrete style

Examples

y <- rpois(100, 20)
yrep <- matrix(rpois(10000, 20), ncol = 100)
group <- gl(2, 50, length = 100, labels = c("GroupA", "GroupB"))

ppc_rootogram_grouped(y, yrep, group, prob = 0.5)

ppc_rootogram_grouped(y, yrep, group, prob = 0.5, style = "discrete")

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.70%. Comparing base (a30a706) to head (654d22c).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #419      +/-   ##
==========================================
+ Coverage   98.66%   98.70%   +0.04%     
==========================================
  Files          35       35              
  Lines        5857     5889      +32     
==========================================
+ Hits         5779     5813      +34     
+ Misses         78       76       -2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mhollanders]

Fantastic, thank you! I will doubtlessly have a use case for this in the next two weeks, and I will report if there's any issues.

[mhollanders]

@behramulukir, sorry for the basic question, but how do I get to the function? I don't see it when I install the latest development version of bayesplot.

[behramulukir]

Did you install the correct branch @mhollanders? As the changes are not merged to the master branch yet, you should install the branch I am working on. You can do it with devtools::install_github("stan-dev/bayesplot", ref = "ppc-rootogram-grouped")

[mhollanders]

Thank you, and pardon for my ignorance. It works great!

[behramulukir]

No worries! I'm glad that you liked it

[jgabry]

Thanks @behramulukir and thanks @mhollanders for testing it out!

Something to note here is that in ppc_bars_grouped, force_axes_in_facets function is used to make sure all facets have axes. This was resulting in errors when used with ppc_rootogram_grouped(style="discrete") since -Inf results in NaN when used with scale_y_sqrt(). It was still used in other rootogram styles, but so far, I hadn't managed to get it working for discrete style. I'll mark this PR as ready once I find a solution to that.

Does 0 work instead of -Inf? If 0 works here we could add an argument to force_axes_in_facets() to change between -Inf and 0 (leaving the default as -Inf so it stays the same for the other plots).

[mhollanders]

It's tangentially related but I might as well piggy-back on to this post: are there any PPC tips for zero-inflated count data? This output is for occupancy models, where zero-inflation is modeled at the level of the site and species. So rare species will have lots of 0s. I'm attaching two plots, one for the aggregated (predicted) counts per species and site (Q), and one for the predicted counts per species, site, and replicate survey. (y) The excess zeros make them hard to interpret!

[behramulukir]

Does 0 work instead of -Inf? If 0 works here we could add an argument to force_axes_in_facets() to change between -Inf and 0 (leaving the default as -Inf so it stays the same for the other plots).

The issue here is that scale_y_sqrt adds a padding before the x-axis line so that it is plotted at the negative y range. This exists in the original ppc_rootogram(style="discrete") as well. It is normally fine and looks visually good. The issue comes up when we try to plot the y-axis line in the negative y range. Since we cannot assign a negative y value, we can only plot the y-axis line up to y=0. This creates a gap between the end of the y-axis line and the x-axis line. Trying different values in force_axes_in_facets() yielded me results that work only in some cases and fail in other cases. For example, here is a half-working solution I was able to come up with:

force_axes_in_facets_sqrt <- function() {
  thm <- bayesplot_theme_get()
  annotate(
    "segment",
    x = c(-Inf, -Inf), xend = c(0, -Inf),
    y = c(0, 0), yend = c(0, Inf),
    color = thm$axis.line$colour %||% thm$line$colour %||% "black",
    linewidth = thm$axis.line$linewidth %||% thm$line$linewidth %||% 0.5
  )
}

As you can see in the following examples, there is a gap between the y-axis line and the x-axis line, and it fails in cases where we have multiple rows, since it doesn't plot an x-axis.

ppc_rootogram_grouped(y, yrep, group, prob = 0.5, style = "discrete")
ppc_rootogram_grouped(y, yrep, group, prob = 0.5, style = "discrete", facet_args = list(nrow=2))

If we try to fix the issue of the lack of an x-axis line in multiple rows of plots, then we get two x-axis lines:

force_axes_in_facets_sqrt <- function() {
  thm <- bayesplot_theme_get()
  annotate(
    "segment",
    x = c(-Inf, -Inf), xend = c(Inf, -Inf),
    y = c(0, 0), yend = c(0, Inf),
    color = thm$axis.line$colour %||% thm$line$colour %||% "black",
    linewidth = thm$axis.line$linewidth %||% thm$line$linewidth %||% 0.5
  )
}

This is simply because we cannot have any negative value at y = c(0, 0), yend = c(0, Inf). The only viable solution I was able to find was removing the padding between the y=0 and the x-axis line by using scale_y_sqrt(expand = expansion(mult = c(0, 0.05))).

If we decide to use this solution, we need either to change the look of the discrete rootogram slightly or to have a small difference between the grouped and normal versions of the discrete rootogram. I am not sure which one is better, though.

[behramulukir]

It's tangentially related but I might as well piggy-back on to this post: are there any PPC tips for zero-inflated count data?

Among the existing PPC count plots, I think ppc_rootogram(style = "suspended") is worth trying if you haven't tried it already. Since it plots the differences between expected and observed counts, it might be easier to interpret where the observed count is zero, and the expected count is non-zero (or vice versa).

You can also check Visualizing Count Data Regressions Using Rootograms paper by Christian Kleiber and Achim Zeileis. The paper contains an example case where "a well-known data set from ethology, for which excess zeros and, more generally, overdispersion require treatment".

If you find a solution that you think is worth implementing in bayesplot, please let me know! I would be interested in implementing it (unless you don't want to do it yourself, of course)

[avehtari]

@mhollanders

I would highly recommend using Teemu's discrete rootogram. Using square scaling on both x and y axis improves readability in case of high proportions and counts, and discrete nature of the data and model is more emphasized.

Also for zero-inflation part you should use binary target calibration plots as illustrated, for example, in Nabiximols case study

[mhollanders]

Thanks! The rootograms in bayesplot are those rootograms, right?

[behramulukir]

Yes, ppc_rootogram(style = discrete) is an adaptation of Teemu's discrete rootogram to the general style of bayesplot.

[jgabry]

@behramulukir I had a chance to come back to this and I have two different ideas for handling the force_axes_in_facets issue. Let me know what you think.

1. First idea is to define this helper function:

force_axes_in_facets_sqrt <- function() {
  thm <- bayesplot_theme_get()
  annotate(
    "segment",
    x = c(-Inf, -Inf), xend = c(Inf, -Inf),
    y = c(0, 0), yend = c(0, Inf),
    color = thm$axis.line$colour %||% thm$line$colour %||% "black",
    linewidth = thm$axis.line$linewidth %||% thm$line$linewidth %||% 0.5
  )
}

Then inside the plotting function do something like this:

  if (style == "discrete") {
    y_scale <- if (!is.null(dots$group)) {
      scale_y_sqrt(expand = expansion(mult = c(0, 0.05)))
    } else {
      scale_y_sqrt()
    }
...

and then use y_scale instead of calling scale_y_sqrt() directly.

Then later:

if (style == "discrete") {
  g <- g + force_axes_in_facets_sqrt()
} else {
  g <- g + force_axes_in_facets()
}

2. Another possibility: In ggplot2 3.5.0+ (we currently require 3.4.0+ but we could increase), facet_wrap() can draw interior axes natively via axes / axis.labels, so we can maybe avoid the annotation hack entirely for this case. Here is some rough code:

bars_group_facets <- function(facet_args,
                              scales_default = "fixed",
                              force_axes = FALSE,
                              axis_labels_default = "all") {
  facet_args[["facets"]] <- "group"
  facet_args[["scales"]] <- facet_args[["scales"]] %||% scales_default
  
  if (force_axes) {
    facet_args[["axes"]] <- facet_args[["axes"]] %||% "all"
    facet_args[["axis.labels"]] <- facet_args[["axis.labels"]] %||%
      axis_labels_default
  }
  
  do.call("facet_wrap", facet_args)
}

# Then inside the plotting function: 
use_native_axes <- style == "discrete" && fixed_y(facet_args)
if (!use_native_axes) {
   g <- g + force_axes_in_facets()
}
g + bars_group_facets(
   facet_args,
   force_axes = use_native_axes,
   axis_labels_default = "margins"
)

One issue with this approach is we can't control the style of the interior axes, they just follow the general theme:

[behramulukir]

Thank you so much, @jgabry! I was also thinking about this and wasn't able to come up with anything new. It's great that you had fresh ideas.

I think your first idea is similar to what I proposed in the previous message, which didn't fully convince me.

The only viable solution I was able to find was removing the padding between the y=0 and the x-axis line by using scale_y_sqrt(expand = expansion(mult = c(0, 0.05))).

I personally prefer to have a gap between the dots and the x-axis line, and I felt like removing it led to worse results aesthetically and felt rather cluttered. It also broke the consistency of the graphs by having an unignorable visual difference between the grouped and non-grouped versions of the same plot.

However, I quite liked the second idea. I wasn't aware of the possibility of drawing interior axes natively, and I liked it. I think it looks nice in the image you shared. Compared to the interior axes of other grouped plots, it seems like the only difference is the existence of tick marks in the axis lines. I think that's a smaller difference than the gap between dots and the x-axis line, and the one which I would prefer.

What's your opinion, though? Do you have a preference between the two options?