Confidence intervals for fitted model parameters

confint() returns confidence intervals for a fitted drmTMB model. Wald intervals are fast and are returned for fixed-effect coefficients and direct response-scale parameter targets by default. Direct Wald targets include constant residual-scale, random-effect standard-deviation, random-effect correlation, and constant residual-correlation rows when the fitted TMB parameter and TMB::sdreport() covariance are available. Correlation Wald intervals are computed on the fitted TMB correlation-link scale, equivalent to a guarded Fisher z/atanh transform, and then returned on the correlation scale. Bootstrap intervals simulate and refit direct targets. For positive scale and SD targets, percentile endpoints are taken on the fitted log scale before back-transforming to the response scale. Profile-likelihood intervals are slower because nuisance parameters are re-optimized; this first public profile path supports explicit fixed-effect, constant distributional-scale, random-effect standard-deviation, random-effect correlation, bivariate phylogenetic q=2 mean-mean correlation, block-diagonal bivariate phylogenetic mu1/mu2 and sigma1/sigma2 correlations, and constant residual-correlation targets. For predictor-dependent scale, residual-correlation, or currently supported corpair() formulae, supply newdata with parm = "sigma", parm = "rho12", or the fitted corpair(...) dpar to profile the fitted response-scale value for each supplied row. Row-specific rho12 profile intervals from a regression rho12 formula are computable this way, but no simulation has measured their coverage; treat them as interval-feasible, not as a calibrated level interval (tracked as issue #802).

Usage

# S3 method for class 'drmTMB'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = c("wald", "profile", "bootstrap"),
  newdata = NULL,
  trace = FALSE,
  profile_precision = c("default", "fast"),
  profile_maxit = NULL,
  profile_engine = c("auto", "endpoint", "tmbprofile"),
  profile_endpoint_max_eval = NULL,
  R = 199L,
  seed = NULL,
  parallel = c("none", "multicore"),
  workers = NULL,
  refit_control = NULL,
  sd_boundary = 1e-04,
  rho_boundary = 0.98,
  small_sample_df = c("location", "none", "group"),
  bias_correct = c("location", "none", "group"),
  bootstrap_re_form = NULL,
  ...
)

Arguments

