Formula-Based Regression Models via RTMB • frmtmb

frmtmb fits regression models that you specify with a brms-style formula grammar. Estimation is maximum likelihood, with the Laplace approximation for latent effects. The backend is RTMB. The package generates each model’s objective as an R closure and differentiates it on the TMB AD tape. There is no MCMC, no Stan, and no compilation at run time. A refit re-tapes in milliseconds. The bootstrap, influence, and multi-start machinery build on that speed.

Documentation: https://aforren1.github.io/frmtmb/

Statement of need

brms gave R a formula grammar that covers a very wide model space. It reaches that space only through Bayesian estimation in Stan. A user who wants the same grammar under maximum likelihood must divide one model across several packages, accept a smaller model, or wait for a sampler. frmtmb closes that gap. It fits the brms grammar by maximum likelihood with the Laplace approximation, and it compiles nothing at run time.

The gap is real because the existing frequentist packages each stop at a different place.

lme4 fits mixed models in a small set of GLM families. Only the mean gets a formula, and each grouping factor gets one unstructured covariance.
glmmTMB adds many families, structured random-effect covariances, and separate formulas for dispersion and for zero inflation. The predictor grammar is not brms’s: there are no nonlinear formulas, no multivariate responses with residual correlation, no monotonic or measurement-error terms, and no correlation of effects across formulas.
gamlss, and its successor gamlss2, give every distributional parameter its own additive predictor, which is the part lme4 and glmmTMB lack. Random effects are additive terms there rather than an lme4 grammar with structured covariances, and the spelling is not brms’s.
drmTMB comes closest. It fits distributional regression by maximum likelihood on a compiled TMB template. Its formulas resemble brms: bf(), one formula for each parameter, and |p| labels. It fits one or two responses, and it refuses a pair of different families. For the beta and NB2 families, it refuses effects correlated across formulas. It also refuses random effects in some parameters, for example the beta precision and the Student-t scale. It has no nonlinear formulas, smooths, or custom families. Its REML covers the gaussian and binomial families only. These statements describe drmTMB 0.7.0 on CRAN.

frmtmb supplies the combination: distributional regression on every parameter of the family, nonlinear formulas, multivariate responses with rescor, effects correlated across formulas with |ID|, custom families written as plain R log-densities, and the lme4 random-effect grammar with structured and spatial covariances, all under maximum likelihood or REML. The package generates each model’s objective as an R closure and differentiates it on the RTMB tape, so it does not select from a fixed set of likelihoods. Features that would each need a bespoke implementation become ordinary code paths, and they compose with each other.

The second need is migration. brms code ports by changing brm() to frm(). The priors can stay where they are: frm(prior = ) takes brms’s own spelling, prior() builds the specification set_prior() does, and a prior object brms itself built is translated. The fit is then MAP, so a prior is a penalty on the likelihood rather than a posterior. A measured audit of the brms vignettes puts about 7 of 10 of their model calls through that transform unchanged; over every call the vignettes make, post-processing included, about 4 of 10 run unchanged, and most of the rest need one spelling change or get a refusal that names the replacement. A user can therefore screen models at millisecond speed, then return to brms for the final Bayesian fit with the same formula. vignette("brms-migration") maps the features and states, under “When you still want brms”, the cases that belong in brms.

A scoped core, with the rest in companion packages. The core package fits models and reports on them. Everything that needs another engine or another literature ships beside it: NUTS sampling, ordinary differential equations, discrete latent states, and drift-diffusion response times each live in their own package in this repository. That keeps the core’s dependencies on CRAN, keeps its promise of no compilation exact, and makes the part that matters reviewable on its own. The seam is public: frmtmb_family(), frmtmb_structure() and the registration functions are the same interface the companion packages use, so a family written outside this repository reaches the core the same way hmm() does.

Audience. Applied statisticians and quantitative researchers who already write brms formulas, and who need a likelihood-based answer: for a maximum-likelihood workflow, for a model that must fit in seconds inside a simulation or a bootstrap, for a report that asks for confidence intervals and likelihood-ratio tests, or for a teaching setting where no toolchain can be installed. Fields where this comes up include psychology and psychophysics, ecology, pharmacometrics, and meta-analysis.

