library(lmebayes)

ch09 <- readRDS(system.file(
  "extdata", "Chapter09_airbnb_poisson.rds", package = "lmebayes"
))

Chapter 09: Poisson GLMMs

Everything so far has been Gaussian. This chapter switches to counts, which changes the sampler substantially even though the interface barely moves: replace lmerb() with glmerb() and add family = poisson().

The data is airbnb from bayesrules (Johnson et al. 2022): listings in Chicago neighborhoods, with the number of reviews as the count outcome. The fit takes about half a minute because every Block 1 draw runs an envelope accept-reject sampler, so it is precomputed by data-raw/make_vignette_outputs.R.

What changes for non-Gaussian responses

	Gaussian	Poisson
Block 1 draw	Conjugate normal, closed form	Envelope accept-reject
Joint posterior	Multivariate normal (when $\tau^2$ fixed)	Never normal
Exact iid route	Available	Not available
Convergence rate	Certified by theorem	Local-Gaussian heuristic
Pilot stage	Only when $\tau^2$ sampled	Always

The last three rows matter. For a Gaussian target the rate $\lambda^*$ of Chapter 06 is a theorem; for a Poisson likelihood it is evaluated at IRLS-style curvature weights taken at the posterior mode, which makes it a good practical guide rather than a certificate. The engine compensates by always running a pilot stage. The distinction is spelled out in vignette("Chapter-B02", package = "lmebayesCore"), under “Beyond the Gaussian case”.

Log-concavity is what makes the Block 1 draw possible at all. The Poisson log-likelihood with a log link is log-concave in the linear predictor, so the likelihood-subgradient envelope of Nygren and Nygren (2006) is a valid bound and accept-reject sampling is exact — not a Metropolis approximation.

The model

data(airbnb, package = "bayesrules")
ab <- airbnb[, c("neighborhood", "reviews", "rating", "room_type")]
ab <- ab[complete.cases(ab), ]
ab$neighborhood <- droplevels(factor(ab$neighborhood))
ab$rating_c <- ab$rating - mean(ab$rating)

form <- reviews ~ rating_c + (1 + rating_c || neighborhood)
model_setup(form, data = ab, family = poisson())

#> Call:
#>   model_setup(formula = form, data = ab, family = poisson())
#> --- Level 1 (Within-Group) ---
#>   reviews ~ (Intercept) + rating_c
#>   Observations             : 1561
#>   Group-effect coefficients: 2
#>   Group                    : neighborhood  [43 levels]
#>   Full-rank D_j: 35 of 43 groups
#>     rank-deficient: Beverly, Brighton Park, Burnside, Edison Park, Englewood, O'Hare, West Englewood, West Lawn
#>   Full-rank & estimable: 35 of 35 full-rank (35 of 43 total neighborhood)
#> --- Level 2 (Across-Group) ---
#>   (Intercept) ~ 1
#>   rating_c    ~ 1
#> --- Level 2 Rank ---
#>   (Restricted to 35 estimable neighborhood)
#>   (Intercept)  groups=35   predictors=1   full-rank
#>   rating_c     groups=35   predictors=1   full-rank
#>   Rank-deficient flags:
#>     (Intercept)  FALSE
#>     rating_c     FALSE
#>   popef.rank_ok: TRUE  -- model rank looks OK

1561 listings across 43 neighborhoods. Each neighborhood gets its own baseline review rate and its own rating effect. Rating is centered so the intercept means “log expected reviews for an average-rated listing”.

Eight neighborhoods are rank-deficient. For a count model that usually means too few listings, or no variation in rating among them. They stay in the fit and still receive coefficients — with nothing to learn from in the deficient direction, they are shrunk to the population value — but they are set aside when the level-2 rank condition is checked and when priors are calibrated.

Priors and fitting

ps <- Prior_Setup_GLMM(form, data = ab, family = poisson(), pop.pwt = 0.01)

set.seed(2024)
fit <- glmerb(
  form,
  data         = ab,
  family       = poisson(),
  pfamily_list = pfamily_list(ps),
  n            = 200L,
  progbar      = FALSE
)
fit

