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. Binomial mu formulas may include standard R offset() terms on the logit-event-probability scale. Poisson, ordinary negative-binomial, and zero-truncated negative-binomial mu formulas may include standard R offset(log(exposure)) terms for exposure or effort, 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,
  ...
)

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", 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(). The first Poisson-response route supports one fixed-effect binary missing predictor with a Bernoulli/logit impute_model(), complete count responses, and no zero-inflation, random, or structured response terms.

engine

Computational engine. The default "tmb" uses the native drmTMB TMB backend. The "julia" compatibility bridge is halted and deferred for future work; it is retained only so existing objects and code can be inspected, not as a current fitting route.

REML

Logical; use restricted maximum likelihood where the selected engine supports it. Native engine = "tmb" keeps REML Gaussian-only and restricts the likelihood by marginalising the admitted fixed-effect mean coefficients. Validated univariate 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. The halted engine = "julia" compatibility bridge is not a supported estimator or REML route. Use native engine = "tmb" for fitting and use REML = FALSE for likelihood-ratio tests, AIC/BIC comparisons across different fixed-effect formulas, non-Gaussian models, and currently unsupported extensions.

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.

...

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