Related work: glmmTMB and drmTMB are the closest relatives. glmmTMB is a mature and fast TMB-based mixed-model package. frmtmb matches its fits where the models overlap and follows its conventions in several places. drmTMB is the closest in scope. In 27 comparisons of models that both packages fit, from beta and NB2 mixed models to animal-model and bivariate fits, the two log-likelihoods agree to within 3e-8. Each comparison states the map between the two parameterizations (tests/testthat/test-drmtmb-agreement.R). drmTMB also does work that frmtmb does not do. It gives a random-effect standard deviation its own formula, sd(g) ~ w, in gaussian models, where w must be constant within the group; frmtmb reaches the same model only through a nonlinear formula. It takes a phylogeny, a pedigree, or a pair of trees directly, where frmtmb needs the covariance matrix. It accepts a full known sampling covariance for meta-analysis. It fits the zero-truncated NB2 family. It imputes a binary missing predictor, and it draws worm plots and centile charts. Its documentation grades each model route by the evidence behind its estimate and its interval. BayesRTMB is a Bayesian-first, Stan-like modeling layer on the same RTMB backend. qbrms also reads brms syntax and sends it to a different engine, INLA, which fits approximate posteriors for latent Gaussian models; frmtmb instead maximizes the likelihood and takes families and nonlinear bodies that INLA’s model class does not cover. brms itself is the right tool when you want priors and full posterior inference; the shared grammar is meant to make movement between the two easy.

Installation

frmtmb is not on CRAN yet. Install the development version from GitHub:

# install.packages("remotes")
remotes::install_github("aforren1/frmtmb")

Eight companion packages live in this repository, each installed the same way and each documented on its own part of the site.

package	adds
frmtmb.sample	NUTS sampling through tmbstan, and the posterior method surface
frmtmb.latent	discrete latent states: `hmm()`, `lca()` and `mixture_mvn()`
frmtmb.ode	ordinary differential equation dynamics: `frm_ode()`
frmtmb.eam	evidence accumulation models: `wiener()`, `gddm()`, `lba()`, `rdm()` and `wiener_gng()`
frmtmb.survival	survival families: brms’s `cox()` with an M-spline baseline, and the flexible parametric `royston_parmar()`
frmtmb.spatial	spatial random fields: brms’s `car()` and the SPDE Matern field `spde()`
frmtmb.learn	reinforcement-learning families: `bandit2arm_delta()`, `bandit2arm_dual()`, `prl_fictitious()`, `bandit4arm2_kalman_filter()`, `ts_par7()`, `igt_pvl_delta()`, `igt_orl()` and `rlddm()`, each parameter an ordinary distributional parameter with its own linear predictor
frmtmb.coupling	coherence and phase between two signals: the complex Wishart family `cross_wishart()`, whose coherence and phase are distributional parameters with their own linear predictors

remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.sample")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.latent")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.ode")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.eam")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.survival")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.spatial")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.learn")
remotes::install_github("aforren1/frmtmb", subdir = "extensions/frmtmb.coupling")

frmtmb has no package that installs everything. To install every extension at once:

# 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))
}

Example

library(frmtmb)
data(sleepstudy, package = "lme4")

fit <- frm(bf(Reaction ~ Days + (Days | Subject)), family = gaussian(),
           data = sleepstudy)
summary(fit)

# distributional regression: model the residual SD too
fit2 <- frm(bf(Reaction ~ Days + (Days | Subject),
               sigma ~ Days + (1 | Subject)), family = gaussian(),
            data = sleepstudy)
anova(fit, fit2)

hypothesis(fit, "sd_Subject__Days^2 / (sd_Subject__Days^2 + sigma^2) > 0",
           class = NULL, method = "boot")

vignette("habit") is a longer worked example: it replicates a published response-preparation model, fits every participant in a few seconds, and then fits the hierarchical version that the original per-participant procedure could not express.

Status

Pre-release. The goal is a CRAN release. Validation has three layers:

Every model class is compared with an exact external reference: glmmTMB, drmTMB, lme4, mgcv, nlme, MASS, survival, nnet, GLMMadaptive, quantreg, mice, closed-form marginals, or hand-written maximum likelihood. The core suite and each companion package’s suite run on every change, and each companion package has its own check workflow.
The model-building layer is compared with brms itself. Design matrices, random-effect structures, and special-term data agree with brms::make_standata() to near machine precision. An opt-in tier verifies that our log-likelihood equals the Stan program’s log density at the estimate.
A pairwise grammar fuzzer sweeps feature combinations against metamorphic invariants. The resulting compatibility map is queryable with frm_compat(), which the companion packages contribute their own rows to, and is published as the feature compatibility article.