#> Call:
#>   glmerb(formula = form, data = ab, family = poisson(), pfamily_list = pfamily_list(ps), 
#>     n = 200L, progbar = FALSE)
#> Bayesian generalized linear mixed model  [poisson; 200 draws, two-block Gibbs]
#> Formula: reviews ~ rating_c + (1 + rating_c || neighborhood) 
#> Random effects (variance components fixed at glmer estimates):
#>  Groups         Name        Std.Dev.
#>  neighborhood   (Intercept) 0.5135  
#>  neighborhood.1 rating_c    1.2326  
#> Number of obs: 1561,  groups: neighborhood, 43
#> --- Block 2 hyperparameters (ICM posterior mode; MCMC means when simulated) ---
#>   RE component  parameter        ICM mode   popef.means    draws SD
#>   ------------  -----------  ------------  ------------  ----------
#>   (Intercept)   (Intercept)        3.1713        3.1595      0.0864
#>   rating_c      (Intercept)        0.2312        0.2597      0.1927

Two interface differences from lmerb(). There is no dispersion_ranef argument — the Poisson family fixes the dispersion at one, so there is no $\sigma^2$ to specify. And the reference fit is glmer() rather than lmer(), which is what Prior_Setup_GLMM() calibrates against.

The column header reads ICM mode rather than fixef.mode: for a non-Gaussian likelihood the pre-sampling point estimate is the posterior mode, found by iterated conditional modes, and the mode is not the mean. That is precisely why a pilot stage is needed — a chain started at the mode begins at an unknown positive distance from the target mean, so the mean-error term of the convergence bound is no longer zero.

Interpreting on the count scale

Coefficients are on the log scale, so exponentiate to get rate ratios.

pm <- unlist(ch09$popef_means)
round(rbind(log_scale = pm, rate_ratio = exp(pm)), 3)
#>            (Intercept).(Intercept) rating_c.(Intercept)
#> log_scale                     3.16                0.260
#> rate_ratio                   23.56                1.297

The intercept exponentiates to the expected number of reviews for an average-rated listing in an average neighborhood — about 24. The rating_c coefficient is a multiplicative effect per unit of rating: each additional rating point multiplies expected reviews by roughly 1.3.

But look at the posterior standard deviation on that slope — it is nearly as large as the coefficient. The population-average rating effect is not well determined, and the summary explains why.

The full summary

summary(fit)