object: A drmTMB fit.
parm: Optional character or integer vector selecting interval targets. NULL selects all direct Wald-ready targets for Wald intervals. Profile intervals require explicit target names or target-set shortcuts. Supported shortcuts are "fixed_effects", "random_effects", "variance_components", and "correlations".
level: Confidence level.
method: Interval method: "wald", "profile", or "bootstrap". If newdata is supplied and method is omitted, method = "profile" is used.
newdata: Optional data frame for response-scale profile intervals for predictor-dependent sigma, sigma1, sigma2, rho12, or fitted corpair() values. Each row is profiled separately by profiling its fixed-effect linear predictor and then transforming the interval to the response scale.
trace: Logical; passed to TMB::tmbprofile() when the tmbprofile profile engine is used.
profile_precision: Profile-control shortcut. "default" leaves TMB::tmbprofile() controls unchanged. "fast" supplies ystep = 0.5 and ytol = 2 unless the caller supplies those controls in ..., giving a quicker first-pass profile for long variance-component or correlation targets.
profile_maxit: Optional positive whole number passed to TMB::tmbprofile() as maxit when method = "profile". Use this as a per-target adaptive-step budget for long or exploratory profile runs.
profile_engine: Profile engine for direct fitted-object targets. "auto" uses a scalar endpoint solver for direct scale, SD, and correlation targets when no TMB::tmbprofile() controls are supplied, and otherwise uses TMB::tmbprofile(). "endpoint" requires the scalar endpoint solver, while "tmbprofile" preserves the previous full-profile route for comparison and debugging.
profile_endpoint_max_eval: Optional positive whole number limiting constrained endpoint evaluations per endpoint side when the scalar endpoint engine is used. This is a diagnostic escape hatch for long variance-component or correlation profiles; when the budget is reached the row is returned with conf.status = "profile_failed" and missing endpoints.
R: Number of parametric-bootstrap refits when method = "bootstrap".
seed: Optional seed for bootstrap simulation.
parallel: Profile or bootstrap backend: "none" or Unix "multicore". For profile intervals, targets or newdata rows are split across workers. For bootstrap intervals, refits are split across workers.
workers: Requested profile or bootstrap workers. If NULL and parallel = "multicore", drmTMB uses about half the detected CPU cores. Multicore execution is capped at the number of jobs and at 10 workers.
refit_control: Optional drm_control() object used for bootstrap refits. The default skips TMB::sdreport() and drops the TMB object because bootstrap intervals use refit point estimates.
sd_boundary, rho_boundary: Boundary thresholds used by method = "wald" to flag intervals where the symmetric Wald interval is unreliable: a variance-component standard deviation within sd_boundary of zero, or a correlation within rho_boundary of +/-1, is returned with conf.status = "wald_at_boundary" and a warning pointing to method = "profile". Defaults match check_drm() (1e-4, 0.98).
small_sample_df: Small-sample reference distribution for the method = "wald" interval. "location" (the default) references a t-quantile with df = g - 1 – a group-based, Satterthwaite-style between-group degrees of freedom – for each location-axis (mu, mu1, mu2) structured random-effect SD target (phylo, spatial, animal, relmat) with a resolvable group count g, and keeps the ordinary normal quantile z = qnorm((1 + level) / 2) for every other target. This widens only the location-axis structured-RE SD intervals, which under-cover at small g under the normal quantile; the t-quantile lifts their coverage and converges back to z as g grows. Dispersion (sigma, sigma1, sigma2) structured SD targets are deliberately left at the normal quantile because they already over-cover; non-structured and fixed-effect targets are never widened; and a plain labelled covariance block such as (1 + x | p | id) is also left at the normal quantile by default, because the correction magnitude was calibrated for structured blocks only. "none" is the opt-out: it uses the normal quantile for every target and is byte-identical to the pre-default behaviour. "group" widens every resolvable SD target – structured and labelled-covariance, location and dispersion axis alike; use it only when you deliberately want to t-reference the dispersion axis or a labelled covariance block too. Note that the t-quantile only widens the interval; it does not move its centre. The residual small-g gap is ML shrinkage bias in the variance-component point estimate (a job for REML), not a quantile problem, so the centre shift is handled by bias_correct.
bias_correct: Small-sample point-estimate bias correction for the method = "wald" interval. "location" (the default) adds a simulation-calibrated shift log(g / (g - 1)) to the log-scale point estimate – before back-transforming – of each location-axis (mu, mu1, mu2) structured random-effect SD target with a resolvable group count g (the same targets small_sample_df = "location" widens), and leaves every other target unshifted. The ML estimate of a structured-RE variance-component SD is biased low on the log scale by about log(g / (g - 1)) – a downward shrinkage that REML, or simply a larger g, removes asymptotically; this shift moves the interval centre up to counter that shrinkage. Used together with small_sample_df = "location", the centre shift and the t(df = g - 1) width act independently, giving exp((log(sigma_hat) + log(g / (g - 1))) +/- qt(p, g - 1) * se_log) and lifting small-g coverage of location-axis structured-RE SD targets to nominal. "none" is the opt-out: it leaves every centre at the ML estimate and is byte-identical to the pre-default behaviour. "group" shifts every resolvable SD target – structured and labelled-covariance, location and dispersion axis alike. The default deliberately excludes the dispersion axis (sigma SD intervals already over-cover, so the upward shift would push them further conservative) and plain labelled covariance blocks (whose correction magnitude is not yet simulation-calibrated).
bootstrap_re_form: re.form passed through to stats::simulate() when method = "bootstrap". The default, NULL, redraws random effects for every replicate (marginal simulation), matching what a parametric bootstrap of the fitted data-generating process requires. Passing NA instead holds random effects fixed at their conditional-mode estimate; this is a documented escape hatch for random-effect structures whose marginal draw is not yet implemented, and it makes the resulting intervals anticonservative because between-group variability is never resampled – it is not a neutral alternative to the default.
...: Additional arguments passed to TMB::tmbprofile() when method = "profile" and the tmbprofile profile engine is used. drmTMB supplies the profiled obj, name, lincomb, and trace arguments internally; set the profile target with parm.

Value

