Distributional regression with drmTMB

Install the pre-CRAN release

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:

install.packages("pak")
pak::pak("itchyshin/drmTMB")

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.

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: 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.406861

For 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().

Learning path

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.

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.

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.

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.