#> Call:
#>   glmerb(formula = form, data = ab, family = poisson(), pfamily_list = pfamily_list(ps), 
#>     n = 200L, progbar = FALSE)
#> Bayesian linear mixed model fit  [200 draws, two-block Gibbs]
#> Formula: reviews ~ rating_c + (1 + rating_c || neighborhood) 
#> Random effects (variance components fixed at glmer estimates):
#>  Groups         Name        Std.Dev.
#>  neighborhood   (Intercept) 0.5135  
#>  neighborhood.1 rating_c    1.2326  
#> Number of obs: 1561,  groups: neighborhood, 43
#> === Block 2 dispersion (RE variance tau^2_k) ===
#> Prior and glmer reference:
#>               Prior 1/E[1/tau2] E[tau2] sqrt(E[tau2]) disp_lower disp_upper
#> (Intercept) dNormal      0.2637  0.2637        0.5135         NA         NA
#> rating_c    dNormal      1.5194  1.5194        1.2326         NA         NA
#>              glmer glmer SD
#> (Intercept) 0.2637   0.5135
#> rating_c    1.5194   1.2326
#> Overview:
#>             Post.Mode Post.Mean Post.Sd MC Error Mean SD Pr(Prior_tail)
#> (Intercept)    0.2637    0.2637  0.0000   0.0000  0.5135             NA
#> rating_c       1.5194    1.5194  0.0000   0.0000  1.2326             NA
#>             Cand/draw
#> (Intercept)         1
#> rating_c            1
#> Distribution percentiles (tau^2):
#>               1.0%   2.5%   5.0% Median  95.0%  97.5% 99.0%
#> (Intercept) 0.2637 0.2637 0.2637 0.2637 0.2637 0.2637 0.264
#> rating_c    1.5194 1.5194 1.5194 1.5194 1.5194 1.5194 1.519
#> SD credible interval (sqrt(tau^2) draws; glmer SD for reference):
#>               2.5% Median  97.5% glmer SD
#> (Intercept) 0.5135 0.5135 0.5135    0.513
#> rating_c    1.2326 1.2326 1.2326    1.233
#> === Block 2: Level-2 fixed effects (hyperparameters) ===
#> Prior and glmer reference:
#>                          Prior Mean Prior.sd  glmer glmer.se
#> (Intercept)::(Intercept)     3.1946   0.8256 3.1645    0.083
#> rating_c::(Intercept)        0.0000   2.0697 0.2394    0.208
#> Overview:
#>                          fixef.mode fixef.means Post.Sd MC.Error Pr(Prior_tail)
#> (Intercept)::(Intercept)    3.17126     3.15954 0.08644    0.006          0.320
#> rating_c::(Intercept)       0.23124     0.25969 0.19266    0.014          0.085
#>                           
#> (Intercept)::(Intercept)  
#> rating_c::(Intercept)    .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Distribution percentiles:
#>                              1.0%     2.5%     5.0%   Median    95.0%    97.5%
#> (Intercept)::(Intercept)  2.99686  3.00845  3.02642  3.15813  3.30024  3.33337
#> rating_c::(Intercept)    -0.15794 -0.09201 -0.03112  0.26369  0.58579  0.68493
#>                          99.0%
#> (Intercept)::(Intercept) 3.379
#> rating_c::(Intercept)    0.736
#> === Block 1: Random effects (group-level) ===
#> Summary of posterior mode (ranef.mode) across groups:
#>                Mean      SD     Min      Q1  Median      Q3     Max MCMC.mean
#> (Intercept)  3.1711  0.4874  2.1303  2.8718  3.2571  3.3901  4.2406     3.164
#> rating_c     0.2331  1.1223 -1.8120 -0.6153  0.2560  0.9860  2.2906     0.243
#> Per-group random effects: inspect fit$groupef.mode or fit$groupef,
#>   or call summary(fit, groups = <level ids>) for selected groups.

The random-effects standard deviation for rating_c is larger than the population mean effect itself. Neighborhoods genuinely differ in how rating relates to review volume — in some, better-rated listings get many more reviews; in others the relationship is flat or reversed. Reporting only the population coefficient would describe none of them well. This is the case for a random slope rather than a fixed one, and it is visible directly in the variance components.

There is no Residual row in the dispersion table. Poisson has no free dispersion parameter; the variance is determined by the mean. If your counts are overdispersed relative to that, the random intercept absorbs some of it, but a genuine overdispersion model is outside this package’s model class.

Convergence

c(m_convergence = ch09$m_convergence)
#> m_convergence 
#>             3

3 sweeps per stored draw, and a pilot stage ran before the main chains. Because the rate is heuristic here rather than certified, the sensible check is empirical: the sweep history should show the cross-chain mean flat and the cross-chain variance settled by the last sweep, exactly as in Chapter 06. fit$sweep_history$main carries it.

Cost

The Gaussian fits in earlier chapters took a second or two. This one takes about thirty. The difference is entirely Block 1: every group, on every sweep, on every chain, runs an envelope construction and an accept-reject loop. fit$groupef.iters reports how many candidates were needed per accepted draw — around eight here, which is normal for a Poisson model and means the envelope is doing real work.

Practical consequences: develop with small n and raise it for the final run; prefer tv_tol at its default rather than tightening it (each extra sweep is a full pass over all groups and chains); and remember that glmbayes Chapter-A10 documents optional OpenCL offload for envelope construction, which this package inherits.

References

Johnson, Alicia A., Miles Q. Ott, and Mine Dogucu. 2022. Bayes Rules! An Introduction to Applied Bayesian Modeling. CRC Press. https://www.bayesrulesbook.com.

Nygren, K. N., and L. M. Nygren. 2006. “Likelihood Subgradient Densities.” Journal of the American Statistical Association 101 (475): 1144–56. https://doi.org/10.1198/016214506000000357.

Chapter 09: Poisson GLMMs

Chapter 09: Poisson GLMMs

What changes for non-Gaussian responses

The model

Priors and fitting

Interpreting on the count scale

The full summary

Convergence

Cost

See also

References