Fit a distributional regression model with TMB

drmTMB() is the main model-fitting entry point. The current implementation supports univariate Gaussian location-scale models, univariate Student-t and skew-normal location-scale-shape models, lognormal location-scale models, Gamma mean-CV models for positive responses, Tweedie mean-scale-power models for non-negative semicontinuous responses, beta mean-scale models for strict proportions, zero-one beta mean-scale-boundary models for continuous proportions with structural exact zeroes or ones, beta-binomial mean-overdispersion models for success counts, fixed-effect Bernoulli/binomial event-probability models, fixed-effect cumulative-logit ordinal location models, fixed-effect Poisson mean, zero-inflated Poisson, negative-binomial mean-dispersion, zero-inflated negative-binomial mean-dispersion, zero-truncated negative-binomial mean-dispersion, and hurdle negative-binomial mean-dispersion models for counts. Student-t, lognormal, Gamma, beta, ordinary Poisson, ordinary negative-binomial, beta-binomial, and zero-truncated negative-binomial mu formulas support ordinary unlabelled random intercepts and independent numeric slopes where documented. Every univariate mu formula may include a standard R offset() term, which enters the location linear predictor as a known constant. Read it on the family's own link scale: for log-link families (Poisson, ordinary negative-binomial, Gamma, Tweedie) offset(log(exposure)) is the usual exposure or effort rate model; for identity-link families (Gaussian, Student-t, lognormal, skew-normal) an offset is a known additive shift of the mean; and for logit-link families (Bernoulli/binomial, beta-binomial, beta, zero-one-beta, ordinal cumulative-logit) it is a known log-odds shift, which is a calibration term rather than an exposure. A zero-one-beta offset shifts only the interior beta component, leaving zoi and coi unchanged, and an ordinal offset shifts the latent location against fixed cutpoints. Zero-truncated and hurdle negative-binomial responses, every bivariate family, and Gaussian sufficient-statistic aggregation reject offsets: the first two renormalise their observed mean over a restricted support, so an exposure term would not scale the reported mean, and the others have no per-response offset contract yet. Offsets remain restricted to mu; other distributional parameters reject them. drmTMB also supports Gaussian random intercepts, independent numeric random slopes, and labelled or unlabelled correlated numeric random intercept-slope blocks in the location formula, known sampling covariance through meta_V(V = V) with deprecated meta_known_V(V = V) retained as a compatibility alias, residual-scale random intercepts and independent numeric random slopes in the scale formula, labelled mu/sigma random-intercept covariance blocks, and one or more group-level random-effect scale formulae such as sd(id) ~ x_group, plus phylogenetic random intercepts, one numeric phylogenetic random slope, and sd_phylo(species) ~ x_species direct-SD models in univariate Gaussian location formulas, Gaussian mu animal-model and user-supplied relatedness random intercepts and one numeric random slope, matching bivariate Gaussian mu1/mu2 location formulas, and matching labelled bivariate Gaussian mu1/mu2/sigma1/sigma2 phylogenetic location-scale blocks, coordinate-based spatial random intercepts and one numeric coordinate-spatial slope in univariate Gaussian mu, fixed-effect bivariate Gaussian distributional models, and matched labelled bivariate Gaussian mu1/mu2, sigma1/sigma2, and same-response mu/sigma random-intercept covariance blocks, including the first matching slope-only mu1/mu2 covariance block, the first all-four q=4 ordinary random-intercept covariance blocks, and predictor-dependent q=2 ordinary or phylogenetic corpair() regressions. Bivariate Gaussian location formulas may be written explicitly as mu1 = y1 ~ ..., mu2 = y2 ~ ..., or with mvbind(y1, y2) ~ ... shorthand when both responses share the same location predictors.

Usage

drmTMB(
  formula,
  family = stats::gaussian(),
  data,
  weights = NULL,
  control = list(),
  impute = NULL,
  missing = miss_control(),
  engine = c("tmb", "julia"),
  REML = FALSE,
  penalty = NULL,
  estimator = c("ml", "mspl"),
  ...
)

Arguments

formula

A drm_formula object created by drm_formula() or bf().

family

A response family, such as stats::gaussian(), student(), skew_normal(), lognormal(), stats::Gamma() with link = "log", tweedie(), beta(), zero_one_beta(), beta_binomial(), stats::binomial() with link = "logit", "probit", or "cloglog", cumulative_logit(), stats::poisson() with link = "log", nbinom2(), truncated_nbinom2(), or biv_gaussian(). Adding zi ~ predictors to a Poisson or nbinom2() model fits the corresponding zero-inflated count model. Adding hu ~ predictors to a truncated_nbinom2() model fits a hurdle count model whose nonzero counts use the zero-truncated NB2 component. The current bivariate Gaussian engine also accepts family = c(gaussian(), gaussian()) and family = list(gaussian(), gaussian()).

data

A data frame.

weights

Optional non-negative likelihood weights. These are row log-likelihood multipliers, not known sampling variances. For meta-analytic sampling variance or covariance, use meta_V() in the model formula instead.

