drmTMB fits distributional regression models for one or
two responses. Start here if you have one response and want to ask
whether predictors change its average, its residual variability, or
both. This article fits a small Gaussian model, explains the output, and
shows what to check before interpretation.
The central workflow is simple: write one formula for each parameter you want to estimate, fit the simplest model that answers the question, and check the fit before interpreting it. If you instead have a phylogenetic tree, begin with Phylogenetic mixed models. If your response is an effect size with known sampling variance or covariance, begin with Mean effects and residual heterogeneity. The broader map, including other response types, is What can I fit today?.
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.
Install
This vignette describes drmTMB 0.7.1, the current
development version. drmTMB 0.7.0 is the first CRAN
release. Install it from CRAN:
install.packages("drmTMB")Or install the current development source from GitHub with
pak:
install.packages("pak")
pak::pak("itchyshin/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.
Fit your first model
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: 16 checks>
#> ok: 16; notes: 0; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> observations_per_parameter ok
#> fixed_effect_collinearity ok
#> dropped_rows ok
#> positive_scale ok
#> fixed_effect_design_size ok
#> value
#> 0
#> converged
#> iterations=23; function=35; gradient=23
#> 144.3
#> <NA>
#> max=0.000000001490; component=beta_sigma[1]
#> ok
#> TRUE
#> min_eig=37.38; cond=8.408
#> range=[0.07859,0.1484]
#> n_inflated=0; max_se=0.1484; reference_median=0.09549
#> n_obs=120; n_par=5; ratio=24.00
#> max_abs_r=0.05479
#> 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 convergence and uncertainty diagnostics are consistent with a proper interior optimum.
#> 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.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> No fixed-effect design column pair exceeds the collinearity note threshold.
#> 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().
Choose your next guide
The first fit above is the right next step when your response is continuous and you want to model its average and residual variation. Choose another guide only when your data require it.
| If your question is… | Continue with | Key distinction |
|---|---|---|
| Which response family fits continuous, count, proportion, robust, or zero-heavy data? | Choosing response families | A family determines what mu, sigma, and
any extra parameters mean. |
| Do related species, sites, animals, or a supplied relationship matrix share deviations? | Structural dependence overview | A structured deviation is not the same as residual variation or residual correlation. |
| Do two responses remain associated after their means and residual SDs are modelled? | Changing residual coupling with rho12 |
rho12 is residual coupling, not a correlation among
groups or species. |
| Do effect sizes come with known sampling variances or covariance? | Mean effects and residual heterogeneity | This is a Gaussian known-variance analysis, not a family for raw observations. |
| Is an estimate or interval safe to report? | Can I fit and report this model? | A model that fits can still have weak diagnostics or unsupported uncertainty. |
Use What can I fit today? only when you
need a compact route-and-limitation map. For any model, run
check_drm() before interpreting coefficients or
intervals.
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.
Check before interpreting
After fitting any model, run check_drm() before
interpreting the estimates:
check_drm(fit)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.
Keep correlation layers separate. A bivariate residual
rho12 describes how two responses vary together within an
observation after their means and residual SDs are modelled. A
random-effect correlation instead describes how group-level deviations
vary together. Phylogenetic and spatial structure are further kinds of
group-level pattern. They answer different biological questions and
should not be reported as the same quantity. If correlation is your
scientific question, continue with the structural dependence overview
before adding it to this introductory model.
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.