frmtmb for brms users

frmtmb reimplements a documented subset of the brms grammar with identical spelling, so brms model code ports mechanically: change brm() to frm(), keep or drop the priors, and get a fit in milliseconds-to-seconds instead of an MCMC run. Estimation is ML (or REML) with the Laplace approximation for latent effects; inference is frequentist (Wald and profile intervals, likelihood-ratio tests, AIC).

Installing frmtmb and its extensions

brms holds every family and term in one package. frmtmb holds the general ones in core and the ones for one domain in extension packages: for example cox() is in frmtmb.survival, car() in frmtmb.spatial and wiener() in frmtmb.eam. Without the package, such a spelling stops with a message that names it. frmtmb has no package that installs everything. To install frmtmb and every extension:

remotes::install_github("aforren1/frmtmb")
# every extension, after frmtmb itself; frmtmb.eam goes before
# frmtmb.learn, which imports it
for (p in c("eam", "latent", "ode", "survival", "spatial", "coupling",
            "sample", "learn")) {
  remotes::install_github("aforren1/frmtmb",
                          subdir = paste0("extensions/frmtmb.", p))
}

Then load the extensions a model needs, for example library(frmtmb.survival).

What ports directly

brms	frmtmb	notes
`bf(y ~ x + (x \| g))`	same	plus `us()`, `diag()`, `homdiag()`, `cs()`, `ar1()` wrappers
`bf(y ~ x, sigma ~ z + (1 \| g))`	same	any dpar, full grammar
`bf(y ~ x, sigma = 1)`	same	fixed via TMB `map`
`bf(y ~ x, sigma1 = "sigma2")`, `lf(sigma1 = "sigma2")`	same	one parameter, `sigma2`; `variables()` and `summary()` list `sigma1` with its value, as brms does; brms’s rules and messages; `theta1 = "theta2"` is refused, where brms stops with an internal error
`bf(y ~ 0 + g, cmc = FALSE)`, `lf(..., cmc = FALSE)`	same	treatment contrasts without the intercept column, in the population-level and the group-level terms, `mm()` included; refused on the covariance structures that read one coefficient per level (`ar1()`, `cs()`, `toep()`, `ou()`, …), which brms does not have
`y ~ .`	same	expanded against `data` before the formula is read, as in brms; each parameter formula on its own, so `sigma ~ .` includes the response; `update(fit, newdata = )` keeps the expanded formula, as in brms
`family = list(gaussian(), poisson())`	same	multivariate models, one family per response in formula order; an entry fills only a response that has no family
`update(fit, bf(~ ., family = acat()))`	same	a `bf()` whose formula is a delta: the location formula is updated, parameter formulas and constants are pooled, and its family wins
`s(x)`, `t2(x, z)`	same	mgcv smooths in any dpar; `te()`/`ti()` unsupported
`(1 \| p \| g)`	same	cross-formula RE correlation
`(x \|\| g)`, `(x \| gr(g, cor = FALSE))`	same	one term coded as `(x \| g)` codes it, with one SD per column and no correlation: `(1 + f \|\| g)` has `Intercept`, `fb` and `fc`, where lme4 gives `f` a column per level
`mvbind(y1, y2) ~ x`, `mvbf()`, `set_rescor()`	same	per-response families, ordinal ones included; `rescor` all gaussian or all student, with one shared `nu`; `mvbf() + family` fills only responses without a family (below)
`y \| trials(n)`, `weights(w)`, `cens(c)`, `trunc(lb=, ub=)`	same	`cens`/`trunc` need a CDF-carrying family; a DISCRETE bound is read as brms reads it (below)
`nl = TRUE`	same	a located prior places the start, a body the link cannot evaluate at 0 moves to the response mean, else provide `start` (`par_template()` names it); se.fit on the nonlinear mu not yet
`lf()`, `nlf()`	same	`nlf()` on any dpar, bodies chain to any depth (below)
`y ~ 0 + Intercept + x`, `bf(center = FALSE)`	same	the intercept is class `"b"` at zero, not class `"Intercept"` at the predictor means; refused in an ordinal family, as in brms
`custom_family()`	same idea	the lpdf is plain R over RTMB advectors, not Stan code
`cumulative()`, `multinomial()`, `categorical()`	same	the categories come from the response, and `refcat` names the reference category, as in brms
`dirichlet_multinomial()`, `brmsfamily("negbinomial2")`	same	`dirichlet_multinomial()` is `multinomial()` with a precision `phi`; `negbinomial2` is `negbinomial()` with `sigma = 1 / shape` and brms’s `gamma(0.4, 0.3)` default on `sigma`
`xbeta()`, `zero_inflated_beta_binomial()`, `hurdle_cumulative()`	same	brms’s density, links and `posterior_epred()`; `hurdle_cumulative()` takes brms’s threshold structures, `thres(gr = )` and `cs()` (see What changes)
`car(M, gr = g, type =)`	same, with `library(frmtmb.spatial)`	all four types; the term is in the extension package frmtmb.spatial
`(1 \| mm(g1, g2))`, `mmc(x1, x2)`	same	multi-membership, `weights =` and `scale =` included (below)
`me(x, sdx)`, `me(x, sdx, gr = g)`, `set_mecor()`	same	latent values integrated by Laplace; new data draws them from `N(x, sdx)`, as brms does (below)
`(x \| gr(g, by = f))`	same	one set of SDs and correlations per level of `f`, named `sd_g__Intercept:fa` as in brms; not with `cov =`
`ar(t, g)`, `ma()`, `arma()`, `cosy()`, `unstr()`	same	gaussian/student only; both `cov = FALSE` (brms’s default) and `cov = TRUE` (below)
`fcor(M)`	same	gaussian/student only; `M` in `data2`, `sigma` scales it, and a predicted `sigma` gives `D M D`, as in brms; `newdata` only on the fitted rows, as brms without `newdata2`
`sar(M, type = "lag")`, `sar(M, type = "error")`	same, with `library(frmtmb.spatial)`	gaussian/student only, constant `sigma`, as in brms; `rho` on `(1 / min eigenvalue, 1 / max eigenvalue)`, where brms takes the eigenvalue range itself
`data2 = list(W = W)`	same	also resolves compound expressions (below)
`bernoulli(link = "probit")`, and every other brms link name	same	the whole brms 2.23.0 roster (below)
`categorical()` on a numeric response	same	the distinct values are the categories, in numeric order: `11:13` is three categories with `11` the reference
`dirichlet()`, `logistic_normal()`, `brmsfamily("dirichlet2")`	same	the response is `cbind(y1, y2, y3)` or a matrix column, `refcat` as in brms; every row is held to 1.5e-8 of one, where brms holds the mean over rows; `lncor` takes class `"lncor"`
`mixture(..., order = "none")`	same	brms’s `order = "mu"`, which orders the components’ intercepts for a sampler, is refused, except on ordinal components, where it gives them one shared threshold vector, as in brms

