Extract standard model-fit quantities

These methods expose drmTMB fits to standard base-R model summary and comparison helpers.

Usage

# S3 method for class 'drmTMB'
vcov(object, ..., type = "model", robust = FALSE)

# S3 method for class 'drmTMB'
logLik(object, ...)

# S3 method for class 'drmTMB'
AIC(object, ..., k = 2)

# S3 method for class 'drmTMB'
BIC(object, ...)

# S3 method for class 'drmTMB'
nobs(object, ...)

# S3 method for class 'drmTMB'
df.residual(object, ...)

# S3 method for class 'drmTMB'
deviance(object, ...)

Arguments

object

A drmTMB fit.

...

Reserved for future extractor options.

type

For vcov(), "model" (the default) for the model-based covariance matrix. "robust" is not implemented and errors rather than silently falling back to the model-based matrix.

robust

For vcov(), a TRUE/FALSE alias for type = "robust"/type = "model"; also errors when TRUE.

k

Numeric penalty per parameter for AIC(); the default is 2.

Value

logLik() returns an object of class "logLik". vcov() returns a numeric covariance matrix. nobs(), df.residual(), and deviance() return numeric scalars.

Details

logLik() returns a "logLik" object with df and nobs attributes so stats::AIC() and stats::BIC() use the fitted likelihood, optimized top-level parameter count, and fitted-row count consistently. nobs() returns the number of fitted rows after complete-case filtering. df.residual() returns nobs - df, where df is the number of optimized top-level parameters recorded in logLik(). deviance() returns -2 * logLik; for these likelihood-based distributional models this is an absolute negative twice log-likelihood value, not a saturated-model GLM deviance. vcov() returns the fixed-effect covariance matrix from TMB::sdreport() when the fit contains an sdreport object with pdHess = TRUE; the same object is available as fit$sdr and fit$sdreport. If sdreport() was skipped, failed, or returned pdHess = FALSE, Wald standard errors and Wald confidence intervals are unavailable while point estimates remain usable. vcov() intentionally does not include random-effect conditional modes or derived response-scale quantities.

Experimental MSPL fits

For a fit made with estimator = "mspl", vcov() behaves differently from an ordinary maximum-likelihood fit in two ways worth knowing before you use it.

First, the shape differs. It returns the full outer-parameter covariance after Laplace marginalisation — the fixed effects and the covariance parameters in their frozen Cholesky coordinates (log_sd_mu, and eta_cor_mu for a q = 2 block) — not the fixed-effect-only block an ordinary fit returns. Code that assumes vcov(fit) has the same dimension across estimators will break.

Second, it is a standard error, not an interval. The matrix inverts the Hessian of the unpenalized Laplace log-likelihood evaluated at the MSPL estimate: the penalty obtains a finite estimate but does not describe sampling variability. Do not form coef ± 1.96 * se. Kosmidis and Firth show that Wald intervals in this setting fail to cover regardless of the nominal level, and that the failure persists even for profile penalized-likelihood intervals; the mechanism is the finiteness of the penalized estimator and its standard error, not separation as such. A second reason applies here: the MSPL estimate maximises the penalized criterion, so the unpenalized score is not zero at it and the usual "evaluate at the maximum likelihood estimate" justification does not transfer. The size of that departure is recorded in fit$mspl$wald$unpenalized_gradient_max_abs.

Accordingly confint(), profile(), logLik(), AIC, BIC, and anova() deliberately error for MSPL fits. When the information matrix is not positive definite — which happens on exactly the strongly separated designs MSPL exists to handle — the standard errors are NA and a drmTMB_mspl_wald_unavailable warning is signalled, rather than a fabricated number.

Examples

set.seed(20260524)
n <- 36
x <- seq(-1.5, 1.5, length.out = n)
dat <- data.frame(
  y = 0.3 + 0.6 * x + rnorm(n, sd = 0.7),
  x = x
)
fit <- drmTMB(bf(y ~ x, sigma ~ 1), data = dat)

logLik(fit)
#> 'log Lik.' -26.64286 (df=3)
nobs(fit)
#> [1] 36
df.residual(fit)
#> [1] 33
deviance(fit)
#> [1] 53.28572
AIC(fit)
#> [1] 59.28572
BIC(fit)
#> [1] 64.03627
vcov(fit)
#>                   mu:(Intercept)          mu:x sigma:(Intercept)
#> mu:(Intercept)      7.145752e-03 -4.290051e-17     -2.479179e-17
#> mu:x               -4.290051e-17  9.012661e-03     -9.728111e-17
#> sigma:(Intercept)  -2.479179e-17 -9.728111e-17      1.388888e-02