First-week intervals: fit, profile, and boundary
Source:vignettes/first-week-intervals.Rmd
first-week-intervals.RmdThis 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)
- Fixed effects and other routine Wald-ready targets: start with
confint(fit). - Random-effect SDs and other direct variance components: prefer
confint(..., method = "profile")afterprofile_targets()shows the row is profile-ready. - Always read
conf.statusandprofile.boundaryon the returned table. - 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_targetWald 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 okRead 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 FALSEWhen not to trust a returned interval
-
profile.boundary == TRUEon a usable interval: treat the interval as indicative of scale, not as a calibratedlevelinterval. 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 inlme4on the same seeds (see?confint.drmTMBBoundary 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 = TRUEhelps 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.statusisprofile_failed,clamp_limited,wald_at_boundary,newdata_required, orderived_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.
What to read next
| 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 |