control

Optional list passed to stats::nlminb(), or a drm_control() object when optimizer settings and fitted-object storage choices should be supplied together.

impute

Optional one-element named list of predictor models for the current missing-predictor routes. Bare formulas such as list(x = x ~ z) define Gaussian models for numeric missing predictors. Use impute_model() for explicit predictor families, such as list(treatment = impute_model(treatment ~ z, family = binomial())) for a binary predictor, list(score = impute_model(score ~ z, family = cumulative_logit())) for an ordered predictor, or list(habitat = impute_model(habitat ~ z, family = categorical())) for an unordered predictor, or list(cover = impute_model(cover ~ z, family = beta())) for a strict proportion predictor in (0, 1), or list(cover = impute_model(cover ~ z, family = zero_one_beta())) for a boundary proportion predictor in [0, 1], or list(cover = impute_model(success ~ z, family = beta_binomial(), trials = trials)) for a denominator-aware success/trial proportion predictor, or list(abundance = impute_model(abundance ~ z, family = poisson())) or list(abundance = impute_model(abundance ~ z, family = nbinom2())) for a count predictor, including list(abundance = impute_model(abundance ~ z, family = truncated_nbinom2())) for a positive zero-truncated count predictor, or list(biomass = impute_model(biomass ~ z, family = lognormal())) for a positive continuous predictor, or list(biomass = impute_model(biomass ~ z, family = Gamma(link = "log"))) for a Gamma positive continuous predictor, or list(biomass = impute_model(biomass ~ z, family = tweedie())) for a non-negative semi-continuous predictor with exact zeros. Grouped Gaussian covariate models use syntax such as list(x = x ~ z + (1 | group)); structured Gaussian covariate models use explicit syntax such as list(x = x ~ z + relmat(1 | line, Q = Q)). Most fitted routes use a univariate Gaussian formula containing one mi(x) location term and missing = miss_control(predictor = "model"). The first non-Gaussian response route also supports family = poisson() with one fixed-effect binary mi() predictor modelled by family = binomial().

missing

Missing-data policy created by miss_control(). The default keeps the existing complete-case behaviour. In the current fitted slices, missing = miss_control(response = "include") is implemented only for univariate Gaussian response masks and bivariate Gaussian partial-response rows without dense known covariance. missing = miss_control(predictor = "model") is implemented for one mi() missing predictor in a univariate Gaussian location model: numeric Gaussian predictors may use a fixed-effect, one random-intercept, or one intercept-only structured Gaussian impute formula; binary, ordered categorical, unordered categorical, strict beta/proportion, zero-one beta boundary-proportion, beta-binomial denominator-aware proportion, Poisson, negative-binomial, or zero-truncated negative-binomial count, lognormal positive continuous, Gamma positive continuous, and Tweedie semi-continuous predictors may use one fixed-effect family-aware impute_model(). Poisson, binomial, nbinom2, beta, Gamma, lognormal, beta-binomial, and student responses each support one fixed-effect binary missing predictor with a Bernoulli/logit impute_model(), complete responses, and no random or structured response terms. nbinom2() also admits one fixed-effect Gaussian impute_model() for a continuous missing predictor (k = 1; no grouped/structured predictor model, no k = 2).

engine

Computational engine. The default "tmb" uses the native drmTMB TMB backend. "julia" routes an admitted model to the DRM.jl bridge: Gaussian, bivariate Gaussian, Student-t, lognormal, Poisson, NB2, Gamma, Beta, and Binomial on the fixed-effect route (several with an additional large-p phylo() or coupled location-scale phylo route), plus truncated_nbinom2(), zero_one_beta(), tweedie(), and beta_binomial() on the fixed-effect route only. Each admitted family carries a same-target parity receipt against engine = "tmb" (coefficients and log-likelihood within 1e-4, Wald SEs within 1e-3 relative); see vignette("julia-engine", package = "drmTMB") for the full family-by-route table, which capability-ledger row documents each receipt, and current boundaries. engine = "julia" needs a local DRM.jl checkout and the optional JuliaCall package, and it is not a drop-in replacement for every native-TMB workflow.

REML

