Testing how accurate is parameter recovery for Poisson and binomial models for a large N=1000. Combinations of zeros (0 hits of 0 trials) do not pose a problem.

recovery_pois_binomial.R

library(tidyverse)
library(brms)
library(tidybayes)

n <- 1000
lambda <- 2
p <- 0.50

x <- rpois(n, lambda)
y <- rbinom(n, x, p)

d <- tibble(x, y)

mx <- brm(
  x ~ 1, data = d, family = poisson
)
summary(mx)
px <- as_draws_df(mx)
px$b_Intercept |> exp() |> mean_hdci()

#          y     ymin     ymax .width .point .interval
# 1 1.982956 1.895036 2.072006   0.95   mean      hdci

my <- brm(
  y | trials(x) ~ 1, data = d, family = binomial
)
summary(my)
py <- as_draws_df(my)
py$b_Intercept |> inv_logit_scaled() |> mean_hdci()
#           y      ymin      ymax .width .point .interval
# 1 0.4867739 0.4661131 0.5087191   0.95   mean      hdci