Links

Every link name brms 2.23.0 accepts is accepted here under the same spelling, and computes the same function: logit, probit, probit_approx, cauchit, cloglog, softit, identity, log, log1p, logm1, inverse, 1/mu^2, sqrt, softplus, squareplus and tan_half. power12 is frmtmb’s own, for the tweedie power parameter.

Two cautions when porting a link:

stats::binomial(), stats::poisson() and stats::Gamma() validate their link string through stats::make.link() before frm() sees it, so a name stats has never heard of (probit_approx, softit, softplus, squareplus) is refused by the constructor rather than by frmtmb. Use frmtmb’s own family constructor instead: bernoulli(link = "softit"), Beta(link = "probit"), negbinomial(link = "softplus").
brms’s probit_approx is two different functions inside brms itself. Its Stan program uses Phi_approx(), the logistic of 0.07056 x^3 + 1.5976 x, while its R-side inv_link() answers pnorm(). frmtmb reproduces the Stan form, because that is the likelihood a brms fit was actually computed with.

The ordinal families take brms’s sets: cumulative() takes logit, probit, probit_approx, cloglog, cauchit and softit, and sratio() and cratio() take the same without softit, and acat() takes all six. Off the logit, acat() uses brms’s second density, as brms does.

What changes

brm(...) becomes frm(...); there are no chains or warmup, and a plain call is maximum likelihood rather than a posterior. REML = TRUE replaces priors as the small-sample correction for variance components.

A prior = c(prior(...), prior(...)) argument is kept as it stands, and so is prior = get_prior(...). frm() spells the argument prior, as brms does, and prior(), prior_() and prior_string() build the same specification set_prior() does; a prior object that brms itself built is translated row by row, so the argument works whichever package’s prior() was in scope.

A row applies whatever its prior string says. That is brms’s own rule, and it means a table edited in place, the ordinary gp <- get_prior(...); gp$prior[i] <- "normal(0, 20)" workflow, is honored. Earlier releases keyed the decision on the table’s source column, which brms does not update after such an edit, so the edited row was dropped.

The prior surface has brms’s names. default_prior() is the table of slots, and get_prior() is its alias, as in brms. validate_prior() checks a prior against the model and returns that table with the prior filled in; frm(prior = ) takes the result, and so does a default_prior() table edited in place. set_prior() takes vectors, set_prior("normal(0, 2)", class = c("b", "sd")), and a prior object prints as b ~ normal(0, 2). empty_prior() and as.brmsprior() exist too. The prior object itself is not a brms brmsprior: it holds parsed densities, and pr$class and the other brms columns are read from it.

What the fit MEANS changes. frm(prior = ) is MAP: the density is a PENALTY on the likelihood and the answer is one mode, not a posterior, so the reported log likelihood, AIC and anova() are penalized quantities. The estimates land close where the priors are weak: a kidney model with brms’s own priors gives sd(patient) 0.38 against brms’s posterior mean of 0.40. But “close” is a measurement, not a guarantee, and the two numbers are not the same quantity. Read a prior here as regularization.

Where a prior LANDS is brms’s placement, class by class. A density on sigma, shape, phi, nu or another distributional parameter is a density on that parameter itself, as in brms, rather than on its logarithm; a class = "Intercept" density is evaluated at the intercept at the MEAN of the predictors, which is the one brms constrains; and an ordinal family’s class = "Intercept" reaches the thresholds. On four of the shapes measured against Stan, frmtmb’s penalized objective is now brms’s posterior density up to one log(2) per lower-bounded parameter, a constant that moves no mode.

Which mode. The optimum is the mode in the INTERNAL coordinates, which is Stan’s target in unconstrained space, and not the mode of the parameter’s own density. A density on sigma enters the objective as log p(sigma) + log|dsigma/deta|, and what is maximized over eta = log sigma is that whole sum, log-Jacobian included.

Measured on a gaussian y ~ x fit with normal(2.5, 0.4) on sigma, over 150 rows drawn as set.seed(150), x <- rnorm(150), y <- rnorm(150, 1 + 0.5 * x, 2.2): frmtmb reports sigma 2.133344, the mode of L * p in the LOG-SIGMA coordinate is 2.133358, and the mode of L * p in the SIGMA coordinate is 2.126533. frmtmb matches the log-sigma mode, and it does so deliberately: that coordinate is Stan’s unconstrained target, and matching it is what keeps the two packages’ objectives one function, which is what everything above rests on. The gap to the sigma-coordinate mode is the same kind of remark as the log(2), and about the same size.

Two rows frmtmb accepts and brms refuses. A distributional class written without resp on a MULTIVARIATE model applies to EVERY response, following frmtmb’s own convention that a class-wide prior is class-wide; brms refuses such a row and asks for resp. Write resp to address one response.

