Purpose

This article explains how to use rtmb_glmer() for hierarchical models and generalized linear mixed models in BayesRTMB.

rtmb_glmer() is the main wrapper for models that combine:

• fixed effects written with R formulas;
• random effects such as (1 | id) and (x | id);
• non-Gaussian likelihoods such as Bernoulli, Poisson, negative binomial, ordered, and sequential models;
• Bayesian priors, MAP estimation, MCMC, variational inference, and classical estimation where applicable.

The function is designed to feel familiar to users of lme4::glmer(), while returning a BayesRTMB RTMB_Model object. From the same model object, you can call:

fit_mcmc <- mdl$sample()
fit_map  <- mdl$optimize()
fit_vb   <- mdl$variational()
fit_cl   <- mdl$classic()

BayesRTMB’s functional priority is MCMC, MAP, VB, and then classical estimation. For final Bayesian inference, use MCMC whenever feasible. MAP is useful for fast checks and Laplace approximation. VB is useful for approximate inference in larger models. Classical estimation is useful for flat-prior wrapper models when you want familiar frequentist summaries.

1. A Minimal Random-Intercept Model

The basic syntax is close to lme4.

library(BayesRTMB)
data(debate)

mdl <- rtmb_glmer(
  sat ~ talk + perf + (1 | group),
  data = debate,
  family = "gaussian",
  prior = prior_normal()
)

This code has:

• fixed effects for talk and perf;
• a random intercept for group;
• a Gaussian likelihood;
• normal and exponential priors generated by prior_normal().

At this point, the model has been created but not fitted.

fit_map <- mdl$optimize()
fit_map

## Starting RTMB optimization...
## 
## 
## Call:
## MAP Estimation via RTMB
## 
## Negative Log-Posterior: 403.89
## Approx. Log Marginal Likelihood (Laplace): -412.12
## Note: Random effects are stored in $random_effects
## 
## Point Estimates and 95% Wald CI:
##    variable  Estimate  Std. Error  Lower 95%  Upper 95% 
## Intercept     2.11240     0.25310    1.61650    2.60820 
## b[talk]       0.26980     0.05520    0.16160    0.37800 
## b[perf]       0.14850     0.03120    0.08730    0.20960 
## sigma         0.77090     0.03910    0.69780    0.85150 
## sd[Int]       0.51020     0.06780    0.39310    0.66210 

Random effects are stored separately in fit_map$random_effects.

2. Choosing an Inference Method

fit_mcmc <- mdl$sample(
  sampling = 1000,
  warmup = 1000,
  chains = 4
)

fit_mcmc$summary()

## variable     mean    sd      map     q2.5    q97.5  ess_bulk  ess_tail  rhat
## lp        -407.11  2.34  -405.02  -412.54  -403.91       812       963  1.00
## b[talk]      0.27  0.06     0.27     0.16     0.38      1210      1314  1.00
## b[perf]      0.15  0.03     0.15     0.09     0.21      1084      1172  1.00
## sigma        0.78  0.04     0.77     0.70     0.86       940      1021  1.00
## sd[Int]      0.53  0.08     0.51     0.39     0.70       731       884  1.01

fit_map <- mdl$optimize(laplace = TRUE)

fit_vb <- mdl$variational(
  method = "meanfield",
  iter = 5000
)

mdl_flat <- rtmb_glmer(
  sat ~ talk + perf + (1 | group),
  data = debate,
  family = "gaussian",
  prior = prior_flat()
)

fit_cl <- mdl_flat$classic()

3. Formula Syntax and Visualization

rtmb_glmer(y ~ x + (1 | id), data = dat, family = "gaussian")

rtmb_glmer(y ~ x + (x | id), data = dat, family = "gaussian")

rtmb_glmer(y ~ x + (1 | subject) + (1 | item),
  data = dat,
  family = "bernoulli"
)

mdl_int <- rtmb_glmer(
  sat ~ talk * perf + (1 | group),
  data = debate,
  family = "gaussian",
  prior = prior_normal()
)

fit_int <- mdl_int$sample()

ce <- conditional_effects(fit_int, effect = "talk:perf")
plot(ce)

## Conditional effects for: talk:perf
## focal variable: talk
## moderator: perf
## moderator values: mean - 1 SD, mean, mean + 1 SD

simple_effects(fit_int, effect = "talk:perf")

## --- Simple Effects Analysis ---
## moderator perf          term estimate  lower upper
## perf       2.930 Slope of talk    0.038 -0.118 0.191
## perf       4.690 Slope of talk    0.266  0.161 0.369
## perf       6.450 Slope of talk    0.494  0.354 0.631

4. Families

rtmb_glmer() supports several likelihood families.

Example:

mdl_bin <- rtmb_glmer(
  y ~ x + (1 | id),
  data = dat,
  family = "bernoulli",
  prior = prior_normal()
)

