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.67%. Comparing base (042663b) to head (8f8ee2c).

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #419   +/-   ##
=======================================
  Coverage   98.66%   98.67%           
=======================================
  Files          35       35           
  Lines        5864     5904   +40     
=======================================
+ Hits         5786     5826   +40     
  Misses         78       78           

☔ 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?

[jgabry]

You're right my first option is basically one of the options you already mentioned. If we like that option but don't like that the points at y=0 are cut off (that was the problem right? was there something else too?) then we could potentially nudge those points up a little bit by adding a small constant to those:

What do you think about that? I think it looks good superficially, but is it maybe too confusing that the points aren't exactly at 0? Actually I think it should be possible to shift the y=0 tick mark and label up by the nudge amount to make it look like the points are exactly at 0! It would definitely be hack but it could work.

For option 2, I think it works great and is not a hack. We just can't make the style of the added axes match the style of other similar plots that use force_axes_in_facets().

So I guess we have to compromise either way.

[behramulukir]

I think it looks good superficially, but is it maybe too confusing that the points aren't exactly at 0? Actually, I think it should be possible to shift the y=0 tick mark and label up by the nudge amount to make it look like the points are exactly at 0! It would definitely be hack, but it could work.

I believe it would look quite confusing without the little hack you suggested. However, since using that hack would solve the issue, I am positive about this workaround. That is especially if you feel like it wouldn't result in undesired consequences in edge cases.

For option 2, I think it works great and is not a hack. We just can't make the style of the added axes match the style of other similar plots that use force_axes_in_facets().

To be honest, I am more positive about this solution. The main question I have is what changes do we apply to the style of the added axes compared to the general theme?

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

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

Looking at the code above, it seems like the force_axes_in_facets is getting the colour and linewidth of the added axes from the general theme, right? If that's the case, the added axes following the style of the general theme shouldn't be an issue. Or am I missing something?

[jgabry]

Looking at the code above, it seems like the force_axes_in_facets is getting the colour and linewidth of the added axes from the general theme, right? If that's the case, the added axes following the style of the general theme shouldn't be an issue. Or am I missing something?

You're right, but when you look at the plots of the two options I put above the axes look different. The interior axis in option 1 is thinner than the main axis. I'm not entirely sure why, but this is what always happens with the interior axes when using the default bayesplot theme (maybe other themes too) with force_axes_in_facets.

If we go with option 2 then the new interior axes have tick marks, which is inconsistent with all of the other plots that do use force_axes_in_facets. I don't think there's a way to remove those. But that's not too big of a deal.

So I think if we accept that the interior axes for the facets will look a little different for this plot than the other plots then option 2 is preferable.

[behramulukir]

You're right, but when you look at the plots of the two options I put above the axes look different. The interior axis in option 1 is thinner than the main axis.

Yes, and that's why I was confused. Looking at the code, it seemed like there shouldn't be any difference, but it's visible that there is.

So I think if we accept that the interior axes for the facets will look a little different for this plot than the other plots then option 2 is preferable.

Personally, I am fine with this as the difference is rather small, and it may be possible in future to tackle this inconsistency. However, if you feel like the consistency here is critical, I think that's a fair argument. In that case, I am also okay with going with the little hack you suggested 2 comments ago.

[jgabry]

I’m ok with the small inconsistency too. I think it’s probably better not to use the hacky solution I came up with.

[behramulukir]

Cool! I just increased the ggplot version to 3.5.0+ and added your solution for ppc_rootogram_grouped. It should pass tests and be ready to get merged