Skip to contents

This page is for an applied ecology or evolution user in the first week with drmTMB. It walks one short path: fit → inventory targets → profile a random-effect SD → read profile.boundary → decide when not to trust the interval. It restates the default uncertainty story in ?confint.drmTMB; it does not claim nominal coverage on every route.

For the longer post-fit toolbox (check_drm(), Wald shortcuts, bootstrap, newdata profiles), read Checking and using fitted models. For which cells have coverage evidence, read Can I fit and report this model?.

Default recipe (read once)

  1. Fixed effects and other routine Wald-ready targets: start with confint(fit).
  2. Random-effect SDs and other direct variance components: prefer confint(..., method = "profile") after profile_targets() shows the row is profile-ready.
  3. Always read conf.status and profile.boundary on the returned table.
  4. A computable interval is not coverage certification.

Fit a small random-intercept model

The example is intentionally tiny and non-phylogenetic so it stays CRAN-safe and fast:

library(drmTMB)
#> 
#> Attaching package: 'drmTMB'
#> The following object is masked from 'package:base':
#> 
#>     beta

set.seed(13)
n_site <- 18
n_per_site <- 6
site <- factor(rep(seq_len(n_site), each = n_per_site))
site_effect <- rnorm(n_site, sd = 0.45)

fish_site <- data.frame(
  site = site,
  temperature = runif(n_site * n_per_site, -1.5, 1.5),
  habitat = factor(
    sample(c("reef", "kelp"), n_site * n_per_site, replace = TRUE)
  )
)

site_mu <- site_effect[as.integer(fish_site$site)]
fish_site$growth <- rnorm(
  nrow(fish_site),
  mean = 1 + 0.7 * fish_site$temperature +
    0.35 * (fish_site$habitat == "kelp") + site_mu,
  sd = 0.35
)

fit_site <- drmTMB(
  bf(growth ~ temperature + habitat + (1 | site), sigma ~ 1),
  family = gaussian(),
  data = fish_site
)

Run check_drm(fit_site) before interpreting estimates. A clean diagnostic table catches optimizer and Hessian problems; it does not certify coverage.

Inventory targets, then profile the RE-SD

profile_targets(fit_site)
#>                                 parm         target_class  dpar        term
#> 1               fixef:mu:(Intercept)         fixed-effect    mu (Intercept)
#> 2               fixef:mu:temperature         fixed-effect    mu temperature
#> 3               fixef:mu:habitatreef         fixed-effect    mu habitatreef
#> 4            fixef:sigma:(Intercept)         fixed-effect sigma (Intercept)
#> 5                              sigma distributional-scale sigma  (constant)
#> 6                   sd:mu:(1 | site)     random-effect-sd    mu  (1 | site)
#> 7 derived:total_variance_share(site)      derived-summary    mu  (1 | site)
#>   tmb_parameter index   estimate link_estimate    scale   transformation
#> 1       beta_mu     1  1.3588788     1.3588788     link linear_predictor
#> 2       beta_mu     2  0.6814230     0.6814230     link linear_predictor
#> 3       beta_mu     3 -0.2416811    -0.2416811     link linear_predictor
#> 4    beta_sigma     1 -0.9603030    -0.9603030     link linear_predictor
#> 5    beta_sigma     1  0.3827769    -0.9603030 response              exp
#> 6     log_sd_mu     1  0.3721615    -0.9884273 response              exp
#> 7          <NA>    NA  0.4859416            NA response   variance_ratio
#>   target_type profile_ready   profile_note
#> 1      direct          TRUE          ready
#> 2      direct          TRUE          ready
#> 3      direct          TRUE          ready
#> 4      direct          TRUE          ready
#> 5      direct          TRUE          ready
#> 6      direct          TRUE          ready
#> 7     derived         FALSE derived_target

Wald is the fast first pass for fixed effects and other Wald-ready rows:

confint(fit_site, parm = "fixed_effects")
#>                      parm level      lower       upper scale   transformation
#> 1    fixef:mu:(Intercept)  0.95  1.1553272  1.56243037  link linear_predictor
#> 2    fixef:mu:temperature  0.95  0.5880442  0.77480167  link linear_predictor
#> 3    fixef:mu:habitatreef  0.95 -0.3939463 -0.08941597  link linear_predictor
#> 4 fixef:sigma:(Intercept)  0.95 -1.1064327 -0.81417334  link linear_predictor
#>   tmb_parameter index method profile.engine conf.status profile.boundary
#> 1       beta_mu     1   wald           <NA>        wald               NA
#> 2       beta_mu     2   wald           <NA>        wald               NA
#> 3       beta_mu     3   wald           <NA>        wald               NA
#> 4    beta_sigma     1   wald           <NA>        wald               NA
#>   profile.message
#> 1            <NA>
#> 2            <NA>
#> 3            <NA>
#> 4            <NA>

For the site random-intercept SD, prefer the profile route and copy the exact target name from profile_targets():

ci_site <- confint(
  fit_site,
  parm = "sd:mu:(1 | site)",
  method = "profile"
)
ci_site
#>               parm level     lower    upper    scale transformation
#> 1 sd:mu:(1 | site)  0.95 0.2554993 0.561249 response            exp
#>   tmb_parameter index  method profile.engine conf.status profile.boundary
#> 1     log_sd_mu     1 profile       endpoint     profile            FALSE
#>   profile.message
#> 1              ok

Read profile.boundary before you trust the interval

Successful profile rows use conf.status = "profile". When the profile lands near a variance-component boundary, profile.boundary is TRUE and confint(method = "profile") warns with class drmTMB_profile_boundary_warning.

ci_site[, c("parm", "lower", "upper", "conf.status", "profile.boundary")]
#>               parm     lower    upper conf.status profile.boundary
#> 1 sd:mu:(1 | site) 0.2554993 0.561249     profile            FALSE

When not to trust a returned interval

  • profile.boundary == TRUE on a usable interval: treat the interval as indicative of scale, not as a calibrated level interval. The D-117 10-group random-effect SD gate measured conditional coverage as low as 0.1021 against a nominal 0.95 when that flag was set (and 0 of 89 in one cell where boundary hits are rare); the same behaviour appears in lme4 on the same seeds (see ?confint.drmTMB Boundary intervals). Prefer more groups when the design allows it.
  • Expect the SD itself to run low. At this design the point estimate was biased 8.3%-15.8% below truth – a property of maximum likelihood at few groups, not of this package. Unconditionally the interval still covered 0.9248 against a nominal 0.95 over 400,000 attempts, so the flagged rows are the bad case rather than the usual one.
  • REML = TRUE helps here, measurably, without fixing everything: on the same design and seeds it moved profile coverage from 0.9248 to 0.9463, roughly halved the SD’s downward bias, and made the misses nearly balanced. The boundary caveat above still applies under REML, and this is measured on this one design only – the default estimator is unchanged.
  • conf.status is profile_failed, clamp_limited, wald_at_boundary, newdata_required, or derived_interval_unavailable: there is no usable calibrated interval on that row — follow the action implied by the status.
  • The capability guide does not name the exact cell as inference-ready: report the point estimate only, or keep the interval experimental.
Next question Go here
Full post-fit toolbox and conf.status table Checking and using fitted models
Which family × effect rows have coverage evidence Can I fit and report this model?
Large tip counts and se_group_sd Working with large data
Reference detail for Wald, profile, and boundary warnings ?confint.drmTMB