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, prints and reads its
summary(), checks the fit with check_drm(),
and points to the next guide for your data.
The workflow has three steps: 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?.
This vignette describes drmTMB 0.7.1, the current
development version; 0.7.0 is on CRAN. Installation instructions,
including the compiler toolchain TMB needs, are on the package home page.
Parameter names used in these guides
Each distributional parameter has one formula name and one plain-English name. The guides use these names consistently:
| Formula name | Name in these guides | What it describes |
|---|---|---|
mu (mu1, mu2 for two
responses) |
location | Where the response distribution sits. For a Gaussian response,
mu is the expected response, and its coefficients read like
lm() or lmer() coefficients. |
sigma (sigma1, sigma2) |
scale | How far observations spread around mu. For a Gaussian
response, sigma is the residual standard deviation (SD),
modelled on the log scale. |
nu and other family-specific names |
shape | Features beyond location and scale, such as Student-t tail weight
(nu). Gaussian models have none. |
zi, hu
|
zero inflation, hurdle | The probability of an extra zero (zi) or of a zero from
a separate hurdle process (hu). |
rho12 |
residual correlation (coscale) | How two responses vary together within an observation after their
mu and sigma models are accounted for. Only
bivariate models have it. |
A family exposes only the parameters it defines: a Gaussian model has
mu and sigma, whereas a Poisson model has
mu, plus zi when you model zero inflation.
If you already use lme4 or glmmTMB
The mean formula is written exactly as in lme4 or
glmmTMB, including random-effect terms such as
(1 | site). Each additional parameter gets its own formula,
named by the parameter. The sigma formula plays the role of
glmmTMB’s dispformula for a Gaussian
model:
# glmmTMB
glmmTMB::glmmTMB(y ~ x + (1 | site), dispformula = ~ x, data = dat)
# drmTMB: the same model
drmTMB(bf(y ~ x + (1 | site), sigma ~ x), family = gaussian(), data = dat)bf() is a short alias for drm_formula();
the name and the one-formula-per-parameter layout follow
brms::bf() (Bürkner 2017). For worked side-by-side fits
against lm(), lme4, glmmTMB,
metafor, and other packages, see Comparing drmTMB with other
packages.
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 its predictability. In the example below, habitat and temperature change mean growth, and habitat also changes the residual SD:
growth_i ~ Normal(mu_i, sigma_i^2)
mu_i = beta_0 + beta_1 grassland_i + beta_2 temperature_i
log(sigma_i) = gamma_0 + gamma_1 grassland_iHere i indexes observations, grassland_i is
1 for grassland and 0 for forest, the beta coefficients
belong to the mu formula, and the gamma
coefficients belong to the sigma formula. The simulation
sets beta = (1, 0.6, 0.4) and
gamma = (-0.5, 0.45), so you can check whether the fit
recovers them:
set.seed(1)
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 model has one formula for mu (the unnamed
growth ~ ... formula) and one for sigma:
fit <- drmTMB(
drm_formula(growth ~ habitat + temperature, sigma ~ habitat),
family = gaussian(),
data = dat
)Check the fit
Run check_drm() before reading any estimate. It reports
one row per check: optimizer convergence, the size of the final
gradient, whether the Hessian is positive definite, whether standard
errors are finite, whether rows were dropped, whether fitted scales are
positive, and, for models with random effects, whether groups are
replicated enough and whether parameters sit on a boundary:
check_drm(fit)
#> <drm_check: 18 checks>
#> ok: 18; 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
#> mle_existence ok
#> observations_per_parameter ok
#> fixed_effect_collinearity ok
#> dropped_rows ok
#> positive_scale ok
#> residual_outliers ok
#> fixed_effect_design_size ok
#> value
#> 0
#> converged
#> iterations=22; function=37; gradient=22
#> 140.7
#> <NA>
#> max=0.0000000009213; component=beta_sigma[1]
#> ok
#> TRUE
#> min_eig=40.18; cond=7.820
#> range=[0.07811,0.1508]
#> n_inflated=0; max_se=0.1508; reference_median=0.09129
#> finite
#> n_obs=120; n_par=5; ratio=24.00
#> max_abs_r=0.002386
#> nobs=120; dropped=0
#> min=0.6156
#> n_beyond_3=0/120
#> 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.
#> The screened fixed-effect geometry is consistent with a finite maximum-likelihood estimate.
#> 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.
#> No excess of quantile residuals beyond |3|.
#> Dense fixed-effect design matrices are modest for this fit.A clean fit looks like the header line above: every check is
ok, with zero notes, warnings, and errors. Read the
message column of any note; resolve a
warning or error before treating the estimates
as stable. The errors, warnings, and
convergence guide explains what to try next.
Read the summary
summary(fit)
#> <summary.drmTMB>
#> estimator: ML
#> estimate std_error
#> mu:(Intercept) 0.9703459 0.07991617
#> mu:habitatgrassland 0.5467444 0.15077857
#> mu:temperature 0.4134070 0.07810551
#> sigma:(Intercept) -0.4851669 0.09128765
#> sigma:habitatgrassland 0.4776546 0.12910105
#> Distributional, random-effect, scale, and correlation parameters:
#> component dpar term estimate minimum
#> fitted:sigma distributional-scale sigma fitted range 0.8040551 0.6155945
#> maximum scale
#> fitted:sigma 0.9925158 response
#> Derived summaries: This fit has a sigma predictor, a non-log link, or a
#> known-dispersion override, so a single scalar residual variance is not defined.
#> Random-effect SD estimates remain available in $sdpars.
#> logLik: -140.7
#> convergence_status: converged
#> convergence: 0The first table has one row per coefficient, labelled by parameter:
- Rows beginning with
mu:are thebetacoefficients, on the response scale.mu:temperatureis the change in expected growth per unit of temperature, read exactly like anlm()slope. - Rows beginning with
sigma:are thegammacoefficients, on the log-SD scale.sigma:habitatgrasslandis the log of the ratio of grassland to forest residual SD.
The second table, fitted:sigma, gives the mean, minimum,
and maximum of the fitted residual SDs across observations, on the
response scale. logLik is the maximised log-likelihood, and
convergence: 0 repeats the optimizer code that
check_drm() already inspected.
Compare the estimates with the simulated values: the mu
rows should be close to 1, 0.6, and 0.4, and the sigma rows
close to -0.5 and 0.45, each within about two standard errors. A
different seed gives different estimates; the truth is a target for the
method on average, not for any single data set.
Because sigma coefficients are on the log scale,
exponentiate them before reporting. exp(gamma_1) is the
ratio of grassland to forest residual SD, and
exp(2 * gamma_1) is the ratio of residual variances:
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.612288 2.599474For a fuller walkthrough of fitted means, residual SDs, residual
variances, and intervals, read When
variance carries signal, Part 1. Continue to Part 2 when a predictor should
change the SD of a random effect. That model uses an sd()
formula: sd(site) ~ x says that the among-site SD of the
(1 | site) intercepts depends on the site-level predictor
x.
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 correlated after their means and residual SDs are modelled? | Residual correlation with rho12 |
rho12 is a residual correlation, 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. |
How does my lme4, glmmTMB, or
metafor model translate? |
Comparing drmTMB with other packages | Side-by-side fits of the same model in both packages. |
Use What can I fit today? when you need a compact route-and-limitation map.
Name the quantity before reporting it
Before reporting a number, name the estimand, the quantity the number
estimates. Four quantities can all involve a predictor x,
but they answer different biological questions:
- A
muslope, fromy ~ x, is the change in the expected response per unit ofx. - A
sigmaslope, fromsigma ~ x, is the change in log residual SD per unit ofx. - A random-slope SD, from
(0 + x | group)in themuformula, is how much thexslope (a reaction norm) varies among groups. - An
sd()slope, fromsd(group) ~ x_group, is how a group-level predictor changes the SD of the group intercepts.
Keep residual and group-level correlations separate
A bivariate residual correlation rho12 describes how two
responses vary together within an observation after their
mu and sigma models are accounted for.
A random-effect correlation describes something else: how group-level deviations, such as site intercepts and site slopes, vary together. Phylogenetic and spatial correlations are further group-level correlations, structured by a tree or by coordinates. Report each as its own quantity. If correlation is your scientific question, continue with the structural dependence overview before adding it to this introductory model.
Next steps
Use Can I fit and report this model? for the current reporting boundary for each model, and the model map for syntax detail. Once this first fit is checked, continue with Checking and using fitted models, or choose a scientific question from the “Choose your next guide” table above.
Working in Julia? The DRModels.jl first-model
page is the Julia starting point, and its R and
Julia vocabulary page maps drmTMB syntax to Julia
syntax. Its first example uses different simulated data from this
article, so compare syntax rather than numbers. Julia is optional for
drmTMB users.
References
Bürkner, P.-C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1-28. https://doi.org/10.18637/jss.v080.i01