Logical; use restricted maximum likelihood where the selected engine supports it. Native engine = "tmb" restricts the likelihood by marginalising the admitted fixed-effect mean coefficients. Validated Gaussian routes include ordinary random effects; mean-side phylo(), spatial(), animal(), and relmat() effects; known sampling covariance through meta_V(); non-unit likelihood weights; and selected scale-side random or structured effects. Mean-side spatial(), animal(), and relmat() REML is deliberately limited to an unlabelled intercept or independent intercept-plus-one-numeric-slope term with a constant residual scale (sigma ~ 1, with no sigma random effect). Slope-only, labelled, multiple-slope, and matched non-phylogenetic mean-scale REML remain outside that route. The two bivariate non-phylogenetic exceptions are matching labelled fixed-covariance q2 location intercepts for spatial(1 | p | site, coords = coords) and relmat(1 | p | id, K = K) in mu1 and mu2, with intercept-only sigma1, sigma2, and rho12, complete response pairs, unit weights, no known meta_V() covariance, and no additional ordinary random effect, direct-SD formula, or corpair() regression. Both have point-fit-recovery evidence only; bivariate relmat(..., Q = Q) REML remains deferred. Aggregation and ordinary direct sd() scale formulae also remain unsupported under REML. For binomial() models, the bounded native route requires exactly one ordinary unlabelled mu random intercept or exactly one independent numeric slope. Fixed-only and multiple-term binomial models, correlated or labelled covariance blocks, structured effects, missing-data engines, and other extensions are not admitted under REML; use REML = FALSE for those models. The binomial route has diagnostic parity and finite-uncertainty evidence, not calibrated interval or coverage evidence. Every other non-Gaussian family rejects REML = TRUE; use REML = FALSE. The package-private AGHQ plus Cox-Reid (O3) estimator is not a drmTMB() argument and is not what REML = TRUE runs. Public cumulative-logit random-slope fits remain maximum likelihood (point_fit_recovery). engine = "julia" fits REML = TRUE on a route-by-route basis, not uniformly: the fixed-effect Gaussian location-scale model, the bivariate q = 4 phylogenetic route, the location-scale-scale sd()/sd_phylo() routes, an ordinary Gaussian mu random intercept, and large-p Poisson phylo() (a Cox-Reid Laplace route native engine = "tmb" does not have) currently fit REML = TRUE; every other bridge route refuses it, some with a raw DRM.jl error rather than a drmTMB-worded one. See docs/design/261-reml-by-route.md for the full route-by-route table before relying on bridge REML for a model not listed there.

penalty

Optional penalty / prior built by drm_phylo_penalty(), or NULL (default) for plain maximum likelihood. A non-NULL penalty switches the fit to a penalized / maximum-a-posteriori (MAP) estimator that regularises a weakly-identified phylogenetic standard deviation; the fit is labeled MAP and logLik() returns the unpenalized data log-likelihood. Native engine = "tmb" only.

estimator

Estimator for the native TMB route. The default "ml" preserves ordinary maximum likelihood. Experimental "mspl" implements the clean-room maximum softly-penalized likelihood criterion for one complete Bernoulli or grouped-binomial model with one ordinary q = 1 or correlated q = 2 grouping block. Supported links are "logit", "probit" and "cloglog". Kosmidis and Firth (2021, Theorem 1 and Section 3.1) prove the Jeffreys penalty gives finite estimates for any link whose working weight vanishes in both tails, which these three do; drmTMB's own TMB-Laplace evidence for them is recorded in docs/dev-log/simulation-artifacts/2026-08-11-mspl-nonlogit-links/. Three limits are worth stating. The soft-penalty scale below is a logit delta-method constant, kept unchanged for all three links because using it for probit and cloglog was measured to cost about 1% of one standard error – measured at q = 1 with two fixed-effect columns, and not at q = 2. Under deep separation with a random slope the standard error is frequently unavailable for every link, logit included (see the NA note below); in the most extreme cells measured this affected the large majority of converged fits. And the evidence behind the two new links used Bernoulli responses with two fixed-effect columns, so grouped-binomial and wider designs inherit the n_eff extrapolation described next without direct measurement, for any link. For grouped binomial data, the soft-penalty scale c_n = 2 * sqrt(p / n_eff) uses n_eff = sum(trials * frequency), the retained total number of Bernoulli trials, in place of the source paper's row count n; the two agree exactly when every row is a single Bernoulli trial (trials = frequency = 1). This is a deliberate grouped-data extension, not a claim that Sterzinger and Kosmidis's (2023) results carry over unchanged. vcov() and the std_error column of summary() are available: they invert the Hessian of the unpenalized Laplace log-likelihood evaluated at the MSPL estimate, so the penalty is used to obtain a finite estimate but not to describe sampling variability. A standard error is reported; an interval is not claimed. Likelihood comparisons (logLik(), AIC, BIC, anova()), profiles, and confint() remain deliberately unavailable. Do not form coef ± 1.96 * se by hand: Kosmidis and Firth show that Wald intervals in this setting fail to cover regardless of the nominal level, a failure that persists even for profile penalized-likelihood intervals. When the information matrix is not positive definite the standard errors are NA with a drmTMB_mspl_wald_unavailable warning rather than a fabricated number.

...

Reserved for future model options.

Value

A drmTMB fit object.

Examples

set.seed(20260525)
dat <- data.frame(
  y = 0.2 + 0.6 * seq(-1, 1, length.out = 24) + rnorm(24, sd = 0.5),
  x = seq(-1, 1, length.out = 24)
)
fit <- drmTMB(bf(y ~ x, sigma ~ 1), data = dat)
fit
#> <drmTMB Gaussian location-scale fit>
#> observations: 24
#> estimator: ML
#> logLik: -18.87
#> convergence: 0