Function map and cheat sheet • drmTMB

Use this page when you know the scientific question but need to find the smallest useful drmTMB workflow. It is a map, not a capability ledger: use What can I fit today? before relying on an unfamiliar random-effect or inference route.

Quick reference

Goal	Start with	Read next
Specify and fit a first model	`bf()`, a family object, `drmTMB()`	Getting started
Decide whether the fit is interpretable	`check_drm()`, `summary()`	Model workflow
Extract a parameter or make a prediction table	`coef()`, `fixef()`, `predict_parameters()`	Distributional outputs and adequacy
Work with repeated, related, or spatially structured data	`ranef()`, `phylo()`, `spatial()`, `relmat()`	Structural dependence overview
Request an interval or profile a target	`confint()`, `profile_targets()`, `profile()`	Checking and using fitted models
Fit when some values are missing	`miss_control()`, `mi()`, `imputed()`	Handling missing data
Confirm an unfamiliar route is supported	`model-map` and `capability-and-limits`	What can I fit today?

The shortest useful workflow

For a continuous response whose expected value and residual variability may both depend on x, write one formula for mu and one for sigma, fit with the default TMB engine, and check the result before interpretation.

set.seed(20260721)
n <- 80
x <- rnorm(n)
dat <- data.frame(
  y = 0.4 + 0.7 * x + rnorm(n, sd = exp(-0.3 + 0.2 * x)),
  x = x
)

fit <- drmTMB(
  bf(y ~ x, sigma ~ x),
  family = gaussian(),
  data = dat
)

check_drm(fit)
#> <drm_check: 12 checks>
#> ok: 12; notes: 0; warnings: 0; errors: 0
#>                      check status
#>      optimizer_convergence     ok
#>           optimizer_budget     ok
#>           finite_objective     ok
#>      logsigma_clamp_active     ok
#>             fixed_gradient     ok
#>            sdreport_status     ok
#>  hessian_positive_definite     ok
#>     standard_errors_finite     ok
#>   standard_errors_inflated     ok
#>               dropped_rows     ok
#>             positive_scale     ok
#>   fixed_effect_design_size     ok
#>                                                                                  value
#>                                                                                      0
#>                                                iterations=11; function=21; gradient=12
#>                                                                                  91.76
#>                                                                                   <NA>
#>                                                 max=0.0002442; component=beta_sigma[1]
#>                                                                                     ok
#>                                                                                   TRUE
#>                                                                range=[0.07841,0.09395]
#>                                        n_inflated=0; max_se=0.09395; median_se=0.08237
#>                                                                     nobs=80; dropped=0
#>                                                                             min=0.5884
#>  total_mb=0.01372; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#>                                                                              message
#>                                                        nlminb convergence code is 0.
#>  Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                             Objective and log-likelihood are finite.
#>                                   The log(sigma) clamp is not active at the optimum.
#>     Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma[1].
#>                                              TMB::sdreport() completed successfully.
#>                                        sdreport reports a positive-definite Hessian.
#>                                         All fixed-effect standard errors are finite.
#>                   No fixed-effect standard error is inflated relative to the others.
#>                   No rows were dropped by model-frame or known-covariance filtering.
#>                                     All fitted scale values are finite and positive.
#>                          Dense fixed-effect design matrices are modest for this fit.
summary(fit)
#> <summary.drmTMB>
#> estimator: ML
#>                     estimate  std_error
#> mu:(Intercept)     0.4806752 0.08550257
#> mu:x               0.6758440 0.09395328
#> sigma:(Intercept) -0.2831471 0.07924446
#> sigma:x            0.1618198 0.07840772
#> Distributional, random-effect, scale, and correlation parameters:
#>                         component  dpar         term  estimate   minimum
#> fitted:sigma distributional-scale sigma fitted range 0.7710237 0.5883647
#>               maximum    scale
#> fitted:sigma 1.081276 response
#> logLik: -91.76
#> convergence: 0

Here mu is the expected response and sigma is the residual standard deviation. The sigma coefficients are on a log-SD scale, so use Which scale are you modelling? before treating them as variance or group-level effects.

Detailed function reference

