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 calibrated level intervals.

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", "cluster"),
  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 = NA,
  ...
)

Arguments

object: A drmTMB fit.
parm: Optional character or integer vector selecting interval targets. NULL selects all direct targets for Wald intervals. Explicit profile intervals usually require target names or target-set shortcuts; for a Gaussian temporal AR1 or OU fit, NULL selects every mean-model coefficient. Temporal profile intervals admit mean coefficients only. 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. If method is omitted without newdata, every row is a Wald interval. Pass method = "profile" to profile random-effect SD rows.
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 the TMB profile step size unchanged and sets ytol to at least 0.5 above the requested likelihood cutoff, with a minimum of 2, when omitted. This gives the profile room to reach the cutoff; endpoint finding can still fail for other reasons. An explicit ytol remains in force. "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", Unix "multicore", or socket "cluster" (PSOCK, any platform including Windows). For profile intervals, targets or newdata rows are split across workers. For bootstrap intervals, refits are split across workers. Every backend returns the serial result; see drmTMB-parallel.
workers: Requested profile or bootstrap workers. If NULL and parallel is not "none", drmTMB uses about half the detected CPU cores. Parallel execution is capped at the number of jobs and at 10 workers. Keep BLAS and OpenMP threads at one while workers run; see drmTMB-parallel.
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). An ordinary (1 | g) SD whose log-scale Wald standard error exceeds 1 is flagged the same way even when its estimate is above sd_boundary; that cut-off is fixed.
small_sample_df: Small-sample reference distribution for the method = "wald" interval. "location" (the default) uses a t-quantile with df = g - 1 for location-axis (mu, mu1, mu2) structured random-effect SD targets (phylo, spatial, animal, relmat) when a unique group count g is available. Other targets use the normal quantile. This changes interval width, not its centre. The retained engine coverage check concerns four q2 phylogenetic and relatedness-matrix mu1:x/mu2:x slope SD cells at g = 8, with a post-hoc check across g = 8, 16, 32; the default's broader application does not establish nominal coverage for every target. Dispersion-axis (sigma, sigma1, sigma2) targets stay at the normal quantile by default. Their coverage depends on the model: a replicated pure sigma intercept over-covered in a pilot, whereas some sigma one-slope intercept targets under-covered at g = 8. Fixed effects and plain labelled covariance blocks are also unchanged by the default. "none" uses the normal quantile for every target and recovers the earlier behaviour. "group" widens every resolvable structured or labelled-covariance SD target on either axis; this broader choice is not covered by the q2 calibration.
bias_correct: Small-sample centre adjustment for the method = "wald" interval. "location" (the default) adds log(g / (g - 1)) to the log-SD estimate of each location-axis structured random-effect target with a unique group count g, before back-transforming; other targets are unshifted. In the retained q2 phylogenetic and relatedness-matrix mu1:x/mu2:x slope SD cells at g = 8 and a post-hoc cross-g check, this shift approximated measured ML log-SD shrinkage. Combined with the t width, it brought coverage near 0.95 in those cells. Applying the same formula to another route is an implemented calculation, not evidence of its bias or interval coverage. "none" leaves every centre unshifted; "group" shifts every resolvable structured or labelled-covariance SD target on either axis. Dispersion targets remain unshifted by default because their coverage is regime dependent; the magnitude is also not calibrated for plain labelled covariance blocks. The shift changes the centre only; it does not alter the fitted model or a profile endpoint.
bootstrap_re_form: re.form passed through to stats::simulate() when method = "bootstrap", with the same lme4 meanings as simulate.drmTMB(). The default, NA, redraws random effects for every replicate (marginal simulation), matching what a parametric bootstrap of the fitted data-generating process requires. Passing NULL 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. Earlier drmTMB versions gave NULL and NA the reverse meaning (#1487).
...: 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. The returned parm values use the fully-qualified target names used by profile_targets(), such as fixef:sigma:z, even when the input parm used a compact label such as sigma:z. Normalize names before joining the result back to user-supplied compact labels. 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. For a plain matrix view of the lower and upper columns, use as.matrix.drm_confint().

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 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.
Gaussian temporal AR1, OU, and homogeneous Toeplitz mean effects: use confint(fit, parm = "mu:<coefficient>", method = "profile") (or omit parm for every mean coefficient). Temporal Wald output remains restricted or deferred, and temporal variance, persistence, decay, bootstrap, and newdata intervals are unavailable. Read check_drm(fit): a non-positive-definite fitted Hessian does not prevent a finite profile endpoint, but it makes the fitted likelihood locally irregular and the interval is not coverage-calibrated.
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 model and parameter combination in Capabilities and limits before reporting its intervals as coverage-validated.

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

Wald intervals for ordinary `(1 | g)` SDs

The default confint(fit) is Wald for every row, including random-effect SD rows. The Wald SD interval undercovers near the boundary. In the sweep-oct16 t9 design (Gaussian location-scale with (1 | g), n_g = 10, 300 replicates per cell) at a true SD of 0.1, Wald coverage was 0.36 with 10 groups and 0.42 with 30 groups, with about half of the replicates returning an infinite upper limit; profile coverage on the same replicates was 0.97 and 0.98 there. At a true SD of 0.5, Wald coverage was 0.97 to 0.98 and profile coverage 0.93 to 0.957. When a SD row is labelled wald_at_boundary (see Boundary intervals), use confint(fit, parm = ..., method = "profile") for that row. Profile is slower than Wald, and it is not a repair at a boundary: it still warns when it reaches one. That sweep does not calibrate intervals for other ordinary-SD routes or different study designs.

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 is unreliable at a variance-component boundary. That warning recommends method = "profile". For an ordinary (1 | g) SD this label is not limited to an estimate at zero: it is also applied when the log-scale Wald standard error exceeds 1. At the default 95% level, this makes the ratio of the upper to lower Wald endpoint exceed about 50. This width does not establish coverage. Treat any such wald_at_boundary SD row as unreliable.

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 simulation measured coverage at 0.1021, 0.2387, and 0.8683 against a nominal 0.95, and 0 of 89 in a fourth scenario 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 simulation 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 model class: eleven of twelve retained ML estimates fell below truth in the sibling nbinom2 provider scenarios (fit-level one-sided sign test p = 0.0032; p = 0.0625 at scenario 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. Only point-estimate recovery has been checked for these routes, and no simulation study 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 default correction is motivated by mixed-model variance-component bias but its log(g / (g - 1)) log-SD magnitude is empirical, not a general consequence of REML. The retained calibration covers q2 phylogenetic and relatedness-matrix mu1:x/mu2:x slope SD cells (engine check at g = 8, post-hoc checks at g = 8, 16, 32); neither other structured routes nor the broader "group" scope inherit that coverage result. In the elementary independent Gaussian variance example, the familiar variance factor is g/(g-1), equivalent to only 0.5 * log(g/(g-1)) on the SD log scale. The cited theory supplies context and direction, not a derivation of this package's full shift. Relevant references:

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.