frmtmb’s own set_prior() takes brms’s words with brms’s meaning, and that is a BEHAVIOR CHANGE for a script written before 0.53.0. A distributional parameter has two spellings, and brms accepts each only on the model the other does not:

the model	the spelling	the density is on
`bf(y ~ x)`	`set_prior(..., class = "sigma")`	sigma itself
`bf(y ~ x, sigma ~ 1)`, `sigma ~ x`	`set_prior(..., class = "Intercept", dpar = "sigma")`	the log-scale intercept of sigma’s predictor

frmtmb says the same and refuses each spelling by name on the other’s model, naming the one that applies. Before 0.53.0, set_prior("student_t(3, 0, 2.5)", class = "Intercept", dpar = "sigma") was accepted on a model with no sigma formula, where it meant a density on LOG sigma. To reproduce that placement, write sigma ~ 1 in bf(), which is the model that HAS a log-scale intercept; to get brms’s placement, drop the dpar = and write the class. The numbers move with the meaning: on a gaussian y ~ x with 200 rows and a maximum-likelihood sigma of 2.0103, normal(0, 0.5) gave 1.9965 under the old link spelling and gives 1.9429 under the class spelling, which is what brms::prior(normal(0, 0.5), class = "sigma") gave all along.

A prior class frmtmb keeps somewhere else (brms’s mixture theta, sdcar, car) is refused by name, with the class that reaches the same parameters where one exists, rather than translated into something else. So is a coef on a sd or cor row: frmtmb addresses a whole random-effect block there, and applying a narrowed row without its coef would put the density on every standard deviation of the block.

Because rows now apply, such a row in a get_prior() table stops the call where it used to be dropped in silence: on a model with s(), mo() or car(), blank that row or write the frmtmb class the message names. Every refused row of a table is named in one message, so one edit round is enough. Of 20 shapes a brms user would actually write, 3 stopped when this was measured: s(), gp() and mo(). A gp() table now carries over whole, since sdgp and lscale are classes here too. The rest translate, zero-inflated, hurdle and negbinomial models included.

A prior with a location DOES place a nonlinear start. brms uses its priors to place its sampler, and frm() reads them the same way for nonlinear parameters: a normal(), student_t() or cauchy() prior on a nonlinear coefficient starts it at the prior location, and a message names what was placed. The insurance-loss growth curve of the brms nonlinear vignette fits from its brms priors alone.