When you need to…	Start with	Then go to
State which distributional parameters have formulas	`bf()` or `drm_formula()`	Formula grammar
Choose a likelihood for continuous, count, proportion, ordinal, or robust data	`gaussian()`, `nbinom2()`, `beta()`, `cumulative_logit()`, `student()`, or another family object	Choosing response families
Fit the model	`drmTMB()`	Getting started
Check convergence and boundary warnings	`check_drm()`	Improving convergence
Read fixed effects and their covariance	`summary()`, `coef(fit, dpar = "mu")`, `fixef()`, `vcov()`	Model workflow
Extract group-level or structured effects	`ranef()`, `structured_effects()`, `corpairs()`	Structural dependence overview
Predict one or more distributional parameters on a covariate grid	`predict_parameters()`	Distributional outputs and adequacy
Inspect fitted rows or residuals	`fitted()`, `residuals()`, `fitted_distribution()`	Distributional outputs and adequacy
Check distributional adequacy	`worm_plot()`, `qq_plot()`, `centile_chart()`, `exceedance()`	Distributional outputs and adequacy
Request an interval or inspect profile targets	`confint()`, `profile_targets()`, `profile()`	Checking and using fitted models
Simulate from the fitted model	`simulate()`	Checking and using fitted models

coef(fit) returns a named list, one numeric block per fitted distributional parameter. Select a block explicitly when writing analysis code:

coef(fit, dpar = "mu")
#> (Intercept)           x 
#>   0.4806752   0.6758440
coef(fit, dpar = "sigma")
#> (Intercept)           x 
#>  -0.2831471   0.1618198

grid <- data.frame(x = c(-1, 0, 1))
predict_parameters(
  fit,
  newdata = grid,
  dpar = c("mu", "sigma"),
  conf.int = TRUE
)
#>   row row_label  dpar            component     type   estimate  std.error
#> 1   1         1    mu             location response -0.1951688 0.11176154
#> 2   2         2    mu             location response  0.4806752 0.08550257
#> 3   3         3    mu             location response  1.1565192 0.14065979
#> 4   1         1 sigma distributional-scale response  0.6408455 0.07385634
#> 5   2         2 sigma distributional-scale response  0.7534090 0.05970348
#> 6   3         3 sigma distributional-scale response  0.8857440 0.09528542
#>     conf.low conf.high conf.level conf.status interval_source  x
#> 1 -0.4142174 0.0238798       0.95        wald            wald -1
#> 2  0.3130932 0.6482572       0.95        wald            wald  0
#> 3  0.8808311 1.4322073       0.95        wald            wald  1
#> 4  0.5112741 0.8032539       0.95        wald            wald -1
#> 5  0.6450268 0.8800023       0.95        wald            wald  0
#> 6  0.7173627 1.0936483       0.95        wald            wald  1

The returned prediction table records its interval provenance in conf.status and interval_source. A Wald interval for a fixed-effect prediction is not automatically an interval for a random-effect SD or a general claim of calibrated coverage.

Common routes

Counts, proportions, and robust continuous responses

Keep the model workflow above, but replace the family object and use the matching parameter names. For example, nbinom2() adds an NB2 dispersion parameter, beta-family models use mu and sigma, and zero-one beta also uses zoi and coi for exact boundary mass. The family guide gives the parameterization and the worked syntax for each route.

Two responses and residual correlation

For a bivariate Gaussian model, name the two location formulas and the residual-correlation formula explicitly. rho12 is residual coupling after the response-specific means and residual scales are modelled; it is not a group-level or phylogenetic correlation.

set.seed(20260722)
pair_dat <- data.frame(x = rnorm(120))
pair_dat$y1 <- 0.5 + 0.6 * pair_dat$x + rnorm(120, sd = 0.7)
pair_dat$y2 <- -0.2 + 0.4 * pair_dat$x + rnorm(120, sd = 0.8)

fit_biv <- drmTMB(
  bf(mu1 = y1 ~ x, mu2 = y2 ~ x, rho12 = ~ x),
  family = biv_gaussian(),
  data = pair_dat
)

rho12(fit_biv)
corpairs(fit_biv)

Read Changing residual coupling with rho12 before interpreting bivariate correlation layers.

Repeated groups, phylogeny, space, and relatedness

Random and structured terms change the estimand. A term such as (1 | id) models ordinary among-group variation; phylo(), spatial(), animal(), and relmat() use supplied structure. Their fitted and inference-backed surfaces are deliberately narrower than the fixed-effect surface, so check the capability and limits page for the exact family, parameter, and evidence boundary.

Keep the boundaries visible

engine = "tmb" is the current workflow. The Julia bridge is halted/deferred future work: it is not required to fit, diagnose, or report the R/TMB models shown here, and it is not an alternate engine for a new analysis.

This sheet helps you choose a first function, not a claim tier. Before a substantive analysis, confirm the supported syntax and inference boundary in What can I fit today? and What can I trust?.