code for converting counts into matrices of counts of counts for `bayesplot::ppc_rootgram`

[geryan]

library(greta)
#> 
#> Attaching package: 'greta'
#> The following objects are masked from 'package:stats':
#> 
#>     binomial, cov2cor, poisson
#> The following objects are masked from 'package:base':
#> 
#>     %*%, apply, backsolve, beta, chol2inv, colMeans, colSums, diag,
#>     eigen, forwardsolve, gamma, identity, rowMeans, rowSums, sweep,
#>     tapply
data("warpbreaks")

X <- as_data(model.matrix(breaks ~ wool + tension, warpbreaks))
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#> 
y <- as_data(warpbreaks$breaks)
#int <- variable()

int <- normal(mean = 0, sd = 2, truncation = c(0, Inf))
coefs <- normal(0, 5, dim = ncol(X) - 1)
beta <- c(int, coefs)

eta <- X %*% beta

distribution(y) <- poisson(exp(eta))


ycalc <- calculate(y, nsim = 10)

yrep <- ycalc$y[,,1] |>
  matrix(ncol = ncol(ycalc$y))

# m <- model(int, coefs)
# 
# opt(m)

library(bayesplot)
#> This is bayesplot version 1.13.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#>    * Does _not_ affect other ggplot2 plots
#>    * See ?bayesplot_theme_set for details on theme setting
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:greta':
#> 
#>     slice
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)

z <- as.numeric(y)
zmax <- max(z)

y_rootgram_tab <- tibble(value = z) |>
  count(value) |>
  full_join(
    y = tibble(value = 0:zmax)
  ) |>
  arrange(value) |>
  mutate(ycount = ifelse(is.na(n), 0, n)) |>
  select(-n)
#> Joining with `by = join_by(value)`


yrepmax <- max(yrep)

counttabs <- apply(
  X = yrep,
  MARGIN = 1,
  FUN = function(x, yrepmax) {
    
    tibble(value = x) |>
      count(value) |>
      full_join(
        y = tibble(value = 0:yrepmax),
        by = "value"
      ) |>
      arrange(value) |>
      mutate(n = ifelse(is.na(n), 0, n))
    
  },
  yrepmax = yrepmax
)

vlist <- lapply(
  X = counttabs,
  FUN = function(x){x$n}
)

yreptab <- do.call(rbind, vlist) |>
  apply(
    MARGIN = 2,
    FUN = mean
  ) |>
  tibble(
    yrep_mean = _,
    value = 0:yrepmax
  )


rg_tab <- full_join(
  y_rootgram_tab,
  yreptab,
  by = "value"
) |>
  mutate(ycount = ifelse(is.na(ycount), 0, ycount))


ppc_rootogram(
  y = rg_tab$ycount,
  yrep = rg_tab$yrep_mean,
  prob = 0.9,
  style = "hanging"
)
#> Error in validate_predictions(yrep, length(y)): is.matrix(predictions) is not TRUE


yrm <- do.call(rbind, vlist)
yrmt <- yrm[,1:(zmax+1)]

ppc_rootogram(
  y = y_rootgram_tab$ycount,
  yrep = yrmt,
  prob = 0.9,
  style = "suspended"
)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the bayesplot package.
#>   Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

ppc_rootogram(
  y = rg_tab$ycount,
  yrep = yrm,
  prob = 0.9,
  style = "suspended"
)

^{Created on 2025-10-17 with reprex v2.1.1}