Without such a prior a nonlinear model still needs start, because frm() optimizes rather than samples and evaluates the objective AT the starting values, where a nonlinear body is rarely defined. par_template(formula, data) lists the parameter names before there is a fit; edit the result and hand it back as start =. brms has no counterpart to it, and needs none: brm(init =) takes Stan program names and brms initializes at random in a bounded range instead.
Posterior summaries become ML estimates: fixef() returns brms’s summary matrix with a standard error and a Wald interval in it, confint() gives Wald or profile intervals, anova() gives likelihood-ratio tests, and AIC()/BIC() replace loo().
posterior_predict() becomes simulate(), or predict(), which summarizes the same draws; posterior_epred() over newdata becomes fitted(newdata = ), or frm_linpred(newdata = , type = "response") for the bare vector. On a cens() response both draw the LATENT uncensored value, as brms does; simulate(censored = TRUE) applies the censoring mechanism instead.
The exact gp(x) is brms’s model, with sdgp * L * zgp for the field. The nugget is relative: 1e-12 * sdgp^2 on the kernel’s diagonal, where brms adds an absolute 1e-12. At sdgp near 1 the two are one number. The model does not depend on the units of y, and a fit depends on them only to the optimizer’s tolerance. brms’s fit does depend on them: at y * 1000 its kernel does not factor at the estimates. frmtmb factors the kernel over a farthest-first order of the distinct positions, where rounding moves the field least, and brms factors it in order of first appearance, so the b[i] columns of frm_sample() draws are not brms’s zgp[i] columns. Read the field through ranef() or the predictive methods.
sample_prior = "only" plus posterior_predict() is the same call on frmtmb.sample::frm_sample(), which samples the prior of the whole model as brms does. Without a sampler, frm_simulate(formula, data, prior = set_prior(...)) draws a parameter vector per simulation and returns it alongside the responses. sample_prior = "yes" and prior_draws() are in frmtmb.sample too.
mo() monotonic effects, mi() one-step imputation of continuous predictors, me(x, sdx) noise-free predictors and mi(sdx) measurement error, cs() (or brms’s other name for it, cse()) category-specific ordinal effects on every ordinal family, gp(x) / gp(x, k =) Gaussian processes, and mixture() families (including group-level latent classes via groups = ~g) all work. mo()/mi() support two-way :/* interactions with numeric terms; discrete predictors cannot be imputed in-model, the same restriction Stan has. Multiply imputed data can instead go through frm_multiple() (Rubin’s rules), the brm_multiple() analog. Being frequentist, that fit also needs pooled tests, which brms gets for free from the posterior: hypothesis() pools any function of the parameters, and anova() on the result compares two nested frm_multiple() fits with the D1, D2 and D3 rules of mice (D3 by default). fixef(), summary() and conditional_effects() of the result pool by Rubin’s rules where brms combines the draws of the imputations; plot(), as_draws_array() and nchains() read draws, so they are refused. Mixture fits are ML: expect multimodality, compare starts (frm_allfit()). pp_mixture() gives the component probabilities at the estimates in brms’s array.
cs() refuses a column written on both sides, where brms fits it. y ~ x + cs(x) is not identified: adding a constant to the population-level coefficient of x and subtracting the same constant from every one of its category-specific coefficients leaves the likelihood unchanged. brms builds both blocks and samples the ridge; maximum likelihood has no prior to hold it, so frmtmb refuses and says to write y ~ cs(x), which fits the same set of distributions. The same holds for y ~ f + cs(f) on a factor and for a constant cs() column. The check is a rank test, so it also covers a column that only a basis spells differently: y ~ poly(x, 2) + cs(x) is refused and y ~ I(x^2) + cs(x) is the edit that keeps the quadratic, while y ~ s(x) + cs(x) has no spelling that keeps both, because it is the smooth’s own unpenalized linear column that is aliased. mo() beside a cs() on the same predictor is refused too, a monotone function of the categories always lying inside the span cs() already gives them. brms builds both blocks in every one of these. mi() and me() are not checked, their design column being a placeholder at assembly, so me(x, sdx) + cs(x) is accepted as it is in brms.

Everything else about cs() is brms’s: a factor or a character predictor becomes treatment-contrast dummies with one coefficient per dummy per category boundary, named bcs_fb[k] the way brms names them, newdata is recoded against the fit’s levels, an unused level’s column is dropped as brms drops it, and a multivariate fit takes cs() in any response.
cs() on cumulative() and hurdle_cumulative() is fitted as brms fits it, with brms’s warning that the effects are experimental: each row reads its own thresholds, the fitted thresholds less the row’s offsets. Those can cross, and then the category between two crossed thresholds has a negative probability. brms’s density is undefined there and so is frmtmb’s. Where brms’s posterior_epred() returns the negative difference, fitted() returns NaN for every category of such a row, and simulate() returns NA.
thres() works on the five ordinal families. y | thres(K) ~ x sets the number of thresholds, and y | thres(gr = g) ~ x fits one threshold vector per level of g, named Intercept[g,k] as in brms. The section “Ordinal thresholds, thres()” of ?frm says where the count and the unobserved categories differ from brms. brms’s deprecated cat(K) is read as thres(K - 1), with brms’s warning.
hurdle_cumulative() takes thres(x = ), thres(gr = ) (each group keeps the hurdle in its density, as in brms), cs() (not together with thres(gr = ), as brms refuses it too) and the threshold structures "equidistant" and "sum_to_zero". The expected category that conditional_effects() shows by default is scored by the category codes 0..K, so the hurdle scores 0; brms scores each column by its position and reads one higher. With the logit link and a modeled disc, brms 2.23.0 cannot evaluate its own model: its Stan function reads a threshold past the last one for the top category, and Stan stops at “index out of range”. frmtmb fits that model.
rate() works as in brms on poisson(), negbinomial() and geometric(): the expected count is mu * d, and on the negative binomial the shape is multiplied by d too. On poisson() with the log link it is the offset(log(d)) model. The log-likelihood equals brms’s compiled program at the same parameters. fitted() reports mu * d, frm_linpred() the predictor without log(d).
subset(), index() and mi(x, idx = ) work as in brms: each response of a multivariate model is fitted on its own rows, with brms’s rule for missing values, and a predictor reads another response’s rows by index. The methods that return one row per observation need resp =, and on newdata they keep the rows where the subset is TRUE, both as in brms, and nobs() counts the data’s rows, or one response’s with resp =. Three differences, all because brms keeps every level of the whole data in a subsetted response and frmtmb keeps the levels the response’s own rows carry: a grouping factor has fewer levels (the marginal likelihood is the same, ranef() lists fewer, and |ID|-linked terms whose rows carry different levels are refused), a factor predictor level that occurs only outside the subset has no column where brms has an all-zero one, and an ordinal category that occurs only outside the subset has no threshold where brms has one. The section “Response subsets” of ?frm has the rest.
binomial(), beta_binomial(), zero_inflated_binomial() and zero_inflated_beta_binomial() without trials() are refused, as in brms, with brms’s message. A 0/1 response takes bernoulli() or y | trials(1); the stats::glm() convention of one implicit trial does not carry over.
ar(), ma() and arma() fit in both forms: brms’s default cov = FALSE and cov = TRUE. cosy() and unstr() have the covariance form only, as in brms. See the section “Within-group residual correlation” below.
glmer’s proportion-response idiom (weights = size) becomes y | trials(size) with either proportions or counts; weights() in the formula stays a frequency weight, as in brms.
se() works as in brms (meta-analysis), including se(x, sigma = TRUE) and the phylogenetic version with gr(g, cov = A). It reaches further than brms does: a family written with frmtmb_family() gets the term by declaring accepts_aterms = c(..., "se"), rather than by being gaussian or student. A family whose scale is not named sigma says which dpar the known standard error replaces with frmtmb_family(se_dpar =), or se_dpar = NA where it replaces none because the known standard error IS its whole scale. A dpar FORMULA does not pin the replaced scale: bf(y | se(s) ~ x, sigma ~ 1) is refused, because it estimates a parameter the density never reads. Pin it with a constant (sigma = 1) or write se(s, sigma = TRUE).
cens() and trunc() on a COUNT read a bound as brms 2.23.0 does. Right censoring at k is P(Y > k), brms’s poisson_lccdf(k | mu); left censoring is P(Y <= k); an interval is P(k < Y <= k2); and trunc(lb = k) is P(Y >= k), the one inclusive bound, through F(k - 1). trunc(lb = 0) is no truncation and is accepted. So a count recorded as “6 or more” is right censored at 5. Through frmtmb 0.72.0 a censoring bound on a count was inclusive instead, so a model written for that reading moves: on 200 poisson draws at lambda = 4 recorded as “6 or more” and right censored at 6, the log-likelihood moves by 30.2 (0.6 per censored row) and the estimate of lambda by 5.5 percent; right censored at 5, the fit is the old one exactly.
mm() ports with its weights = and scale = arguments and with mmc(); its other arguments have different spellings here (below).
A family added to a multivariate formula goes only to the responses that have no family yet. So bf(o ~ x) + cumulative() + bf(y ~ x) + gaussian() keeps o ordinal and makes y gaussian. brms gives the last family to every response, so the same line makes o gaussian there. Any number of bf() formulas can be added together, and lf() and nlf() take brms’s resp = to name the response they belong to.

Nonlinear formulas: `nl = TRUE` and `nlf()`

nl = TRUE says the RESPONSE formula is a nonlinear body, so only mu can be nonlinear. nlf() names the parameter instead, so any of them can be, and the two spellings meet in the same place:

# the same model, both ways
frm(bf(y ~ exp(b * x), b ~ 1, nl = TRUE), data = d)
frm(bf(y ~ a) + nlf(a ~ exp(b * x)) + lf(b ~ 1), data = d)

# and the model nl = TRUE cannot spell: a nonlinear sigma, linear mu
frm(bf(y ~ x) + nlf(sigma ~ a + b * z) + lf(a ~ 1, b ~ 1), data = d)

A parameter’s link is applied to its body’s value, as in brms: the sigma above is exp(a + b * z). Bodies chain to any depth (nlf(a ~ cc * x) + nlf(cc ~ exp(b))), are computed in dependency order, and a cycle is refused by name.

Two deliberate differences from brms, both supersets:

brms requires nl = TRUE on the bf() whenever the response formula names a nonlinear parameter. Here nlf() has already declared the name, so bf(y ~ a) + nlf(a ~ ...) is enough; the flag still works and means the same thing.
A body may name another distributional parameter of the same response, which reads that parameter’s per-row value. brms reads such a name as a data column and refuses the model when no column has it. A column still wins here, so a ported body keeps its meaning. This is the spelling for a variance function of the model’s own mean, which nlme writes varPower(form = ~ fitted(.)):

# sd = exp(ls) * |mu|^th, one joint maximum likelihood
frm(bf(y ~ x) + nlf(sigma ~ ls + th * log(abs(mu))) +
      lf(ls ~ 1, th ~ 1), data = d)

Distribution functions in a body

brms writes a nonlinear body in Stan, where Phi() and the rest belong to the language. Here the body is R evaluated onto an AD tape, and a function the body CALLS is looked up in RTMB first:

# a probit-shaped body: no RTMB:: prefix, no library(RTMB)
frm(bf(y ~ a * pnorm(b * x + c0), a ~ 1, b ~ 1, c0 ~ 1, nl = TRUE),
    data = d, start = list(beta = c(1, 1, 0)))

pnorm, qgamma, plogis, dbeta, besselK, matrix and their relatives all resolve this way, and RTMB wins over a function of the same name that you defined yourself. Write stats::pnorm() when you want the other one. The rule covers only the names RTMB replaces, and only where the body calls them, so your own helpers and the objects a body reads are untouched. The RTMB:: prefix keeps working and gives the identical fit.

Multi-membership

(1 | mm(g1, g2)) means what it means in brms: the observation belongs to several levels of one grouping factor, and its design row averages those levels’ effects. The pooled level set, the 1/J default weights and the scale = TRUE normalization of a supplied weight matrix are brms’s, checked against make_standata()’s J_*, W_* and Z_* arrays, so a ported model builds the same design matrix.

mmc(x1, x2) also ports: it is one random-slope coefficient of the term whose covariate value is member specific, member k taking argument k.

# brms and frmtmb, identical spelling
frm(bf(y ~ x + (1 | mm(school1, school2))) + gaussian(), data = d)
frm(bf(y ~ x + (1 | mm(school1, school2,
                       weights = cbind(share1, share2)))) + gaussian(),
    data = d)
frm(bf(y ~ (mmc(hours1, hours2) | mm(school1, school2))) + gaussian(),
    data = d)

brms’s remaining mm() arguments are spelled the way any other grouping term spells them here, and the parser names the replacement when it refuses one:

brms	frmtmb
`mm(g1, g2, cor = FALSE)`	same, or `(x \|\| mm(g1, g2))`
`mm(g1, g2, id = "q")`	same, or the `\|ID\|` key, `(x \| q \| mm(g1, g2))`
`mm(g1, g2, cov = A)`	same
`mm(g1, g2, by = cbind(f1, f2))`	same, one by-variable column per member
`mm(g1, g2, pw = )`, `dist =`, `by =` with `cov =`	no equivalent yet

?frmtmb-multimembership is the full page.

Noise-free predictors: `me()`

me(x, sdx) means what it means in brms: the column x is a noisy reading of a latent value with known measurement SD sdx, the latent values are normal with an estimated mean meanme and SD sdme, and the model uses the latent value. gr = g gives one latent value per level of g. Several me() terms have correlated latent values unless set_mecor(FALSE) is added, as in brms.

frmtmb integrates the latent values out with the Laplace approximation. For a gaussian response with linear me() terms the marginal likelihood is multivariate normal and the approximation is exact. The names follow brms: the coefficient of me(x, sx) is bsp_mexsx, and summary() reports meanme_mex, sdme_mex and corme__mex__mez, which brms reports only through variables().

On new data, a row has no latent value. predict() and simulate() draw one per replicate from N(x, sdx), as brms’s draws do. fitted() has no draws: it reports the mean of the expected response over that draw, by Gauss-Hermite, and the draw’s variance and quantiles in Est.Error and the Q columns. frm_linpred() stays at the observed x and adds the variance to se.fit. New data must hold sdx, as in brms. Under gr = g, rows at one level read one latent value. Two responses or two dpars read one latent value too, where brms draws one per predictor against its own Stan program.

?frmtmb-me is the full page.

Within-group residual correlation

ar(), ma(), arma(), cosy() and unstr() are written where brms writes them, in the model formula, and mean what brms means. Under cov = TRUE (and always for cosy() and unstr()) the residuals of one group become a single correlated draw, y_g ~ N(mu_g, D R D) with D the diagonal of that group’s sigma values. ?frmtmb-autocor has the full page; five things matter when porting.

The term grammar is brms’s. time is ONE variable name, or absent for the row’s position within its group, and gr is variable names crossed with :. An expression in either position is refused rather than evaluated, exactly as brms refuses it: ar(x + t, g) would index the correlation by the sum, and gr = g1/g2 would group the residual by the quotient. Compute the column first and name it.

cov = FALSE, brms’s default for ar(), ma() and arma(), ports unchanged. It is a different likelihood from cov = TRUE: mu gains a regression on the group’s earlier residuals, and every row keeps the family’s own density, so the likelihood conditions on each group’s first rows. It is brms’s likelihood to 1e-15 relative, and for one series it is the conditional sum of squares of stats::arima(method = "CSS"). As in brms, the coefficients are unconstrained, the lag is counted in rows, and weights(), cens(), trunc(), mi() and rescor = TRUE combine with it. fitted() and predict() give brms’s one-step mean, which reads the observed earlier residuals of newdata’s response. Where that response is NA, or absent, predict() fills it as brms does, with a draw at the row’s one-step mean, and fitted() fills it with its expected value, because a maximum likelihood fit has no draws to carry the fill through (brms’s fitted() fills with draws, and so does frmtmb.sample’s posterior_epred()). simulate() runs the recursion over its own draws. cosy() and unstr() have no cov argument in brms either and port unchanged.

# brms
brm(y ~ x + ar(week, subj), data = d)
brm(y ~ x + ar(week, subj, cov = TRUE), data = d)
# frmtmb
frm(bf(y ~ x + ar(week, subj)) + gaussian(), data = d)
frm(bf(y ~ x + ar(week, subj, cov = TRUE)) + gaussian(), data = d)

sigma is brms’s sigma. Under ar(), ma() and arma() with cov = TRUE it is the innovation SD, as in brms, whose cholesky_cor_ar1() divides by 1 - ar^2; under cosy() and unstr() it is the marginal SD, as in brms. So sigma(), summary(), a prior on sigma and the coefficients of a sigma ~ x model compare with brms’s directly. Pearson residuals divide by the marginal SD, sigma / sqrt(1 - phi^2) for an AR(1).

Other families get brms’s latent residual. gaussian() and student() are the two families with a residual covariance. On poisson, binomial and the rest, brms and frmtmb add a latent residual to the linear predictor, with its own SD, sderr. frmtmb integrates the latents by the Laplace approximation:

# brms
brm(cnt ~ x + ar(week, subj, cov = TRUE), data = d, family = poisson())
# frmtmb
frm(bf(cnt ~ x + ar(week, subj, cov = TRUE)), family = poisson(), data = d)

?frmtmb-autocor, section Latent residuals, has the model. brms’s refusals hold here too: cov = TRUE is of order one on these families, and cov = FALSE has no MA part.

The Laplace approximation of this model is not always good. On binary and ordinal responses, on a binomial with at most 3 trials, on counts with a mean below 1 and on a family whose dispersion trades off against the latent, frm() warns, because sderr and ar are biased: sderr reads 0.39 against an exact 1.10 on bernoulli and 0.26 against 0.65 on beta_binomial(). frm() refuses cov = FALSE on Gamma() or Beta(), where no number from the fit says whether it reads true, and a mean less than 2 marginal sds of the latent from the edge of an identity or inverse link. frm_sample() fits all of these and agrees with brms. ?frmtmb-autocor, section Where the Laplace approximation fails, has the tables.

Higher orders are allowed, and gaps are honored. brms limits cov = TRUE to order one; ar(p = 3) and arma(p = 2, q = 1) work here. And under cov = TRUE the lag between two rows is the distance between their positions in the global set of time levels, so a subject missing week 3 gets cor(week2, week4) = rho^2; brms indexes ar()/ma()/arma()/cosy() by position within the group and treats the missing row as no gap. On complete balanced groups the two agree. Under cov = FALSE the lag is counted in rows, as brms counts it.

The parameters appear in summary() under “Within-group residual correlation” and in confint_varcorr() under brms’s names (ar[1], cosy, cortime__1__2); autocor_matrix() returns the fitted correlation matrix under cov = TRUE. Under cov = TRUE, weights(), cens(), trunc(), se(), mi(), rescor = TRUE, mixtures and quadrature = TRUE are refused, because the likelihood no longer factorizes over rows. brms refuses the same core set. Random effects alongside the correlated residual are the point of the feature and are supported.

Variance functions (nlme `weights = varFunc`)

There are no varFunc objects. The sigma formula is the variance function: sigma uses a log link, so most of the nlme vocabulary is a direct translation, and the translation was verified against nlme::gls() under ML (log-likelihood agreement 1e-10 or better).

nlme	frmtmb
`varIdent(form = ~ 1 \| g)`	`sigma ~ 0 + g`
`varExp(form = ~ v)`	`sigma ~ v`
`varPower(form = ~ v)`	`sigma ~ log(abs(v))`
`varFixed(~ v)`	`sigma ~ offset(0.5 * log(v))`
`varComb(A, B)`	the two terms in one `sigma` formula

The coefficient of log(abs(v)) is nlme’s power theta. The translation is also more general than a varFunc: sigma takes the full predictor grammar, so random effects and smooths in the variance are ordinary spellings.

varPower(form = ~ fitted(.)), variance as a function of the model’s own mean, is spelled with nlf() (see the nonlinear-formulas section): nlf(sigma ~ ls + th * log(abs(mu))) + lf(ls ~ 1, th ~ 1) is sd = exp(ls) * |mu|^th, estimated by one joint maximum likelihood. Expect agreement with gls near 1e-2 rather than 1e-10: gls alternates the mean fit against iteratively updated fitted values and stops without a gradient check, while this is one joint maximization of the same likelihood. For a positive response the coupling often belongs in the family instead: Gamma() has sd proportional to the mean, and tweedie() estimates the power in Var = phi * mu^p directly.

Structural objects: `data2`

Matrices that are not columns of data go in data2, the brms spelling: the adjacency matrix of car() and the mesh triple of spde() (both frmtmb.spatial), and the matrices of gr(prec = ), gr(cov = ) and equalto().

frm(bf(y ~ x + car(W, gr = loc, type = "icar")) + gaussian(),
    data = d, data2 = list(W = W))
frm(bf(phen ~ 1 + (1 | gr(species, cov = A))) + gaussian(),
    data = d, data2 = list(A = A))

Two differences from brms. First, brms accepts a bare name in data2 and nothing else, while frmtmb evaluates the whole expression with data2 in front of the data mask, so gr(g, cov = solve(Q)) finds Q in data2. Second, data2 is optional: a matrix visible from the formula environment is still found, which is what frmtmb did before data2 existed.

Use data2 anyway. Its objects are stored on the fit, so the fit is self-contained: saveRDS() writes the matrices with it, and refit(), update(), drop1(), influence() and frm_multiple() re-assemble in a fresh session where the calling environment that built the matrix is gone. Without data2 those refits depend on the name still resolving where the model was written.

Method conventions

frmtmb takes brms’s argument names everywhere. Where lme4 or glmmTMB and brms disagree on a name, brms decides; the lme4 name is dropped and refused, not kept as an alias. What still follows stats is the SHAPE of what a standard generic returns:

brms-origin, same names and argument spellings: conditional_effects(), conditional_smooths(), make_conditions(), update_adterms(), hypothesis() (a delta-method Wald test here), pp_check() (simulate-based, through bayesplot’s ppc_* functions), prior_summary(), fixef(), ranef(), VarCorr(), ngrps(), variables() (brms’s parameter names, e.g. b_Days and sd_Subject__Days). VarCorr() and hypothesis() return brms’s structures, with the estimate and its standard error in the columns brms fills from draws; fixef(), ranef() and coef() take brms’s summary, robust and probs in brms’s positions and refuse summary = FALSE by name, because a fit has no draws to return. ranef() also takes brms’s pars and groups, and family() its resp. hypothesis() reads brms’s grammar too: every hypothesis states a relation, so write "x = 0" and not "x".
predict(), fitted() and residuals() return brms’s summary matrix: the columns Estimate, Est.Error, Q2.5 and Q97.5, one row per observation. A maximum-likelihood fit has no draws, so Est.Error is a standard error and the Q columns are a Wald interval, except on predict(), whose draws are SIMULATED and whose Q columns are the simulation’s own quantiles. An ordinal or categorical response gives brms’s n by 4 by K array from fitted() and brms’s P(Y = k) proportions from predict(). fitted(scale = "linear") on an ordinal fit with a cs() term gives one layer per threshold, as brms does, because the category-specific part moves the predictor differently at each threshold. residuals() takes brms’s newdata (by name: its second argument is type), which is the response there minus fitted(newdata = ), and answers every response of a multivariate fit in brms’s n by 4 by nresp array. confint() still takes level = rather than probs =.
A factor column of newdata may hold the levels as numbers, as brms reads it: Trt = 0 is the level "0". A number that is not a level is refused as a new level.

A matrix is not a vector, and every scalar summary of it changes. This is brms fidelity, not a defect, and it is the one change in this list that gives a WRONG NUMBER rather than an error. brms’s residuals() is a matrix, so mean(residuals(fit)) averages four columns and not one. Measured on y ~ x + f, n = 150:

expression	before	now
`mean(residuals(fit))`	-2.894331e-06	0.03897951
`sd(residuals(fit))`	0.9701508	0.8689341
`sum(residuals(fit)^2)`	140.2377	453.1844
`length(qqnorm(residuals(fit))$x)`	150	600

The residual vector is residuals(fit)[, "Estimate"], and the same holds for fitted() and predict(). Code written against brms already does this; code written against frmtmb before the move does not. insight::get_residuals(fit) returns the vector, as it does for a brmsfit.

predict() is the PREDICTIVE distribution, not the linear predictor. This is the one call whose meaning moved. brms’s predict() summarizes the posterior predictive distribution with observation noise in it, and frmtmb’s now does the same, by simulating from the fit: each replicate draws the parameters from their asymptotic normal law and then a response from the family’s own simulator. The LINEAR PREDICTOR, which predict() returned before, is frm_linpred(), with the same arguments it always had: frm_linpred(fit, type = "link") is the old predict(fit) and frm_linpred(fit, se.fit = TRUE) the old standard errors. A ported script that passes type or se.fit to predict() is refused by name and told where each one went, so the scale cannot change in silence.
Their ARGUMENT names are brms’s: predict(), fitted(), frm_linpred() and simulate() take re_formula and allow_new_levels, and lme4’s re.form and allow.new.levels are refused with the replacement named in the message. An argument none of them has is an error rather than a silently ignored name.
fixef() is brms’s summary matrix too, with brms’s row names (Intercept, sigma_Intercept), and vcov() covers exactly those rows. The per-predictor list fixef() used to return is fixef_by_dpar(), and fixef(flatten = TRUE) is one vector of estimates in confint()’s spelling. ngrps() is brms’s named list. summary() carries brms’s $fixed, $random, $spec_pars and $cor_pars slots and prints brms’s section headings, and brms’s EMPTY shapes with them: $random is NULL when the model has no grouping factors, and $spec_pars and $cor_pars are frames with no rows when there is nothing to report. On an ordinal fit the thresholds are fixef() rows named Intercept[1], Intercept[2], and vcov() covers them, as they are in brms.
coef() and ranef() name their coefficients as brms names them, Intercept and not (Intercept), so that anything pairing coef() with vcov() by name, lmtest::coeftest() among them, keeps every row. The design-column spelling is fixef_by_dpar().
conditional_effects() takes brms’s surface, too_far, select_points and spaghetti. spaghetti = TRUE needs draws: on a fit it takes band = "boot", where each refit is one curve, and it is refused with any other band. select_points measures the predictors only, where brms also drops the observations whose RESPONSE is far from its mean. A nonlinear predictor gets the Wald band by the delta method through the Jacobian of the body, where it needed band = "boot" before. conditional_smooths() draws each smooth term alone, with brms’s grid, keys and columns; on a fit the band is the delta method, on draws brms’s median and quantiles. Under re_formula = NULL a grid row reads a group’s fitted effects where it sets the group to a level the fit saw, and a new level where it sets none or an unseen one, as brms does with its own allow_new_levels = TRUE. band = "boot" and frm_bootstrap() hold every smooth at its fit and redraw only group-level terms. A row reads only its own f level’s block of gr(g, by = f), and each member of mm(g1, g2) is a grouping variable; the “New group levels” section of ?conditional_effects has the refusals. One invalid effect among valid ones is dropped with brms’s warning.
plot() of a conditional-effects object and of a hypothesis() result takes all of brms’s arguments (rug, stype, mean, the *_args lists, nvariables, ignore_prior, chars and others). frmtmb draws with base graphics, so plot = FALSE returns a list of plot objects that draw when printed, not ggplot objects, and a ggplot2 theme is not applied.
conditional_effects() on an ordinal fit draws one probability curve per response category, which is what brms draws under categorical = TRUE, and keys the effect "x:cats__" as brms does there. It is the default here because brms’s own default warns that it is treating an ordered factor as continuous; categorical = FALSE asks for that summary, the expected category number, with brms’s warning.
coef() is the per-group convention shared by all three packages: fixed effects plus conditional modes per grouping level.
predict(allow_new_levels = TRUE) draws an unseen level’s effect from its block’s estimated covariance, which is brms’s sample_new_levels = "gaussian", so the interval at an unseen level is wider than at a known one. fitted() takes the same sample_new_levels = "gaussian". predict() also takes "old_levels": each unseen level borrows one seen level’s effect, chosen once per call as brms chooses it, from every replicate’s draw of the group effects. brms’s "uncertainty" mixes the two per posterior draw, and a maximum-likelihood fit has none, so it is refused by name. fitted() takes "old_levels" too: the unseen level reads the chosen seen level’s fitted effect, and its standard error that level’s conditional variance.
update() takes brms’s newdata and records its name in attr(fit$data, "data_name"), which summary() prints. It also takes the data = of stats::update(), which brms refuses.
emmeans() takes brms’s dpar, nlpar, resp, epred and re_formula, with brms’s defaults. A multivariate fit without resp stacks its responses as brms’s rep.meas factor. Where brms summarizes draws, a nonlinear mu and epred = TRUE use the delta method through frm_lp_basis(); see ?frmtmb-emmeans.
plot(fit) draws residual diagnostics, not MCMC traces; for simulation-based residual checks use dharma_residuals(). brms’s N, variable and regex choose per-parameter panels of draws, so they are refused by name; plot(frm_sample(fit), variable = , regex = ) draws brms’s panels.
add_criterion() and summary(waic = TRUE) ask for posterior criteria, so on a fit they are refused by name; summary() prints AIC and BIC, and add_criterion(frm_sample(fit), "loo") stores the criterion as brms does.

Seeing the model code

brms users inspect the generated model with stancode(fit), which here stops with the reason and points to this section; so does standata(fit). The frmtmb analog is the objective function itself: an ordinary R closure that RTMB::MakeADFun() tapes, rebuilt deterministically from the fitted frame.

nll <- frmtmb:::build_objective(fit$frame)
nll(fit$estimates)   # joint negative log density at the estimates

Two properties matter when reading it:

The closure takes the named parameter list (the shape of fit$frame$par_template) and returns the JOINT negative log density: random-effect prior terms are inside it. The marginal (Laplace) likelihood the optimizer minimized is what MakeADFun(random = ) makes of it, exposed as fit$obj$fn over the flat outer-parameter vector; at the optimum fit$obj$fn(fit$opt$par) equals -logLik(fit).
All data is baked into the closure’s environment (the reason refits re-tape instead of recompiling), so environment(nll) holds the assembled frame, and the closure evaluates on plain numeric parameter lists without any tape.

RTMB::MakeADFun(nll, fit$frame$par_template, random = "b", map = fit$frame$map) reproduces fit$obj. There is no generated code to print: the closure IS the model, and body(), environment(), and the debugger work on it like on any R function.

The posterior half of the port

This page ports the model grammar and the frequentist workflow built around it. The sampling workflow ports as well, into the companion package frmtmb.sample. frm_sample() takes the same bf() object this page describes, assembles it through the same code path, and runs Stan’s NUTS on the frmtmb objective under brms’s default priors; the draws it returns carry the brms post-processing surface, with the real loo() on it. as_tmbstan() hands the objective to tmbstan directly.

Nothing in the grammar changes when you sample, which is why the two halves split cleanly. What does change is everything around the draws: the argument names of the call, the parameter spellings, and which brmsfit methods have no counterpart. vignette("brms-posterior", package = "frmtmb.sample") is that half, and vignette("posterior-diagnostics", package = "frmtmb.sample") covers the diagnostics of a sampled fit.

When you still want brms

With frmtmb.sample installed, a posterior, brms’s default priors and loo()-based comparison are all available here, and frmtmb.latent adds hidden Markov and latent class families, so the reasons to reach for brms itself are narrower than they look:

a marginal likelihood, and so a Bayes factor between models that are not nested. A point null INSIDE one model is now answered here: hypothesis() on draws made with sample_prior = "yes" reports the Savage-Dickey evidence ratio and posterior probability that brms reports, from the prior draws as brms does, for every parameter with a proper prior except the overall intercept, which brms leaves out too. What still needs brms is bridge sampling, which estimates the marginal likelihood itself;
reloo() or kfold() on a model whose LOO approximation fails, since both need refits that this package does not do;
a family or term neither package here implements, which frm_compat() answers per feature;
Stan code you intend to read, extend or hand to somebody else.

A practical workflow is model screening with frm(), which costs seconds, sampling the survivors with frm_sample(), and a final fit in brms for the questions in that list.