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, and 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. The development branch also includes complete-response Bernoulli spatial location intercepts or slope-only effects with a matching formula such as sd(site, level = "spatial") ~ z, for logit, probit, or complementary-log-log links. For a slope, the random-effect predictor must vary within sites; the SD-formula predictors must be constant within sites. These Bernoulli spatial SD routes have focused recovery and output checks, but their interval coverage and R-to-Julia parity are not established. 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"),
  ...,
  priors = NULL
)

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_family(), 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 named list of predictor models for the current missing-predictor routes. Most routes use one entry. The Gaussian response also supports two entries for two independent mi() predictors: either two Gaussian models or one Gaussian and one binary model, both with fixed effects. This two-predictor route is available with the native engine, not the Julia bridge. 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_family())) 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"). That Gaussian response also accepts two independent fixed-effect predictor models: two Gaussian models or one Gaussian and one binary model, on the native engine only. No other two-predictor family pair is supported. 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. Among currently supported models, missing = miss_control(response = "include") is implemented for univariate Gaussian response masks, bivariate Gaussian partial-response rows without dense known covariance, and the documented count and positive-continuous routes (including Poisson, NB2, Gamma and Binomial ordinary mu random effects with sd() formulas, and Poisson, NB2 and Gamma with one structured animal, relmat, phylo or spatial intercept that has its own sd(group, level = ...) formula, or with one animal intercept-plus-slope block that has one sd() formula per coefficient). missing = miss_control(predictor = "model") is implemented for one mi() missing predictor in the documented family-specific routes. A univariate Gaussian response also supports two independent mi() predictors with fixed-effect models: two Gaussian models or one Gaussian plus one binary model, on the native engine only. No other two-predictor family pair is supported. For one predictor in the 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).

engine

Computational engine. The default "tmb" uses the native drmTMB TMB backend. "julia" routes an admitted model to the DRModels.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(), beta_binomial(), skew_normal(), cumulative_logit(), biv_student(), and biv_lognormal() on the fixed-effect route only. Same-target live tests use coefficients and log-likelihood within 1e-4 and Wald SEs within 1e-3 relative against engine = "tmb". Measured gaps at DRModels.jl source pin 19e5e2fe for skew_normal(), cumulative_logit(), biv_student(), and biv_lognormal() sit inside those bars; see NEWS and tests/testthat/test-julia-family-skew_normal.R, test-julia-family-cumulative_logit.R, test-julia-family-biv_student.R, and test-julia-family-biv_lognormal.R. These comparisons do not establish the same tolerances for every admitted family. See vignette("julia-engine", package = "drmTMB") for the full family-by-route table and current boundaries. engine = "julia" needs a local DRModels.jl checkout and the optional JuliaCall package, and it is not a drop-in replacement for every native-TMB workflow. engine = "julia" is exported but experimental: the first call in an R session prints a one-time notice, and the bridge is not being widened beyond the families and routes listed above.

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 two-axis (q = 2) fixed-covariance 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. Only point-estimate recovery has been checked for both; bivariate relmat(..., Q = Q) REML is unsupported. 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. An internal experimental AGHQ plus Cox-Reid estimator is not a drmTMB() argument and is not what REML = TRUE runs. Public cumulative-logit random-slope fits remain maximum likelihood and have point-estimate recovery evidence only. engine = "julia" admits REML = TRUE only on named routes. Checked Gaussian examples include fixed-effect location-scale, an ordinary mu intercept with numeric-only mean predictors, mean-only and scale-side phylogenetic terms, and selected sd()/sd_phylo() location-scale-scale fits. Documented bivariate Gaussian shapes include the residual-only model with constant sigma1, sigma2, and rho12; matching q = 2 mean intercepts with phylogenetic, spatial, or supplied-relatedness structure; and a q = 4 phylogenetic route. The residual-only and q = 2 cells have same-target bridge tests; the q = 4 route has no same-target R-bridge parity claim. The bridge also admits the measured Poisson phylogenetic-intercept REML cell. It refuses a Poisson ordinary mu (1 | g) before REML dispatch. A separately tested direct DRM.jl fit for that shape does not make it available through the R bridge. Gaussian ordinary mean effects with factor, character, or logical predictors are also refused before Julia starts. These are route-specific point-fit or same-target checks, not blanket interval or coverage permissions. Requests outside the admitted shapes may be refused by the bridge or DRModels.jl; do not infer bridge REML support from a family name alone.

Under native REML = TRUE, admitted mean fixed-effect coefficients are marginalised in the TMB objective rather than estimated as ordinary fixed parameters. summary.drmTMB() and vcov.drmTMB() can still report finite Wald standard errors for those coefficients from the full sdreport covariance, but confint.drmTMB() does not currently form Wald or profile intervals for REML-integrated mean coefficients. Wald rows return NA bounds with conf.status = "wald_unavailable", and profiling such a mean coefficient can fail because the coefficient is not a top-level TMB parameter. Use confint(..., method = "bootstrap") when an interval for a REML mean coefficient is needed.

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) establish finiteness for Jeffreys-prior penalised binomial fixed-effect GLMs under their conditions, including full-rank design and links such as logit, probit, and cloglog. That result motivates the fixed-effect penalty; it does not prove finiteness or interval calibration for a mixed model fitted through TMB-Laplace. drmTMB has directly checked its three link implementations in named Bernoulli mixed-model settings. 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 scenarios measured this affected the large majority of converged fits. The comparisons for probit and cloglog 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: separated binomial GLM results do not calibrate an interval for this mixed-model estimator, and the package has not admitted one. 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.

priors

Optional fixed-effect normal priors, ridge penalties, or a Gamma prior on a random-intercept SD, created by drm_priors(). Fixed-effect priors support Gaussian identity and binomial logit models with complete data and unit weights; Gaussian models may also include one independent mean random intercept. The SD prior supports a complete-data, unit-weight Gaussian ML model with exactly one unlabelled independent mean random intercept, selected by group and (Intercept). It acts on the natural SD scale and can be combined with fixed-effect priors. The Julia bridge currently refuses factor predictors in Gaussian mean random-intercept formulas. drm_resolve_priors() previews selections. Active priors produce penalized/MAP point estimates. Covariance is local penalized-Hessian curvature, not calibrated sampling covariance. Confidence intervals and likelihood comparisons are not implemented. Cannot be combined with REML, estimator = "mspl", or penalty.

...

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_status: converged
#> convergence: 0
#> fixed effects:
#> mu: (Intercept) = 0.2281, x = 0.5928
#> sigma: (Intercept) = -0.6327