A data frame with columns parm, level, lower, upper, scale, transformation, tmb_parameter, index, method, and profile.engine, conf.status, profile.boundary, and profile.message. Successful rows currently use conf.status = "wald", "profile", or "bootstrap". Failed numeric profile rows use "profile_failed" with missing endpoints; profile rows mark intervals that land near a lower SD boundary or correlation boundary. Bootstrap interval results carry a "bootstrap.diagnostics" attribute with one diagnostic row per refit and target, including refit convergence, target availability, draw use, and the refit message.

Details

Modelled missing-predictor bootstrap is not yet implemented: it requires joint predictor simulation and refits that preserve the imputation model and missingness pattern. Such requests fail before simulation. Positive predictor scales use log-scale Wald intervals transformed by exp(); predictor mixture probabilities use logit-scale intervals transformed by plogis(). Public covariance and summary standard errors remain on the reported natural scale.

Target names follow the profile target namespace. For fixed effects, use names such as "fixef:mu:x", "fixef:sigma:(Intercept)", or "fixef:rho12:w". Compact coefficient labels from summary(fit), such as "mu:x", are also accepted. Random-effect SD intervals are reported on the SD scale, and random-effect correlation intervals are reported on the correlation scale. For bivariate Gaussian and exact bivariate lognormal fits with constant residual correlation, parm = "rho12" profiles the residual correlation and reports the interval on the response correlation scale. For fits with constant sigma, sigma1, or sigma2, parm = "sigma" and friends report response-scale intervals.

The fastest routine route is confint(fit), which uses Wald intervals for fixed effects and direct response-scale targets. For long phylogenetic, spatial, animal-model, or relatedness fits, profile only the needed variance-component or correlation rows with the default profile_engine = "auto" first; direct scalar scale, SD, and correlation targets use the endpoint engine when no full-profile controls are supplied. Use profile_engine = "tmbprofile" or profile_precision = "fast" when you want the previous full-curve TMB::tmbprofile() route for comparison, diagnostics, or control tuning.

Default uncertainty story

Use this recipe for ordinary first-week inference; it restates measured behaviour and does not claim nominal coverage on every route.

Fixed effects and other routine Wald-ready targets: start with confint(fit) (method = "wald"). It is the fastest fitted-object route when TMB::sdreport() succeeded.
Random-effect SDs and other direct variance-component targets: prefer confint(fit, parm = ..., method = "profile") after profile_targets(fit) shows the row is profile-ready. Profile re-optimizes nuisance parameters and is slower than Wald.
Always read conf.status and profile.boundary on the returned table. A usable profile interval that lands on a variance-component boundary warns with class drmTMB_profile_boundary_warning; treat that interval as indicative of scale, not as a calibrated level interval (see Boundary intervals below).
A computable interval is not coverage certification. Check the exact cell in Capabilities and limits before reporting a route as inference-ready.

For a short applied walkthrough, see the vignette First-week intervals: fit, profile, and boundary.

Boundary intervals

Every interval method is unreliable when a variance component approaches zero or a correlation approaches +/-1, and each warns about its own case.

method = "wald" flags a row with conf.status = "wald_at_boundary" and warns (class drmTMB_wald_boundary_warning), because a symmetric Wald interval undercovers under boundary (chi-square-mixture) inference. That warning recommends method = "profile".

method = "profile" sets profile.boundary = TRUE and warns (class drmTMB_profile_boundary_warning) when it returns a usable interval that reaches a boundary. A profile interval is not a repair for a boundary. Conditional on that flag, a seeded 10-group random-effect SD gate measured coverage at 0.1021, 0.2387, and 0.8683 against a nominal 0.95, and 0 of 89 in a fourth cell where boundary hits are rare (0.09% of fits). The mechanism is a maximum likelihood random-effect SD biased 8.3%-15.8% low at small group counts, which anchors the interval low so that it misses from above.

Unconditionally – across all fits, flagged or not – the same gate measured 0.9248 against a nominal 0.95 over 400,000 attempts. A flagged interval is the bad case, not the typical one, and the flag is returned in the table so you can tell which case you are in.

Refitting the same design with REML = TRUE improved this without repairing it: over 400,000 paired replicates, profile coverage moved from 0.9248 to 0.9463 against a nominal 0.95, the SD point estimate's downward bias roughly halved (pooled -10.9% under maximum likelihood to -4.6% under REML), and the upper-to-lower miss asymmetry fell from 5.7:1 to 2.0:1. Coverage conditional on the boundary flag improved but stayed well below nominal (0.74 to 0.83), so a flagged interval remains the bad case under either estimator. The default estimator is unchanged; this measurement covers the one design below only.

