Posterior predictive distributions for multilevel Poisson-GLM

General
1

Hi, all,

I have generated posterior predictive distributions (on the outcome scale) for a multilevel varying intercepts Poisson-GLM with five categories along the lines described in the Stan User’s Guide, Sec. 24.5.2 (24.5 Posterior prediction for regressions | Stan User’s Guide). Now, as a form of model validation, I’d like to use bayesplot or ggplot2 to overlay the observed Poisson-like distributions in the five categories with examples of the sampled posterior predictive distributions in each category, employing a combination of geom_bar(), facet_wrap() and an (hopefully existing) analog of geom_jitter() to shift the bars inside each plot slightly so they don’t overlap.

Can anyone point me to existing functions in bayesplot (I am not really making good progress with ppc_bars_grouped() and the likes as they require the observed predictor data as input rather than counterfactual data) or any good reference elsewhere to help me along with my objective? Any suggestions are much appreciated. Many thanks.

Cheers, Henk

1 Like
2

Alright, here’s what I figured out: possibly of interest to others, and surely with potential for improvement.

Cheers.

mod6samples <-
  as.matrix(x = model6Stan,
            pars = "ynew") %>%
  tibble::as_tibble(x = .)

qwe <- rep(x = rep(x = 1:5, each = 11), times = nrow(mod6samples))

yxc <- mod6samples %>%
  tidyr::pivot_longer(data = ., cols = everything()) %>%
  tibble::add_column(.data = .,
                     timeID = qwe,
                     .before = "name")

ggplot(data = yxc, mapping = aes(x = value)) +
  geom_bar(
    mapping = aes(y = ..prop.., fill = factor(name)),
    stat = "count",
    alpha = 0.5,
    width = 0.5,
    position = position_dodge2(
      width = 1,
      preserve = "single",
      padding = 0.5
    ),
    show.legend = FALSE
  ) +
  geom_bar(
    data = d,
    mapping = aes(x = count, y = ..prop.., group = timeID),
    stat = "count",
    fill = "black",
    width = 0.1
  ) +
  coord_cartesian(xlim = c(0, 8)) +
  labs(
    title = "title",
    subtitle  = "subtitle",
    caption = waiver(),
    tag = waiver()
  ) +
  xlab(label = "xlab") +
  ylab(label = "ylab") +
  scale_y_continuous(labels = scales::percent) +
  facet_wrap(facets = ~ timeID) +
  theme_bw()
1 Like
3

Cool that you manage to resolve this. If I understand your code and description correctly than no, this is not available anywhere out of the box, so implementing it yourself is the way to go.

Best of luck with your model!

1 Like