Model grammar

bf() uses the brms spelling, and the family is a separate argument: frm(bf(y ~ x), family = poisson(), data = d). With no family, frm() fits gaussian, as brms and glmmTMB do. See “Statement of need” above for the port from brm().
Every distributional parameter can take its own formula with the full predictor grammar (bf(y ~ s(x) + (1 | g), sigma ~ s(z) + (1 | g))), or a constant.
Random effects use the lme4 syntax and add structured covariances: us, diag, homdiag, cs, homcs, ar1, hetar1, toep, homtoep, ou, spatial exp/gau/mat over num_factor(x, y) coordinates, reduced-rank rr(d =), and known structure gr(g, cov = A) / gr(g, prec = Q) / equalto(). The |ID| syntax correlates effects across formulas, including known-matrix blocks (a multi-trait animal model is one Kronecker block). Multi-membership terms use mm(g1, g2) with weights and member-specific mmc() covariates.
Within-group residual correlation uses brms’s R-side terms: ar(), ma(), arma(), cosy(), and unstr(), for gaussian and student responses, validated against nlme::gls under ML and REML. brms’s default cov = FALSE form of ar(), ma() and arma() is brms’s own likelihood, a regression on earlier residuals. brms’s fcor(M), a fixed residual covariance over the observations, is in frmtmb too, and its sar() in frmtmb.spatial.
brms’s car(M, gr = g, type =) spatial Gaussian Markov random field, with all four of its types (escar, esicar, icar, and bym2 with brms’s scaling convention), and spde(fem, gr = node), a Matern field over a finite-element mesh, are in frmtmb.spatial. They register through core’s covariance-structure seam, frmtmb_register_covstruct().
mgcv smooths s()/t2() work in any linear predictor. Matrix-covariate terms give scalar-on-function, function-on-scalar, and function-on-function regression.
gp(x) fits an exact Gaussian process and predicts with kriging. gp(x, k = 30) uses the Hilbert-space approximation with brms’s exact input convention, so the same call is the same approximation. Up to three dimensions are supported, with per-dimension lengthscales; iso = TRUE shares one.
mo() fits monotonic effects. mi() imputes continuous predictors in one step, and mi(sdx) models measurement error. cs() fits category-specific ordinal effects. mo() and mi() support two-way interactions.
Multivariate models use mvbf, per-response families, and rescor. Nonlinear formulas use nl = TRUE. Inside a nonlinear body, a bare pnorm() or qgamma() is the tape-capable version, so a process model reads as it would on paper. multinomial() and dirichlet_multinomial() take a matrix response.
The families include the usual GLM(M) set plus student, tweedie, compois, beta-binomial, skew-normal, ex-gaussian, weibull, shifted lognormal, and quantile regression (asym_laplace), with zero-inflated and hurdle variants and four ordinal families. mixture() fits finite mixtures, including group-level latent classes (mixture(..., groups = ~g)) with class-specific random effects for growth-mixture models. categorical() fits nominal responses, dirichlet() and logistic_normal() compositions (rows of proportions that sum to one), and von_mises() circular responses. brms’s cox() proportional hazards model is in frmtmb.survival. custom_family() takes a plain R log-density; section 11 of vignette("case-studies") writes one with a data-bounded link and checks it against the built-in family it shadows.
The addition terms are weights(), trials() (counts or proportions), cens(), trunc(), se() (meta-analysis), mi(), and vint()/vreal(). A companion package can register one of its own.
The companion packages extend the same grammar. frmtmb.latent adds hmm(K) for hidden Markov models, with covariate-dependent transitions and forward-backward decoding, and lca(K) for latent class analysis with class-membership regression, and mixture_mvn(K, D, model =), which fits multivariate gaussian components over mclust’s covariance taxonomy, with covariate-dependent means and gating. frmtmb.ode adds frm_ode(), which solves compartment models inside nonlinear formulas (population pharmacokinetics), with repeated dosing, infusions, and estimated bioavailability. frmtmb.eam adds wiener(), the drift-diffusion first-passage density, with a formula for each of the drift rate, boundary separation, non-decision time and starting bias, and Ratcliff’s across-trial variability as three more; gddm(), the generalized drift-diffusion model with time-varying drift and collapsing bounds, solved on a fixed grid; and lba(), the linear ballistic accumulator, which races any number of accumulators and so reaches choices with more than two alternatives.