This is a property of profile intervals near a variance boundary rather than of drmTMB: lme4::lmer on the same data-generating process and seeds agreed on boundary incidence for 4000/4000 replicates and matched the conditional coverage to four decimal places. Treat a flagged interval as indicative of scale, not as a calibrated level interval, and prefer more groups where the design allows it.

These figures come from one measured design: Gaussian, a single random intercept on the mean, 10 groups, 4 or 10 observations per group, maximum likelihood. They describe that corner, and are not a general statement about every family, provider, or group count.

method = "bootstrap" flags a row with conf.status = "bootstrap_at_boundary" and warns (class drmTMB_bootstrap_boundary_warning) when at least 5% of the retained resamples land on the target's bound. Resampling does not repair a boundary: a percentile interval whose draws pile up at zero is reporting the constraint rather than the sampling distribution. Measured on three Gaussian random-intercept fits at R = 200, a true SD of 0 put 43% of draws on the bound and a true SD of 0.25 put 5% there with a lower endpoint of exactly zero, while a true SD of 0.9 put none there. The flag needs at least 20 retained draws to fire, because a share computed from a handful of resamples is noise; bootstrap.n reports how many were retained.

Rows with conf.status = "profile_failed" or "clamp_limited" also carry profile.boundary = TRUE, but return missing endpoints and are not warned about separately; read their conf.status and profile.message instead.

check_drm() does not assess any of this. It reads the fit, so a target whose point estimate sits well clear of zero passes every fit-level check while its interval is still flagged here. check_drm() emits an interval_reliability_scope note saying so; read conf.status before reporting an interval.

Profiling a structured `sigma` random-effect SD

A profile of a sigma-axis random-effect SD under phylo(), animal(), relmat(), spatial(), or phylo_interaction() is computable for zero_one_beta(), nbinom2(), and zi_nbinom2(), but the maximum likelihood point estimate it is built around is biased low for this class of cell: eleven of twelve retained ML estimates fell below truth in the sibling nbinom2 provider cells (fit-level one-sided sign test p = 0.0032; p = 0.0625 at cell level, once the shared data-generating process and random-number stream are respected). Two corrections that exist elsewhere for this bias do not reach these routes: bias_correct shifts only the method = "wald" centre and never a profile endpoint, and native scale-side REML is unavailable here because drm_validate_reml_spec() admits only Gaussian and binomial models. The profile interval therefore inherits the bias of its centre uncorrected. These routes carry point-fit recovery evidence only, and no seeded campaign has yet measured the coverage of a profile interval built on them, so treat the endpoints as a computed profile rather than calibrated inference.

Ordinal cutpoints

For a cumulative_logit() fit, the public ordinal targets are the ordered latent-logistic cutpoints "ordinal:cutpoint:<label>", not the internal "ordinal:theta_ord:<label>" coordinates. Use confint(fit, parm = "ordinal:cutpoint:<label>", method = "profile") to obtain a pointwise likelihood-ratio interval on the cutpoint scale. The constrained engine fixes the cumulative cutpoint while preserving strict ordering of every threshold; it does not change the fitted likelihood. These are not simultaneous bands or category-probability intervals, and they are not yet calibrated coverage claims. Raw theta_ord entries remain visible in profile_targets() only as internal diagnostics and reject interval methods. The fixed-effect coefficients on mu (and on any other fitted distributional parameter) remain available through the usual confint(fit, parm = "fixed_effects") or confint(fit, parm = "fixef:mu:x") calls.

References

The small-sample corrections applied by small_sample_df and bias_correct (whether the default "location" scope or the broader "group" scope) are motivated by established mixed-model theory but are simulation-calibrated, not derived: the log(g / (g - 1)) SD-scale centre shift is about twice the leading-order REML SD correction (0.5 * log(g / (g - 1))), matching the larger ML shrinkage measured on these structured/bivariate cells (their effective df is well below g - 1). The magnitude's authority is the per-model-class simulation in docs/design/219-structured-re-small-sample-bias-correction.md, not a single source. Relevant references:

