drmTMB fits distributional regression models for one or
two responses. This introduction explains the model components, fits a
first model, and shows how to interpret it. If you only need to find a
function or tutorial, use the shorter Get
started function map.
The central workflow is simple: write one formula for each parameter you want to estimate, then check that the fitted model matches the symbolic model you had in mind. The broader implemented map lives in What can I fit today?, and family choice lives in Choosing response families.
The parameter names stay consistent across the site. Location is a
family-specific centre or location parameter; for Gaussian models it is
the expected response, but it need not be the unconditional response
mean. Scale describes residual variability; shape describes distribution
features beyond location and scale; and coscale means residual
correlation such as rho12 between two responses. A family
may expose only a subset of these components.
For example, an applied user might ask: do mean trait values and
residual variability change with an environmental predictor, after
accounting for repeated measures from the same site or species? In
drmTMB, that question is written as one formula for the
mean and one formula for the residual scale.
drmTMB 0.7.0 is the first CRAN-targeted release. At the
time this documentation was built, it had not yet been accepted by CRAN.
Install the current development source from GitHub with
pak:
After CRAN accepts the package, install the released version with
install.packages("drmTMB").
You need R 4.1.0 or newer and a working compiler toolchain because TMB models are compiled during installation. If installation fails while compiling C++, install the usual R build tools for your platform: Rtools on Windows, Xcode Command Line Tools on macOS, or the R development toolchain on Linux.
The core runtime dependencies are installed automatically by
pak: cli, Matrix,
TMB, and the compiled headers from RcppEigen
and TMB. The articles and development checks also use
optional packages such as glmmTMB, lme4,
MASS, metafor, knitr,
rmarkdown, testthat, withr, and
pkgdown.
Start with a Gaussian location-scale model when the response is continuous and the scientific question is about both the expected value and predictability. In the small example below, habitat and temperature affect mean growth, while habitat also changes residual variation:
set.seed(13)
n <- 120
dat <- data.frame(
habitat = factor(rep(c("forest", "grassland"), each = n / 2)),
temperature = rnorm(n)
)
mu <- 1 + 0.6 * (dat$habitat == "grassland") + 0.4 * dat$temperature
sigma <- exp(-0.5 + 0.45 * (dat$habitat == "grassland"))
dat$growth <- rnorm(n, mean = mu, sd = sigma)The fitted model uses one formula for mu and one formula
for sigma:
fit <- drmTMB(
drm_formula(growth ~ habitat + temperature, sigma ~ habitat),
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=23; function=35; gradient=23
#> 144.3
#> <NA>
#> max=0.0005430; component=beta_sigma[1]
#> ok
#> TRUE
#> range=[0.07859,0.1484]
#> n_inflated=0; max_se=0.1484; median_se=0.09549
#> nobs=120; dropped=0
#> min=0.7396
#> total_mb=0.02184; max_cols=3; largest=mu; largest_class=matrix; largest_density=0.8333
#> 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.Read the sigma coefficient as a log residual-SD
contrast. Exponentiating it gives an SD ratio; exponentiating twice the
coefficient gives a residual variance ratio:
sigma_habitat <- coef(fit, "sigma")["habitatgrassland"]
data.frame(
residual_sd_ratio = exp(sigma_habitat),
residual_variance_ratio = exp(2 * sigma_habitat)
)
#> residual_sd_ratio residual_variance_ratio
#> habitatgrassland 1.186112 1.406861For a fuller walkthrough of fitted means, residual SDs, and residual
variances, read When variance carries
signal, Part 1. Continue to Part
2 when predictors model a grouped or phylogenetic random-effect SD
through sd().
Start with the question you want to answer, not with the full list of features. The tutorials are arranged so that each article pairs symbolic equations, R syntax, fitted output, and interpretation.
The opening rows below follow the Get started menu: check the reporting boundary, choose a family, find the smallest useful function sequence, and then check the fitted model. The remaining rows route common modelling and diagnostic questions to the matching guide.
| If your question is… | Read this first | Main parameter |
|---|---|---|
| Can I fit this model, report its point estimate or interval, and what should I try if I cannot? | Can I fit and report this model? | reporting scope and named fallback |
| Which family fits my response — continuous, count, proportion, or robust? | Choosing response families | family-specific mu, sigma,
nu, zi, hu |
| I know my task but not which function or shortest workflow to use. | Function map and cheat sheet | task-to-function route and minimal workflow |
| After fitting, which checks and interval-status columns should I read? | Checking and using fitted models | check_drm(), profile_targets(),
conf.status |
| How do I profile a random-effect SD and know when not to trust the interval? | First-week intervals: fit, profile, and boundary | profile_targets(),
confint(method = "profile"),
profile.boundary |
| Am I modelling residual variation, group-level variation, or likelihood weights? | Which scale are you modelling? | sigma, sd(group),
weights |
| Do group or phylogenetic covariance rows answer a bivariate question beyond residual coupling? | What can I fit today? | corpairs(), corpair(),
sd1(group, level = "phylogenetic"),
sd2(group, level = "phylogenetic") |
| Do predictors change the mean and residual variability of one response? | When variance carries signal, Part 1 | mu, sigma |
| Do predictors change residual variability and the SD of a grouped or phylogenetic random effect? | When variance carries signal, Part 2 | mu, sigma, sd(group) |
| Do counts show extra-Poisson variation or structural zeros? | Count abundance and extra zeros | NB2 mu, sigma, zi |
Are proportions successes out of trials, continuous rates inside
(0, 1), or continuous rates with structural exact
boundaries? |
Proportions and success rates | beta-binomial, beta, and zero-one beta mu,
sigma, zoi, coi |
| Does a predictor change residual coupling between two responses? | Changing residual coupling with rho12 | rho12 |
| Do repeated individuals have correlated average responses after accounting for residual coupling? | Changing residual coupling with rho12 | corpairs(level = "group") |
| Do effect sizes have known sampling variances or covariance? (a specialist route, not a family choice for raw response data) | Mean effects and residual heterogeneity in meta-analysis | meta_V(V = V), sigma (implemented/tested;
tier unregistered, with no interval or coverage claim) |
| Do species trait means remain similar after accounting for shared ancestry, do sites share coordinate-structured deviations, or do individuals have known relatedness? | Structural dependence overview | animal(), phylo(), spatial(),
relmat() |
| I have more than one candidate fit — which one should I report? | Model selection with AIC and BIC | AIC(), BIC() |
| Is the fitted distribution itself adequate, not just the mean structure? | Distributional outputs and adequacy | worm_plot(), qq_plot(),
centile_chart() |
| The optimizer did not converge cleanly — what do I check? | Errors, warnings, and convergence | optimizer_convergence, fixed_gradient |
| My fit is slow or memory-heavy — how do I keep fitting at this scale? | Working with large data | keep_data, keep_model_frame,
se = FALSE, se_group_sd |
Here nu is a shape parameter, zi and
hu are zero-process parameters, and sd(group)
refers to a group-level random-effect standard deviation rather than the
residual scale sigma. Use the model guide “What can I fit
today?” when you need the longer status map before choosing syntax.
For a first applied analysis, fit the simplest model that answers the
question, run check_drm(), and then read the coefficient
table on the parameter scale used by the model. For example,
sigma coefficients are on a log scale in Gaussian
location-scale models, while rho12(fit) returns residual
correlations on the response scale.
For slope and variance-component questions, name the estimand before
reporting the number. A mu slope is an expected-response
effect, a sigma slope is a log residual-SD effect, a
random-slope SD is among-group variation in a reaction norm, and
sd(group) ~ x_group is a model for the SD of a group-level
mean effect. Those four quantities can all involve a predictor, but they
answer different biological questions.
After fitting any model, run check_drm() before
interpreting the estimates:
The diagnostic table checks convergence, gradients, Hessian status,
standard errors, dropped rows, scale values, random-effect replication,
and relevant parameter boundaries. Inspect a note; resolve
a warning or error before treating estimates
as stable. The errors, warnings, and
convergence guide explains what to try next.
The key rule is to keep correlation layers separate. A bivariate
residual rho12 is within-observation coupling after the two
means and residual SDs are modelled. A random-effect correlation is a
group-level quantity. Phylogenetic rows from
corpairs(..., level = "phylogenetic") are structured-effect
quantities, and the first coordinate-spatial mean-mean row is reported
by corpairs(..., level = "spatial") when matching labelled
spatial terms are fitted in mu1 and mu2. The
constant coordinate-spatial q=4 location-scale block is also fitted when
matching labelled spatial terms appear in all four
mu1/mu2/sigma1/sigma2
endpoints. Predictor-dependent spatial corpair() rows
remain planned. Those layers should not be reported as if they were the
same estimand.
Use Can I fit and report this model? for the current reporting boundary and named fallback; use the model map when you need syntax detail. Once this first fit is checked, continue with Checking and using fitted models or choose a scientific question from the Learning path table above.