Estimation and inference

Estimation is ML or REML. quadrature = TRUE uses adaptive quadrature for scalar random effects. set_prior() and prior() give MAP estimation with the brms spelling, nlpar = and resp = included, and the same call carries hard bounds through lb/ub: the prior vocabulary is the only way to bound a parameter, and it reaches every one that brms reaches, including the residual-correlation classes ar, ma, cosy, cortime and rescor. frmtmb_control() offers pluggable optimizers, profile = TRUE for many-coefficient models, sparse_x = TRUE for many-level fixed factors, autoscale = TRUE for badly scaled predictors, and verbose = for timed progress on slow fits. par_template() reports the names and starting values a model takes, before there is a fit to ask.
confint() offers Wald, profile, likelihood-root, and bootstrap intervals. hypothesis() tests expressions of coefficients, with natural-scale sd_/cor_ names for ICC-type quantities. anova() and drop1() give likelihood-ratio tests. frm_bootstrap(), influence() with cooks.distance() and dfbetas(), frm_allfit(), and frm_multiple() (Rubin pooling, including variance components and hypotheses) complete the set.
Simulation goes both ways. simulate() draws from a fit. frm_simulate() draws from a bare design, with natural-scale parameters (sd_g__Intercept = 0.7) or with parameters drawn per replicate from set_prior() specifications. That is the sample_prior = "only" prior-predictive workflow, without MCMC.
Posterior inference lives in frmtmb.sample. Its frm_sample() runs NUTS on the fitted objective, or on a formula with brms’s default priors, and returns a full draws surface (posterior_epred(), posterior_predict(), hypothesis(), pp_check(), loo()). Its check_laplace() measures the Laplace and Wald approximations against the sampler on the same objective.
Diagnostics include one-step-ahead residuals calibrated for censored, truncated, and ordinal responses, deviance residuals across the GLM families, and response-scale se.fit for every family through the joint delta method, which is what fills brms’s Est.Error column on fitted() and residuals(). DHARMa, pp_check(), plot(), and conditional_effects() are supported. Hooks for emmeans, marginaleffects, and insight/easystats are registered.
Curve inference reads any fitted curve on a grid: frm_curve() gives pointwise and simultaneous bands, frm_curve_deriv() its derivatives, and frm_curve_feature() the location of a peak or a crossing, each with a delta-method standard error.

Life cycle

frmtmb is maturing. The package is not yet on CRAN.

The model grammar is stable. It follows brms, so it changes only when brms changes. Formulas that fit today will fit in later versions.
The fitted-object API is stable, with one exception in flight. frm(), the accessor methods (coef(), confint(), vcov(), predict(), and the rest) and the family constructors keep their current behavior. The exception is the move to brms’s RETURN SHAPES: predict() is brms’s predictive summary and the linear predictor it used to return is frm_linpred(); fitted(), residuals(), fixef(), vcov(), ngrps() and summary() return brms’s shapes. NEWS.md says what breaks and what replaces each one.
The extension interface is public and versioned. frmtmb_family(), frmtmb_structure() and the registration functions are documented as one contract, and the companion packages in this repository are written against it with no reach into internals. frmtmb.eam was built that way as a test of it. Each companion package pins the core version it needs, and the pin moves when the interface does.
Internal structure can change. Fields of the fitted object that no exported method reaches are not part of the API. Use the accessors.
Some parts are still moving. Multivariate coverage of the post-fit methods and the mixture families gain features between releases.

Version numbers stay below 1.0 until the CRAN release. Breaking changes are listed in NEWS.md.

The package is actively developed and maintained. Planned work: full multivariate support in every post-fit method, more reference comparisons, and CRAN submission.

Citation

Run citation("frmtmb") in R. Please also cite RTMB and TMB, and tmbstan if you sampled.

Contributing

See CONTRIBUTING.md. Everyone taking part must follow the Code of Conduct.

See SPEC.md for the design and NEWS.md for the changelog.