Restricted maximum likelihood and the variance-component bias that motivates the centre shift (REML debiases the variance by g/(g-1); the shift here is on the SD log scale and ~2x the leading-order REML SD term):

Patterson, H. D., & Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. Biometrika, 58(3), 545-554. doi:10.1093/biomet/58.3.545

Harville, D. A. (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358), 320-338. doi:10.1080/01621459.1977.10480998

Searle, S. R., Casella, G., & McCulloch, C. E. (1992). Variance Components. New York: Wiley. doi:10.1002/9780470316856

Analytic and penalized first-order bias reduction of maximum likelihood estimates (the general framework an additive log-scale shift instantiates):

Cox, D. R., & Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society, Series B, 30(2), 248-265. doi:10.1111/j.2517-6161.1968.tb00724.x

Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27-38. doi:10.1093/biomet/80.1.27

The t-quantile / between-group effective degrees of freedom used by small_sample_df:

Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110-114. doi:10.2307/3002019

Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53(3), 983-997. doi:10.2307/2533558

Boundary regime (a variance component at or near zero), where neither the t-width nor the centre shift restores nominal coverage:

Self, S. G., & Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605-610. doi:10.1080/01621459.1987.10478472

Stram, D. O., & Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171-1177. doi:10.2307/2533455

Parametric-bootstrap bias correction of mixed-model variance components, and the general delicacy of bootstrap bias estimation (a single-level parametric bootstrap does not recover the centre bias for these targets at small g):

Kubokawa, T., & Nagashima, B. (2012). Parametric bootstrap methods for bias correction in linear mixed models. Journal of Multivariate Analysis, 106, 1-16. doi:10.1016/j.jmva.2012.01.011

Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. New York: Chapman & Hall.

Small-sample variance-component (repeatability) interval coverage in ecology and evolution:

Wolak, M. E., Fairbairn, D. J., & Paulsen, Y. R. (2012). Guidelines for estimating repeatability. Methods in Ecology and Evolution, 3(1), 129-137. doi:10.1111/j.2041-210X.2011.00125.x

Examples

dat <- data.frame(y = c(0.2, 0.5, 1.1, 1.4), x = c(-1, 0, 1, 2))
fit <- drmTMB(bf(y ~ x, sigma ~ 1), data = dat)
confint(fit)
#>                      parm level       lower      upper    scale
#> 1    fixef:mu:(Intercept)  0.95  0.51798632  0.6620137     link
#> 2              fixef:mu:x  0.95  0.36120108  0.4787989     link
#> 3 fixef:sigma:(Intercept)  0.95 -3.39479062 -2.0088873     link
#> 4                   sigma  0.95  0.03354758  0.1341379 response
#>     transformation tmb_parameter index method profile.engine conf.status
#> 1 linear_predictor       beta_mu     1   wald           <NA>        wald
#> 2 linear_predictor       beta_mu     2   wald           <NA>        wald
#> 3 linear_predictor    beta_sigma     1   wald           <NA>        wald
#> 4              exp    beta_sigma     1   wald           <NA>        wald
#>   profile.boundary profile.message
#> 1               NA            <NA>
#> 2               NA            <NA>
#> 3               NA            <NA>
#> 4               NA            <NA>
confint(fit, parm = "variance_components")
#>    parm level      lower     upper    scale transformation tmb_parameter index
#> 1 sigma  0.95 0.03354758 0.1341379 response            exp    beta_sigma     1
#>   method profile.engine conf.status profile.boundary profile.message
#> 1   wald           <NA>        wald               NA            <NA>
confint(fit, parm = "sigma", method = "profile")
#>    parm level      lower     upper    scale transformation tmb_parameter index
#> 1 sigma  0.95 0.03817528 0.1645079 response            exp    beta_sigma     1
#>    method profile.engine conf.status profile.boundary profile.message
#> 1 profile       endpoint     profile            FALSE              ok
# Use the full-profile engine when you need the older tmbprofile route:
# confint(fit, parm = "sigma", method = "profile",
#   profile_engine = "tmbprofile", profile_precision = "fast")
# Direct-target parametric bootstrap is available when refit cost is worth it:
# confint(fit, parm = "sigma", method = "bootstrap", R = 99)
# Bootstrap intervals for positive scale and SD targets use link-scale
# percentiles before back-transforming to the response scale.