Check convergence and diagnostic flags for a drmTMB fit
Source:R/check.R, R/julia-diagnostics.R
check_drm.Rdcheck_drm() runs a compact set of model-fit diagnostics. It is intended as
a first-pass guardrail before interpreting distributional models, especially
fits with random effects, known sampling covariance, phylogenetic location
effects, or bivariate residual correlation rho12.
Usage
check_drm(object, ...)
# S3 method for class 'drmTMB'
check_drm(
object,
gradient_tolerance = 0.001,
rho_boundary = 0.98,
sd_boundary = 1e-04,
...
)
# S3 method for class 'drmTMB_julia'
check_drm(object, gradient_tolerance = 0.001, ...)Arguments
- object
A
drmTMBfit.- ...
Reserved for future diagnostic options.
- gradient_tolerance
Maximum absolute fixed-parameter gradient treated as acceptable.
- rho_boundary
Absolute residual or structured correlation value above which a bivariate Gaussian fit receives a warning.
- sd_boundary
Random-effect standard deviation below which a fit receives a warning that the variance component is near the lower boundary.
Value
A drm_check data frame. Its stable native reader columns are, in
order, check, status, value, and message, all character. status
uses only "ok", "note", "warning", and "error"; its "ok"
attribute is TRUE exactly when no row has status "warning" or
"error".
Details
The current checks cover optimizer convergence, finite objective values,
optimizer evaluation counts, fixed-parameter gradients including the largest
gradient component label, whether
TMB::sdreport() was computed, skipped, or failed, Hessian status from
TMB::sdreport(), the fixed-effect Hessian's minimum eigenvalue and
condition number at the optimum (hessian_conditioning; a comparable, but
not numerically identical, read of the same object underlying pdHess).
hessian_conditioning is computed from
TMB::sdreport()'s already-materialized fixed-effect covariance
(sdr$cov.fixed), never by evaluating the TMB C++ object directly, so it
is available for random-effect and REML fits (where the underlying
Hessian is Laplace-marginalized) and for a fit restored via saveRDS()/
readRDS(), and reports a note with a stated reason, not a bare NA,
only when sdr$cov.fixed itself is absent or incomplete (for example
drm_control(se = FALSE), or a failed sdreport()), finite fixed-effect standard errors, standard errors that
are finite but inflated relative to the others despite a positive-definite
Hessian (a weakly identified, near-flat direction such as a boundary
correlation), dropped rows,
positive scale parameters, random-effect standard deviations near the lower
boundary, bivariate residual-correlation rho12 values near the boundary,
Student-t nu boundary behaviour, skew-normal nu finite-value checks,
known sampling covariance summaries,
dense known-covariance storage scale, dense fixed-effect design size and
density, random-effect replication, and random-slope design variation. If a
univariate Gaussian fit includes one or more matched labelled
mu/sigma random-intercept covariance blocks, check_drm() also reports
group replication and whether either component is tiny relative to its
interpretation scale for each independent block. If a bivariate Gaussian fit
includes one or more matched same-response labelled mu/sigma
random-intercept covariance blocks, check_drm() reports one row per block.
If a bivariate Gaussian fit includes a matched
labelled mu1/mu2 random-intercept covariance block, check_drm() reports
group replication and whether either group-level SD is tiny relative to the
matching residual scale. For a matched labelled sigma1/sigma2 block, it
reports group replication and whether either log-sigma random-effect SD is
tiny. If a bivariate Gaussian fit includes an ordinary all-four q=4
mu1/mu2/sigma1/sigma2 block, it reports group replication, location
SDs relative to residual scales, log-sigma SDs, and whether any latent
correlation is near the boundary. If a bivariate Gaussian fit includes
matching mu1/mu2 phylogenetic location effects, check_drm() also reports
whether the fitted phylogenetic mean-mean correlation is near the boundary,
whether either phylogenetic SD is tiny relative to the matching residual
scale, and whether an ordinary group-level covariance block uses the same
grouping factor. Matching bivariate coordinate-spatial q=2, animal(), and
relmat() q=2 location effects receive the corresponding structured
replication, SD-ratio, and boundary-correlation diagnostics. If a bivariate
Gaussian fit includes a
phylogenetic, coordinate-spatial, animal-model, or relmat() q=4
mu1/mu2/sigma1/sigma2 block, it reports level replication, location
SDs relative to residual scales, log-sigma SDs, and whether any latent
structured correlation is near the boundary. If a univariate
Gaussian fit includes phylo(1 | species, tree = tree) or
phylo(1 + x | species, tree = tree) in mu, it reports species
replication, the fitted phylogenetic SDs, and whether the smallest
phylogenetic SD is tiny relative to the residual scale. If a univariate
Gaussian fit includes spatial(1 | site, coords = coords) or
spatial(1 + x | site, coords = coords) in mu, it reports site
replication, fitted coordinate range, the spatial SDs, and whether the
smallest spatial SD is tiny relative to the residual scale. If a Gaussian
fit includes sd_phylo(species) ~ x_species,
sd_phylo1(species) ~ x_species, or
sd_phylo2(species) ~ x_species, it reports species replication and the
fitted direct-SD surface range. If a univariate Gaussian fit used
drm_control(aggregate_gaussian = TRUE), it reports original rows,
aggregation cells, compression ratio, and largest cell size. If a fit was
stored with
drm_control(keep_tmb_object = FALSE), the
fixed-gradient check is reported as a note because the TMB
automatic-differentiation object is not available. If a fit used
drm_control(se = FALSE), the sdreport_status, Hessian, and
finite-standard-error checks are reported as notes. If sdreport() was
requested but failed, those rows are warnings.
Use check_drm() before interpreting coefficients, fitted values, or
response-scale quantities. A note records something to inspect, such as
dropped rows or a singly observed random-effect level. A warning means the
fitted model may still be useful but needs inspection before inference. An
error means at least one basic diagnostic failed. A Hessian or
sdreport() warning is therefore an inference and identifiability signal,
not automatic proof that fitted point estimates are unusable. For
programmatic checks, the returned object has attr(x, "ok") == TRUE only
when no rows have warning or error status.
For native drmTMB fits, this is a stable reader schema: the
drm_check data frame has, in order, the character columns check,
status, value, and message. status is one of "ok", "note",
"warning", or "error". Its "ok" attribute is TRUE exactly when
no row has status "warning" or "error".
Fits from engine = "julia"
A fit produced by the DRM.jl bridge (engine = "julia", class
drmTMB_julia) dispatches through this same generic and returns the same
four-column drm_check schema, but it reports a deliberately different
and much shorter set of rows. Almost every check listed above reads a
TMB object or a TMB::sdreport() that a bridge fit does not have, so a
clean engine = "julia" table is a narrower claim than a clean
engine = "tmb" table on the same model. The rows are:
engine_routeAlways a
note. Names the engine, the DRM.jl route, and the estimator ("engine=julia; route=<route>; estimator=<ML|REML>"), and its message lists the nativeengine = "tmb"checks that did not run. Being anote, it never flipsattr(x, "ok"): reporting which machinery ran is not itself a fault.optimizer_convergenceWhether DRM.jl reported convergence. This is DRM.jl's own convergence flag, not
nlminb()'sconvergencecode, and the two are not claimed to agree.fixed_gradientRoute-aware. Only some DRM.jl routes store a gradient callback, and the bridge omits the gradient entirely rather than sending a fabricated one on a route without it; on such a route this row is a
notenaming the route, stationarity is not scored, andattr(x, "ok")does not reflect it. When a gradient is present the row reportsmax|gradient|and the largest component. In both cases the value string carriessource=, naming which producer made the number in DRM.jl's owngrad_sourcevocabulary ("locscale","stored","forward","finite","none","unavailable") so that an exact analytic gradient, an automatic-differentiation gradient, and a finite difference accurate only to roughly1e-6relative do not print identically. The current bridge payload has exactly one gradient producer, the fit's own stored callback, so a gradient that crossed the bridge is reported assource=stored. When no gradient crossed, the bridge does not report which of the remaining five cases applies, so the row recordssource=unknown– deliberately not one of DRM.jl's six values – and points you at DRM.jl's owncheck_drm(fit)in Julia, whosegrad_sourcefield tells them apart. This is not TMB'ssdreport()gradient and the two engines are not claimed to agree numerically. A source outside that vocabulary, or one that contradicts the payload (a producer named when no gradient crossed, or"none"/"unavailable"alongside one that did), aborts rather than being printed.bridge_covarianceThe DRM.jl analogue of
sdreport_status: whether the fixed-effect covariance the bridge marshalled back is complete (ok), partial (some distributional parameters finite, others not, which some structured routes produce by design), or unavailable.TMB::sdreport()is not run for a bridge fit.bridge_standard_errorsThe DRM.jl analogue of
standard_errors_finite: whether every fixed-effect standard error implied by that covariance is a finite real number. A complete covariance can still carry a negative variance on the diagonal, whichbridge_covariancecannot see.
The gradient_tolerance argument applies as it does for a native fit; the
rho_boundary and sd_boundary arguments do not, because the
corresponding rows are not computed for a bridge fit.
Examples
set.seed(1)
dat <- data.frame(y = rnorm(40), x = rnorm(40))
fit <- drmTMB(drm_formula(y ~ x, sigma ~ x), data = dat)
check_drm(fit)
#> <drm_check: 15 checks>
#> ok: 15; notes: 0; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> convergence_status ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> hessian_conditioning ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> observations_per_parameter ok
#> dropped_rows ok
#> positive_scale ok
#> fixed_effect_design_size ok
#> value
#> 0
#> converged
#> iterations=12; function=20; gradient=12
#> 50.32
#> <NA>
#> max=0.0000000002422; component=beta_sigma[1]
#> ok
#> TRUE
#> min_eig=37.07; cond=2.308
#> range=[0.1129,0.1491]
#> n_inflated=0; max_se=0.1491; reference_median=0.1337
#> n_obs=40; n_par=4; ratio=10.00
#> nobs=40; dropped=0
#> min=0.7849
#> total_mb=0.007614; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#> message
#> nlminb convergence code is 0.
#> Optimizer convergence and uncertainty diagnostics are consistent with a proper interior optimum.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> The log(sigma) clamp is not active at the optimum.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_sigma[1].
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> Minimum eigenvalue and condition number of TMB's sdreport() fixed-effect covariance (sdr$cov.fixed), inverted. These are a genuinely different read of the fit's conditioning than TMB's internal pdHess flag -- comparable across fits, not claimed to be numerically identical to any raw TMB gradient or Hessian quantity. This fit's Hessian conditioning is within the requested threshold.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> Observations per estimated parameter are at or above the small-samplenote threshold (10).
#> No rows were dropped by model-frame or known-covariance filtering.
#> All fitted scale values are finite and positive.
#> Dense fixed-effect design matrices are modest for this fit.