5. Data Handling

mdl_wide <- rtmb_glmer(
  cbind(y_t1, y_t2, y_t3) ~ cond + (1 | id),
  data = dat_wide,
  family = "gaussian",
  within = list(time = c("t1", "t2", "t3")),
  prior = prior_normal()
)

mdl_fac <- rtmb_glmer(
  y ~ cond * time + (1 | id),
  data = dat,
  family = "gaussian",
  factors = c("cond", "time"),
  prior = prior_normal()
)

mdl_cwc <- rtmb_glmer(
  y ~ x + x_cwc + (1 | group),
  data = dat,
  family = "gaussian",
  cwc = list(cluster = "group", pars = "x"),
  prior = prior_normal()
)

6. Priors

mdl_flat <- rtmb_glmer(
  y ~ x + (1 | id),
  data = dat,
  family = "gaussian",
  prior = prior_flat()
)

mdl_norm <- rtmb_glmer(
  y ~ x + (1 | id),
  data = dat,
  family = "gaussian",
  prior = prior_normal()
)

mdl_weak <- rtmb_glmer(
  y ~ x + (1 | id),
  data = dat,
  family = "gaussian",
  prior = prior_weak(y_range = c(1, 5))
)

mdl_rhs <- rtmb_glmer(
  y ~ x1 + x2 + x3 + (1 | id),
  data = dat,
  family = "gaussian",
  prior = prior_rhs(y_range = c(1, 5))
)

mdl_jzs <- rtmb_ttest(
  y ~ group,
  data = dat,
  prior = prior_jzs()
)

7. Residual Correlation

rtmb_glmer() can represent residual correlation structures such as AR(1).

mdl_ar1 <- rtmb_glmer(
  y ~ time + cond,
  data = dat,
  family = "gaussian",
  resid_corr = "ar1",
  resid_time = "time",
  resid_group = "id",
  prior = prior_normal()
)

Here, resid_time defines the ordering or lag variable within a residual-correlation block, and resid_group defines the independent blocks.

If the model has no random effect term such as (1 | id), you still need resid_group = "id" so that the AR(1) correlation is applied within each participant rather than across the whole data set.

If a random effect grouping factor is already present and resid_group is omitted, BayesRTMB may use the first grouping factor internally. However, it is usually clearer to specify resid_group explicitly.

Avoid adding (1 | id) reflexively when the purpose is only to model AR(1) residual dependence. A random intercept and a residual correlation structure can both be valid, but they answer different questions and can become redundant in simple repeated-measures settings.

8. ANOVA-Style Classical Workflows

When you want an ordinary analysis-of-variance workflow, use rtmb_lmer() or rtmb_glmer() with prior_flat() and then call classic().

mdl_aov <- rtmb_lmer(
  y ~ cond * time + (1 | id),
  data = dat,
  prior = prior_flat()
)

fit_aov <- mdl_aov$classic()
anova(fit_aov)

## Analysis of Variance Table
##           Df  Chisq Pr(>Chisq)
## cond       1  5.812     0.0159
## time       2 12.406     0.0020
## cond:time  2  6.134     0.0465

Estimated marginal means can be obtained with lsmeans().

emm <- lsmeans(fit_aov, specs = "cond")
emm
plot(emm)

##  cond estimate    SE lower.CL upper.CL
##  0       3.21 0.082     3.05     3.37
##  1       3.48 0.084     3.32     3.65

Use this route when your goal is close to conventional ANOVA, but you still want the model to be represented inside the BayesRTMB wrapper system.

9. Inspecting Generated Code

Wrappers are convenient, but they still create rtmb_code() internally.

mdl$print_code()

Use print_code() when you want to understand the exact likelihood, priors, random-effect declarations, transformed quantities, or generated quantities created by a wrapper.

10. Model Comparison

For Bayesian model comparison, fit with MCMC and use bridge sampling.

fit_full <- mdl$sample()
fit_full$bridgesampling()

## Bridge Sampling Converged: LogML = -412.36 (Error = 0.0184, ESS = 91.2)

Bayes factors can be calculated by comparing against models with parameters fixed.

bf <- fit_full$bayes_factor(fixed = list("b[talk]" = 0))
bf

## --- Bayes Factor Analysis (Bridge Sampling) ---
## Bayes Factor (BF12) : 18.42
## Log Bayes Factor    : 2.913
## Evidence            : Strong evidence for Model 1

Because Bayes factors depend strongly on priors, use carefully designed priors. prior_weak() is often a reasonable starting point when you want weak but not implausibly diffuse priors.

11. Related Articles

• Introduction: the overall BayesRTMB workflow.
• Quick Start: a minimal first model and basic plots.
• Wrapper Functions: overview of standard wrapper functions.
• Writing Model Codes: writing custom rtmb_code() models.
• RTMB Internals and Inference Algorithms: internal model representation and inference pipeline.