Model fitting and post-fit tools

Status: Reference

Mirrors drmTMB's Model fitting and post-fit tools (23 in drmTMB). The fitting verb is drm; the post-fit accessors below cover coefficients, fitted scale / correlation, random-effect estimates, predictions, simulation, inference, and convergence diagnostics.

Fitting

DRModels.drm Function

drm(formula::DrmFormula, family; data) -> DrmFit

Fit a distributional regression model by maximum likelihood. A formula bundle has one linear predictor per distributional parameter:

fit = drm(bf(@formula(y ~ x1), @formula(sigma ~ x1)), Gaussian(); data = dat)

Univariate Gaussian fits support fixed effects plus the structured-effect markers documented under phylo, spatial, animal, relmat, and meta_V.

ids = … labels the rows of the relatedness matrix A of an animal(1 | g) term, so A is matched to the data by label instead of by order of first appearance (see animal). Likewise ids = … labels the rows of coords of a spatial(1 | site) term (or coords is a Dict keyed by site label); see spatial.

algorithm and sparse — solver selection

algorithm (default :auto) and sparse choose how the model is fit:

• :auto (default) — uses the all-node sparse L-BFGS route for the Gaussian phylogenetic-mean cell (phylo(1 | g) on mean with sigma ~ 1). For phylogenetic location-scale-scale models (sd(species, phylogenetic) ~ z), :auto selects the O(p) sparse augmented GMRF engine when G > 500 species and the dense scaled-covariance engine for smaller trees. Other Gaussian cells keep their cell-specific default fitters.

• :gls, :lbfgs — legacy dense leaf-covariance fitters for the Gaussian phylogenetic-mean cell and aliases for the usual default fitters elsewhere.

• :em — force the all-node sparse conjugate-EM route for the Gaussian phylogenetic-mean cell. It reaches the same MLE as the dense GLS fit (same β, residual σ, and marginal logLik) via closed-form E/M steps with exact O(p) Takahashi traces. Any other model cell raises a clear ArgumentError. The EM path has no coefficient vcov (the M-steps are closed-form), so vcov(fit) is filled with NaNs. Refit with :gls for dense-fit Wald inference. re_sd(fit) reports the EM's Brownian phylo SD σ_phy (a different scale from the GLS fit's correlation-matrix σ_s).

• :sparse — force the verified sparse structured-Gaussian route where one is implemented, including the two-structured phylo + animal/relmat sparse path and the O(p) sparse phylogenetic location-scale-scale engine.

• :sparse_lbfgs — force the default all-node sparse L-BFGS route for the Gaussian phylogenetic-mean cell, or the O(p) augmented-state Takahashi selected inverse engine for phylogenetic location-scale-scale models.

• sparse = true — keyword alias to select sparse solvers (e.g. for whole-tree phylogenetic LSS or two-structured Gaussian models).

fit = drm(bf(@formula(y ~ x + phylo(1 | sp)), @formula(sigma ~ 1)), Gaussian();
          data = dat, tree = tree)

Bivariate Gaussian fits use BivariateDrmFormula; with no structured marker they fit the residual rho12 model, and with shared phylo(1 | group) markers on mu1, mu2, sigma1, and sigma2 they route to the verified q=4 phylogenetic engine.

method — ML (default) or REML

method (default :ML) selects the estimator. :REML is opt-in and is implemented for: (a) the fixed-effect Gaussian location–scale cell, (b) a single Gaussian mean random intercept (1 | g) using the Woodbury identity, (c) Location–Scale–Scale (LSS) models (sd(g) ~ z, sd(species, phylogenetic) ~ z, and multi-component LSS models), and (d) the bivariate q=4 PLSM Laplace engine (reml_q4), and (e) the Gaussian phylo(1 | g)-on-sigma location–scale routes (scale-only, separate, and phylo_coupled = true), where one joint Laplace approximation integrates the phylo effects and BOTH the mean and scale fixed effects — the restricted likelihood native drmTMB maximises.

Ordinary σ-RE, random slopes, and non-Gaussian REML stay rejected. REML likelihoods are not comparable across fixed-effect structures.

marginal: how a random effect on sigma is integrated

marginal = :LA (the default) leaves every route unchanged; a random intercept on sigma, sigma ~ 1 + (1 | g), is then integrated by 32-node Gauss–Hermite quadrature. marginal = :Laplace integrates it by the Laplace approximation instead, which is what drmTMB computes for this model; the fit is tagged fit.marginal === :Laplace. :Laplace requires a fixed-effect mean and maximum likelihood and is refused on every other Gaussian model.

The two names are not synonyms. :LA means "the default integrator for this route", which is not always the Laplace approximation (here it is Gauss–Hermite quadrature; wherever the Gaussian marginal is closed-form it is exact). :Laplace always means the Laplace approximation drmTMB computes. Bootstrap refits of a :Laplace fit use :Laplace too.

Missing response handling

Incomplete responses (missing or NaN in y) are supported under the observed-rows pattern (matching response = "include" in the R bridge). Inf and -Inf are not missing values: they raise ArgumentError naming the response column and the row. A non-finite predictor on an observed row, or a non-finite offset or coordinate, raises the same way, naming that argument. A non-finite predictor on a row whose response is missing or NaN is dropped with that row. predict on newdata that does not carry the response checks every row, so a NaN predictor there is an error rather than a NaN prediction. A returned fit never has converged = true when its log-likelihood is non-finite or a failed-fit sentinel (≤ -1e15, including -1e18 and -floatmax) or when any coefficient is non-finite. For Location-Scale-Scale models, the group index and scale design Z_g are parameterised over all G levels, while the likelihood is evaluated on observed rows. With a (… | g) random effect or meta_V(...) and no structured term, rows with a missing response are dropped, with a warning that gives the count.

Missing predictor handling

Rows with a missing predictor value are dropped before the fit, with a warning that gives the count, so the fit uses the complete cases (R's na.omit, which is what drmTMB does). nobs(fit) counts the rows used. A relmat/animal/spatial term whose level order the drop would change is refused instead, because its matrix is matched to the levels by position. An mi(x) joint formula models its missing predictors rather than dropping them.

source

drm(formula::BivariateDrmFormula, Gaussian(); data, tree = nothing,
    K = nothing, A = nothing, coords = nothing, g_tol = 1e-8,
    q4_g_tol = 1e-3, q4_iterations = 300, q4_n_newton = 40,
    q4_vcov = true, method = :ML) -> DrmFit

Fit a bivariate Gaussian distributional regression model.

With no structured-effect marker, this is the residual-correlation model: mu1, mu2, sigma1, sigma2, and residual rho12 each have their own fixed effect formula.

With matching structured markers on mu1 and mu2 only, this routes to the q=2 exact-Gaussian point-fit cell, method = :ML (default) or method = :REML (Patterson–Thompson; marginalises beta_mu1/beta_mu2 only — see reml_q2.jl). The q=2 route currently accepts phylo(1 | group) with tree, relmat(1 | group) with K, animal(1 | group) with A, or spatial(1 | group) with coords; it requires complete responses, identical mean fixed-effect designs, and intercept-only sigma1, sigma2, and rho12. spatial(1 | group) supplies a covariance exactly as relmat/animal do: the exponential kernel exp(-d / rho) at a FIXED range rho (keyword spatial_range, else the mean off-diagonal pairwise distance), the same rule the q=4 structured route uses. It does NOT estimate the range jointly the way the univariate spatial route (_fit_spatial_gaussian) does, so the two are not the same model and must not be compared as one. The range in force is recorded as fit.ranef.spatial_range.

Row k of coords is the k-th level of group in the order the levels first appear in data, not sorted label order, which is the same contract K and A carry on this route. Only the row COUNT is checked, so coordinates supplied in the wrong order fit a different spatial structure without error.

With the same phylo(1 | group) marker on all four location/scale predictors (mu1, mu2, sigma1, and sigma2) and a supplied tree, this routes to the verified q=4 phylogenetic location-scale engine. The residual rho12 formula remains the residual correlation; the group-level 4×4 covariance Σ_a is stored as fit.ranef.Sigma_a, with axes (:mu1, :mu2, :sigma1, :sigma2). Population parameter prediction skips the internal :phylocov coefficient block.

fit = drm(
    bf(mu1 = @formula(y1 ~ x + phylo(1 | species)),
       mu2 = @formula(y2 ~ x + phylo(1 | species)),
       sigma1 = @formula(sigma1 ~ 1 + phylo(1 | species)),
       sigma2 = @formula(sigma2 ~ 1 + phylo(1 | species)),
       rho12 = @formula(rho12 ~ 1)),
    Gaussian();
    data = dat,
    tree = phy,
)

source

drm(f::BivariateDrmFormula, ::Student; data, g_tol = 1e-8, method = :ML)

Fit a bivariate Student-t scatter-correlation model — drmTMB's biv_student().

mu1/mu2 are locations (identity link), sigma1/sigma2 are scale parameters (log link) — not standard deviations; for ν > 2 the marginal SD = σ·sqrt(ν/(ν−2)) — rho12 is the scatter correlation (guarded atanh link), and nu is the shared degrees of freedom on the logm2 scale ν = 2 + exp(η), so ν > 2 and the variance is finite.

nu is shared across both responses by construction, and zero rho12 does not mean independent margins at finite ν.

Structured markers (phylo/relmat/animal/spatial) — NOT implemented

This route stays residual-only: no random effects, no structured markers, and no REML. This is a deliberate rejection, not a missing port. The bivariate Gaussian structured routes are exact (q=2) or a verified Laplace approximation with a hand-derived, analytically differentiated per-leaf density (q=4); both depend on the conditional response being Gaussian, which a Student-t density is not. A Gaussian group-level random effect under a heavy-tailed, per-row scale-mixture likelihood has no closed-form marginal, and there is no verified non-Gaussian engine in this codebase to reuse for it — building one for this route alone would be inventing a parallel, unverified numerical design rather than mirroring the established one, on exactly the family (correlated structured scale) this project has already gotten a scale convention silently wrong once. drmTMB has the same limit, and that is the load-bearing half of this rejection, so it is written down reproducibly rather than asserted. In native drmTMB, the fixed-effect fit below is accepted and adding the structured marker is refused:

library(drmTMB); library(ape)
tr <- compute.brlen(stree(8, type = "balanced"), method = "Grafen")
tr$tip.label <- paste0("s", 1:8)
d <- data.frame(y1 = rnorm(40), y2 = rnorm(40), x = rnorm(40),
                sp = factor(rep(paste0("s", 1:8), each = 5)))
# fixed effects only -> ACCEPTED
drmTMB(bf(mu1 = y1 ~ x, mu2 = y2 ~ x, sigma1 = ~1, sigma2 = ~1,
          nu = ~1, rho12 = ~1), family = biv_student(), data = d)
# add a structured marker -> REFUSED
drmTMB(bf(mu1 = y1 ~ x + phylo(1 | sp, tree = tr), mu2 = y2 ~ x,
          sigma1 = ~1, sigma2 = ~1, nu = ~1, rho12 = ~1),
       family = biv_student(), data = d)
#> `biv_student()` currently allows fixed-effect formulas only; random and
#>  structured effects are deferred.

So there is no reference implementation on either side of the port. For a parity goal that matters more than the numerical argument above: a parity gap cannot be closed against a capability the reference package does not have, and implementing one unilaterally would mean inventing the answer this port exists to mirror.

f = bf(mu1 = @formula(y1 ~ x), mu2 = @formula(y2 ~ x),
       sigma1 = @formula(sigma1 ~ 1), sigma2 = @formula(sigma2 ~ 1),
       nu = @formula(nu ~ 1), rho12 = @formula(rho12 ~ 1))
fit = drm(f, Student(); data = dat)
2 + exp(coef(fit, :nu)[1])   # estimated shared degrees of freedom

source

drm(f::BivariateDrmFormula, ::LogNormal; data, tree = nothing, K = nothing,
    A = nothing, coords = nothing, spatial_range = nothing, g_tol = 1e-8,
    q4_g_tol = 1e-3, q4_iterations = 300, q4_n_newton = 40, q4_vcov = true,
    method = :ML)

Fit a bivariate lognormal residual-correlation model — drmTMB's biv_lognormal().

Both responses must be strictly positive. log(y1), log(y2) are modelled as bivariate normal, so mu1/mu2 are means on the log scale (identity link) and rho12 is the log-residual correlation — not the Pearson correlation of the raw responses. sigma1/sigma2 are log-scale SDs (log link).

Structured markers (phylo/relmat/animal/spatial)

phylo(1 | group), relmat(1 | group), animal(1 | group), and spatial(1 | group) markers are supported through the same delegation that already gives this family its residual fit: log(Y) is exactly bivariate Gaussian, so a structured call here literally runs drm(f, Gaussian(); data = log.(data), …) and shifts the reported log-likelihood by the (parameter-free) Jacobian. There is no separate lognormal structured engine to build or verify — the q=2 exact-Gaussian route (matching markers on mu1/mu2 only) and the q=4 sparse-Laplace PLSM route (matching markers on mu1, mu2, sigma1, and sigma2) are exactly the ones documented under drm(::BivariateDrmFormula, ::Gaussian), run on logged data. Because the Jacobian does not depend on any parameter, theta/vcov/ranef and the fitted objective's gradient carry over untouched; only the log-likelihood value (and everything derived from it: aic/bic/deviance) shifts.

Scale of a structured SD. The group-level covariance (fit.ranef.Sigma_a for q=4, the random-intercept variance for q=2) is on the log-response scale for the mu1/mu2 axes, and on the log(SD of log-response) scale for the sigma1/sigma2 axes — two log transforms deep for the scale axes, mirroring the fixed-effect sigma1/sigma2 convention above. ranef(fit) and vc(fit) inherit this without change because they are the untouched Gaussian output on log(y).

Matching drmTMB, method = :REML is not implemented for this family: the residual-only cell has no random effects to integrate out, and REML is also unavailable for the structured models; use method = :ML (the default).

d = (; y1 = exp.(0.4 .+ randn(200)), y2 = exp.(0.1 .+ randn(200)), x = randn(200))
fit = drm(bf(mu1 = @formula(y1 ~ x), mu2 = @formula(y2 ~ x),
             sigma1 = @formula(sigma1 ~ 1), sigma2 = @formula(sigma2 ~ 1),
             rho12 = @formula(rho12 ~ 1)),
          LogNormal(); data = d)
coef(fit)      # mu1/mu2 on the log scale
corpairs(fit)  # log-residual rho12

source

drm(f::BivariateDrmFormula, fams::Tuple; data, kwargs...)

Cross-family bivariate fit — two responses from different families coupled by a latent-scale correlation rho. The Julia twin of drmTMB's family = c(gaussian(), poisson()) route.

fit = drm(bf(mu1 = @formula(y1 ~ x), mu2 = @formula(y2 ~ x),
             sigma1 = @formula(sigma1 ~ 1), sigma2 = @formula(sigma2 ~ 1)),
          (Gaussian(), Poisson()); data = dat)
fit.rho_latent            # the latent-scale correlation
mf_summary(fit)

Each axis takes its own mu and (where the family has one) sigma formula. sigma is ignored for dispersionless families (Poisson, Binomial).

rho12 is not a formula here

In the Gaussian bivariate route rho12 ~ x models a per-observation residual correlation. In the cross-family route the correlation is a latent scalar — the two responses live on different scales, so there is no common residual to correlate per observation. Supplying rho12 is an error rather than being silently ignored.

Returns the fit_mixed_family result; see mf_summary, mf_coef. This route is experimental — see the cross_family_latent capability row.

source

DRModels.DrmFit Type

DrmFit

A fitted distributional regression model. Accessors: coef, coef(fit, :mu), vcov, loglik, nobs, fixef.

source

Optimization controls

DRMControl selects fitting controls; control_info reports what the fit used. The returned DRMControlReport is an Experimental interface, as defined in API stability.

DRModels.DRMControl Type

DRMControl(; maxit = nothing, start = nothing, multi_start = 1, se = true,
             seed = 1, jitter = 0.5, optimizer = :default)

Optimiser and inference settings for drm, passed as one keyword: drm(f, family; data, control = DRMControl(...)). The field names follow drmTMB's drm_control().

Field	Meaning
maxit	Iteration cap of the outer optimiser (nothing = the route's default cap).
start	Starting values: a vector on the internal parameter scale, in coef(fit) order of a default fit. fit.coefnames groups its labels by parameter block. nothing = the route's own start.
multi_start	Number of starts. 1 = a single start. With k > 1, the first start is start (or the route's own start), and starts 2:k add reproducible perturbations to it. The fit with the lowest objective is kept.
se	false skips the standard errors. vcov(fit) is then all NaN.
seed	Seed of the perturbations for multi_start (a DRModels.jl extension).
jitter	Perturbation size for multi_start. Start j > 1 is x0 .+ jitter .* max.(abs.(x0), 1) .* z, where z is standard normal (a DRModels.jl extension).
optimizer	:default preserves the route's current method and restart policy. :lbfgs or :bfgs runs that Optim.jl method directly on qualified routes.

Which routes honour which field

• se = false works on every route. Fixed-effect fits of every univariate family skip the Hessian. Routes that already had an se keyword skip it too. On any other route, the Hessian is still computed and then discarded. control_info(fit).se tells you which case applied (:skipped or :discarded).

• maxit, start and multi_start are honoured on the fixed-effect routes of every univariate family. This includes the main Gaussian location-scale route. A model with random, structured (phylo, relmat, spatial, …) or meta_V terms, or a bivariate model, refuses these fields with an ArgumentError before fitting. A route that does not honour a field it was given never ignores it silently: it throws an ArgumentError that names the field.

• optimizer = :lbfgs and optimizer = :bfgs work on the fixed-effect ML routes above and Gaussian ML with one mean random intercept (1 | g). That random-intercept route uses its own conditional two-start restart; maxit, start, multi_start, seed and jitter remain unavailable there. Other latent, structured, bivariate and REML routes refuse an explicit optimizer. Gaussian selection also requires algorithm = :auto and no alternate integration, imputation or profiling options. This control selects the primary fit only. Fit-based profile_result, confint(...; method = :profile) and bootstrap_* refuse an explicitly selected fit until their refits can forward the method. Wald intervals remain available when standard errors were computed.

The default DRMControl() changes nothing: a fit with it is identical to a fit without control.

See control_info for the requested and actual method, attempted stages, native stopping settings, the kept start and the se outcome.

Example

f = bf(@formula(dist ~ 1 + speed), @formula(sigma ~ 1 + speed))
fit = drm(f, Gaussian(); data = cars, control = DRMControl(se = false))
all(isnan, vcov(fit))                                  # no standard errors

fit5 = drm(f, Gamma(); data = dat, control = DRMControl(multi_start = 5, maxit = 2000))
control_info(fit5).start_logliks                       # one log-likelihood per start

source

DRModels.control_info Function

control_info(fit::DrmFit) -> Union{DRMControlReport, Nothing}

What drm(...; control, weights) did for this fit (a DRMControlReport), or nothing when neither keyword was given. Use it to see the log-likelihood reached from each multi_start start (start_logliks), which start was kept (kept_start), whether the starts agree (starts_agree), and whether se = false skipped the Hessian or discarded it (se).

Example

fit = drm(f, Gamma(); data = dat, control = DRMControl(multi_start = 5))
info = control_info(fit)
info.start_logliks, info.kept_start, info.starts_agree

source

DRModels.DRMControlReport Type

DRMControlReport

What drm(...; control, weights) did, returned by control_info.

Fields:

• control::DRMControl: the settings used.

• start_logliks::Vector{Float64}: Optim.jl's reported log-likelihood from each start, in start order (NaN for a start that failed). Empty when the optimiser settings were not used.

• kept_start::Int: index of the start that was kept (0 when not used).

• starts_agree::Bool: true when every successful start reached the kept log-likelihood to within 1e-6 * (1 + |loglik|).

• se::Symbol: :computed, :skipped (the Hessian was not computed) or :discarded (the route computed it, and the wrapper dropped it because se = false).

• weights: the prior weights, or nothing.

• requested_optimizer / actual_optimizer: the request and the method run on the selected stage (:unrecorded when this route does not expose its optimizer to this report).

• optimizer_runs: immutable per-attempt records in execution order with method, stage, start, native iterations/convergence/objective/gradient residual, stopping settings and finite/error status. The objective and gradient residual on fixed-effect routes are Optim.jl's reported values; Gaussian random-intercept entries recompute both at the returned minimizer.

• selected_optimizer_stage: the stage that supplied the returned fit.

• native_controls: a per-fit snapshot of maxit, start, multi_start, seed, jitter and se.

source

Coefficients and (co)variance

DRModels.fixef Function

fixef(fit) -> Vector{Pair}

Fixed-effect coefficients per distributional parameter, with their names.

source

DRModels.re_sd Function

re_sd(fit; scale = :native, tree = nothing) -> Dict{Symbol,Float64}

Estimated random-effect standard deviations, keyed by grouping factor. A mean-axis random intercept (y ~ x + (1|g)) or independent slope (y ~ x + (0+x|g)) is keyed by the bare group name and is on the mean's linear-predictor (link) scale, which is the response scale only for an identity link such as Gaussian's. A scale-axis random intercept (sigma ~ 1 + (1|g)) is keyed <group>_logsigma because that SD lives on the log-σ scale — the two are NOT directly comparable, and the suffix keeps them distinct so a side-by-side read is not silently mixing scales. A correlated random intercept+slope block ((1 + x | g)) is keyed <group>_intercept and <group>_slope, consistent with sqrt.(diag(vc(fit)[:g])); the correlation itself is not returned here — use vc for the full 2×2 covariance.

scale (#732, twin drmTMB#1272)

For a phylo(1 | g) grouping fitted on the default raw branch-length covariance (fit.phylo_scale === :covariance: the sparse Gaussian-mean route on a non-ultrametric tree and all non-Gaussian Laplace/GLMM phylo routes; the sparse Gaussian-mean route on an ultrametric tree is on the correlation scale since #979), re_sd's default scale = :native returns σ on that raw branch-length scale (tip variance = the tree's height h). drmTMB instead reports the phylogenetic SD on the tip-correlation scale (ape::vcv(tree, corr = TRUE), tip variance 1 regardless of h); the two quantities differ by the exact factor sqrt(h) (sd_drmTMB == re_sd(fit)[:g] * sqrt(phylo_tree_height(augmented_phy(tree))), confirmed against drmTMB 0.7.1 to within optimiser tolerance on both a Gaussian and a non-Gaussian (CumulativeLogit) phylo fit — see test/test_twin_gap_732.jl).

Pass scale = :drmtmb and the SAME tree (or newick string) given to drm(...) to get drmTMB's number directly instead of doing that conversion by hand:

re_sd(fit; scale = :drmtmb, tree = tree)

For a fit whose phylo/relmat term was instead built on the tip-correlation matrix already (fit.phylo_scale === :correlation), scale = :drmtmb is a no-op (that route's raw re_sd already matches drmTMB) and tree is not required. scale = :drmtmb on a fit with no random effects returns the empty Dict. scale = :native (the default) is unchanged from before this option existed.

source

DRModels.vc Function

vc(fit) -> Dict{Symbol,Matrix{Float64}}

Random-effect covariance summary per grouping factor.

• Correlated random-effect block ((1 + x | g)): a 2×2 covariance matrix — sqrt.(diag(vc(fit)[:g])) are the intercept/slope SDs, the off-diagonal their covariance.

• Scalar / structured variance components ((1 | g), relmat/animal/phylo): a 1×1 matrix holding that component's variance σ². A fit with two structured components (e.g. phylo(1|species) + relmat(1|id)) reports both, keyed by grouping factor.

source

StatsAPI.coeftable Function

coeftable(fit::DrmFit; level = 0.95) -> StatsBase.CoefTable

Coefficient table across every parameter block: columns Estimate, Std.Error, z, Pr(>|z|), and a level confidence interval (Lower/Upper). Values are on each block's working scale (μ on the response scale; σ on log σ; ρ12 on atanh ρ12; random-effect SDs on log σ_b). Row names are prefixed with the block (e.g. "mu: (Intercept)") so they stay unique across blocks.

The z and Pr(>|z|) columns are the Wald z = Estimate / Std.Error and its two-sided normal p-value for rows with a finite standard error, except on fixed-prior MAP fits where z, p and interval columns are withheld. The null is always "coefficient = 0 on the working scale of that block":

block	null tested by the z and p columns
:mu, :mu1, :mu2	coefficient = 0
:rho12	atanh ρ12 = 0, i.e. ρ12 = 0
:sigma, :sigma1, :sigma2	coefficient = 0 on log σ: an intercept tests σ = 1 (unit-dependent), a slope tests equal dispersion
:resd, :resid	log σ_b = 0, i.e. σ_b = 1 — not the σ_b = 0 variance boundary
:recov, :phylocov	Cholesky entry = 0 (no single interpretable null)

The table object carries no statement of these nulls (show(fit) prints them under each block heading). To test a variance component against 0 use a boundary-corrected likelihood-ratio test (lrt_boundary), not the :resd row. A boundary / singular direction (non-finite SE) reports NaN z / p rather than a spurious z = 0, p = 1.

When the fit is not converged (is_converged is false) or its Hessian is not positive definite, the Std.Error, z, p and interval columns are all NaN (#1193), as drmTMB refuses vcov() there; show(fit) prints the reason. stderror and vcov still return the Hessian-based values.

source

StatsAPI.coef Function

coef(fit::DrmFit)

Estimated coefficients, all parameter blocks concatenated (coef(fit, :mu) returns one block). Extends StatsAPI.coef.

source

All raw fitted coefficients in the prepared-engine order.

source

coef(fit::JointDrmFit, parameter)

:mi_<variable> maps to the internal predictor-model block. For a Gaussian predictor, :sigma_mi_<variable> returns the natural predictor SD, while :logsd_mi_<variable> returns its raw log-covariance coordinate.

source

All raw fitted coefficients in the two-predictor prepared-engine order.

source

All raw fitted finite-state coefficients in beta, delta, alpha, cutraw order.

source

StatsAPI.vcov Function

vcov(fit::DrmFit)

Variance–covariance matrix of the estimated coefficients. Extends StatsAPI.vcov.

source

StatsAPI.nobs Function

nobs(fit::DrmFit)

Number of observations. Extends StatsAPI.nobs.

source

Fitted scale and correlation

DRModels.sigma Function

sigma(fit)

Fitted scale / dispersion. Returns the per-observation fitted scale(s) computed at the MLE — drmTMB's sigma().

• Univariate location–scale with a single scale (:sigma): the σ_i vector (exp.(Xσ·β̂_σ) on the response scale; for meta-analysis √(vᵢ + τ²)).

• Bivariate co-scale: Dict(:sigma1 => …, :sigma2 => …).

• Families with no separately fitted scale stored (e.g. Poisson, whose variance is fixed by the mean): returns an empty Dict — there is no free dispersion.

For the BetaBinomial sigma-axis random-intercept route, this is the conditional plug-in scale at the fitted log-sigma random-effect modes. ranef(fit) exposes those modes; predict_parameters() reports the population-level sigma predictor and the separate :sd prediction for its random-effect SD.

For the population-level scale coefficients (on the working/log scale) use coef(fit, :sigma) instead.

source

DRModels.corpairs Function

corpairs(fit)

Fitted between-response residual correlation(s) — drmTMB's corpairs(). For a bivariate co-scale model this is the per-observation ρ12 = tanh(Xρ·β̂_ρ) (constant when rho12 ~ 1, varying when rho12 ~ x). Univariate models have no between-response correlation and return an empty Dict.

For random-effect (within-group) correlations, see vc.

source

DRModels.rho12 Function

rho12(fit)

Fitted residual correlation ρ12 for a bivariate model (bf(mu1=…, mu2=…, rho12=…)), on the response scale (ρ12 ∈ (-1, 1)), one value per observation. Mirrors drmTMB's rho12. Errors for univariate fits, which have no residual correlation.

source

DRModels.coevolution_cor Function

coevolution_cor(fit) -> NamedTuple

Among-axis correlation matrix of a q=4 phylogenetic bivariate location–scale ("coevolution") fit: the 4×4 group-level correlation between the four shared-phylogenetic axes (mu1, mu2, sigma1, sigma2),

R = D^{-1/2} Σ_a D^{-1/2},   D = Diagonal(Σ_a),

back-transformed from the log-Cholesky parameterisation of the stored covariance fit.ranef.Sigma_a. The off-diagonals are the coevolutionary correlations — e.g. R[1, 2] = ρ_a(mu1, mu2) is the among-species correlation of the two trait means (the "coevolution of means" signal), R[3, 4] = ρ_a(sigma1, sigma2) the coevolution of the two log-scales, and the mean↔scale entries the lability couplings.

Returns a NamedTuple:

• cor — the 4×4 correlation matrix (symmetric, PD, unit diagonal);

• axes — the axis labels (:mu1, :mu2, :sigma1, :sigma2), the row/column order.

For the residual between-response correlation ρ12 see corpairs; for the raw covariance see coevolution_vc or vc.

Example

fit = drm(bf(mu1 = @formula(y1 ~ x + phylo(1 | species)),
             mu2 = @formula(y2 ~ x + phylo(1 | species)),
             sigma1 = @formula(sigma1 ~ 1 + phylo(1 | species)),
             sigma2 = @formula(sigma2 ~ 1 + phylo(1 | species)),
             rho12 = @formula(rho12 ~ 1)),
          Gaussian(); data, tree = phy)
R = coevolution_cor(fit)
R.cor[1, 2]        # ρ_a(mu1, mu2): coevolution-of-means correlation

source

DRModels.coevolution_vc Function

coevolution_vc(fit) -> NamedTuple

Per-axis phylogenetic variance components of a q=4 coevolution fit: the diagonal of the group-level covariance Σ_a (the among-species variance on each of the four shared-phylogenetic axes) and the matching standard deviations.

Returns a NamedTuple:

• axes — the axis labels (:mu1, :mu2, :sigma1, :sigma2);

• variance — Dict(:mu1 => σ²_a(mu1), …), the per-axis phylo variances (diag(Σ_a));

• sd — Dict(:mu1 => σ_a(mu1), …), their square roots;

• cov — the full 4×4 covariance Σ_a (== fit.ranef.Sigma_a).

The variances are positive by construction (the log-Cholesky parameterisation keeps Σ_a positive-definite). For the correlations see coevolution_cor; vc(fit)[group] returns the same Σ_a keyed by grouping factor.

source

DRModels.coevolution_summary Function

coevolution_summary(fit) -> NamedTuple

Tidy summary of a q=4 coevolution fit's among-axis structure, combining coevolution_vc and coevolution_cor into long-form vectors convenient for printing or assembling a table.

Returns a NamedTuple:

• axes — the axis labels (:mu1, :mu2, :sigma1, :sigma2);

• variance — per-axis phylo variances, in axes order;

• sd — per-axis phylo SDs, in axes order;

• pair — the 6 unordered axis pairs as Tuple{Symbol,Symbol} (upper triangle of the correlation matrix);

• correlation — the among-axis correlation for each pair, in matching order;

• covariance — the among-axis covariance for each pair, in matching order;

• cor — the full 4×4 correlation matrix;

• cov — the full 4×4 covariance matrix Σ_a.

source

Random-effect estimates

DRModels.ranef Function

ranef(fit) -> Dict{Symbol,...}

Per-level conditional random-effect estimates (BLUPs), keyed by grouping factor. These are the posterior means of the random effects at the fitted variance components — drmTMB's ranef().

• Scalar random intercept (1 | g): a Vector of length n_levels(g).

• Correlated (1 + x | g): an n_levels × 2 matrix ([intercept slope]).

• Multiple components (1 | g) + (1 | h): one entry per factor.

Populated for the Gaussian closed-form RE paths (exact GLS conditional means) and structured routes that retain conditional modes. Returns an empty Dict for a fit with no random-effect block at all. Most non-Gaussian Gauss–Hermite/Laplace marginal routes (for example, ordinary Binomial or Poisson (1 | g) fits) integrate random effects without retaining conditional modes, so ranef() throws for those fits. The Tweedie related-matrix mean-intercept SD route is an exception: its sparse Laplace calculation retains the fitted conditional mode.

source

Prediction and simulation

StatsAPI.predict Function

predict(fit::JointFiniteDrmFit, newdata; dpar = :mu, type = :response, se = false)

Known-state plug-in prediction for a formula-fitted ordinal or categorical missing predictor. Mean predictions require every new-data marker value to be an observed fitted level and are evaluated with the retained formula schema; they never condition on a supplied new response. Sigma predictions use their own retained schema and therefore do not require the marker column. Missing new-data marker values and unknown levels remain outside this bounded route.

source

DRModels.predict_parameters Function

predict_parameters(fit, newdata; type = :response, se = false)
    -> Dict{Symbol,Vector{Float64}}  (se = false)
    -> Dict{Symbol,NamedTuple}        (se = true)

Population-level prediction of every distributional parameter at newdata (a NamedTuple / column table), with random / structured effects set to zero (exactly like predict: the average group, not an average over groups). The returned Dict has one entry per distributional parameter the model carries — always :mu and (when the family uses it) :sigma, plus any family extras present (:nu, :zi, :hu, :zoi, :coi). For supported Binomial, BetaBinomial, Poisson, NB2, zero-truncated NB2, Gamma, Beta, Student, or Tweedie random-effect SD fits with sd(group) ~ …, it also contains :sd: the random-effect SD evaluated for each row's group, repeated for rows in the same group.

type = :response (default) applies each parameter's inverse link. Student nu returns the full degrees of freedom 2 + exp(η); SkewNormal nu returns the slant on its identity link; Tweedie nu returns 1 + logistic(η). ZeroOneBeta mu uses the logistic link and returns the interior beta-component mean, rather than the zero/one-inflated mixture mean. Do not assume in-sample equality with marginal_parameters or fitted: stored values can include conditional random effects, a mixture mean, or another scale. This includes LogNormal routes delegated to Gaussian fits that store conditional means, and a meta_V sigma stored as √(v + σ²). type = :link returns the linear predictor Xβ̂ per parameter (the working scale).

For a univariate fit the parameters are :mu, (:sigma) plus family extras; for a bivariate fit they are :mu1, :mu2, :sigma1, :sigma2, :rho12 (each from its own fixed-effects RHS, with the σ links exp and the ρ12 link tanh).

se = false (default) returns Dict{Symbol,Vector{Float64}} of point values. se = true returns Dict{Symbol,NamedTuple} with p => (; value, se) per parameter: each se is the delta-method standard error using that parameter's own coef range r_p from fit.blocks (V_p = vcov(fit)[r_p, r_p]), the response scale multiplying by that parameter's inverse-link derivative at η̂ (:sigma→exp, :rho12→1−ρ², etc.). value matches the se = false point prediction.

Example

x = randn(200)
y = 0.5 .- 0.8 .* x .+ exp.(-0.3 .+ 0.4 .* x) .* randn(200)
data = (; y, x)
fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data)

p = predict_parameters(fit, data)              # Dict(:mu => …, :sigma => …)
p[:mu]    ≈ fit.means[:mu]                      # in-sample reproduction
p[:sigma] ≈ fit.scales[:sigma]

predict_parameters(fit, data; type = :link)[:mu]   # == Xβ̂ (predict link scale)

source

predict_parameters(fit::DrmFit{Tweedie}, newdata;
                   type = :response, se = false, K = nothing)

Predict fitted distributional parameters for new observations. For a related-matrix SD target, :sd is the multiplier of the supplied kernel. Pass K to also return :marginal_sd, using the covariance matrix ordered by the first-seen group labels in newdata. With se = true, both entries include delta-method standard errors. Without K, only the SD multiplier is returned.

source

DRModels.marginal_parameters Function

marginal_parameters(fit) -> Dict{Symbol,Vector{Float64}}

In-sample fitted per-observation distributional parameters, read straight from the stored fit — a cheap accessor with no recomputation. Returns the mean(s) from fit.means (:mu, or :mu1/:mu2 for a bivariate fit) and every per-observation scale / correlation parameter from fit.scales (e.g. :sigma, :nu, :zi, :hu, :zoi, :coi; :sigma1/:sigma2/:rho12 for a bivariate fit).

Despite the name, this accessor does not integrate over random effects or recompute predictions. Stored values depend on the fitted route: some include conditional random effects or use a mixture mean or another scale (see fitted). In particular, Gaussian sparse phylogenetic routes and LogNormal routes delegated to them can store conditional means. Use predict_parameters for fixed-effect predictions with random and structured effects set to zero; do not assume in-sample equality on every route. For BetaBinomial fits with a random effect on sigma, :sigma here uses the conditional log-sigma modes, while predict_parameters() reports the population-level fixed-effect sigma predictor.

Example

fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data)
m = marginal_parameters(fit)
m[:mu]    == fit.means[:mu]
m[:sigma] == fit.scales[:sigma]

source

DRModels.prediction_grid Function

prediction_grid(reference::NamedTuple; n::Int = 50, kwargs...) -> NamedTuple

Build a newdata column table for predict / predict_parameters by sweeping one or more predictors over supplied value ranges (their Cartesian product) while holding every other predictor fixed at a reference value.

Pure data — no fitted model is needed, so it is trivially testable and composes directly with predict_parameters(fit, prediction_grid(...)).

Arguments

• reference::NamedTuple: the predictor columns to hold constant. Each held value is reduced to a scalar by this rule:

  • if reference[col] is an AbstractArray of numbers → its mean;

  • if reference[col] is any other AbstractArray → its first element;

  • otherwise → the value itself (already a scalar).

  Held columns are broadcast to the product length.

• n::Int = 50: reserved for future default-range generation; currently unused (every swept predictor supplies its own explicit values via kwargs).

• kwargs...: each predictor = values gives a vector/range of values to sweep for that predictor. The output rows enumerate the full Cartesian product of the swept predictors. A swept predictor named in reference overrides (replaces) the held value.

Returns

A NamedTuple of equal-length column vectors. With no swept predictors the grid has a single row at the reference; with one swept predictor it is just that range (others held); with several, the full Cartesian product.

The swept columns appear first (in kwargs order), then the remaining held columns (in reference order).

Example

g = prediction_grid((; x = randn(100)), x = range(-2, 2; length = 25))
length(g.x) == 25                       # a 25-row sweep over x

# Two swept predictors → Cartesian product (5 × 3 = 15 rows), z held at its mean:
g2 = prediction_grid((; x = [0.0], z = [1.0, 2.0, 3.0]), x = -2:1.0:2)
length(g2.x) == 5
all(==(2.0), g2.z)                       # z held at mean([1,2,3])

# Composes with a fit:
preds = predict_parameters(fit, g)       # Dict(:mu => …, :sigma => …) of length 25

source

DRModels.simulate Function

simulate(fit; nsim = 1, rng = default_rng(), newdata = nothing,
         re = :marginal, data = nothing, K = nothing, A = nothing,
         tree = nothing, coords = nothing)

Draw parametric replicate response(s) from the fitted model, the building block of a parametric bootstrap and posterior-predictive checks. Each draw uses the fitted fixed-effect mean and the fitted dispersion / scale parameters for the family.

Random effects (re, #1092):

• re = :marginal (the default) draws FRESH random effects for every replicate from their fitted distribution, u* ~ N(0, sd^2 K), adds them to the fixed-effect linear predictor and then draws the response. This is the same marginal simulator the parametric bootstrap uses (#459, #462), and it matches drmTMB's default simulate(fit) (re.form = NA). A random-effect fit needs data =, the table it was fitted to (and tree / K / A / coords for a structured term, as in bootstrap_ci), because the fit does not keep the grouping columns. Without data, or for a random-effect structure the marginal simulator does not cover, simulate throws an ArgumentError rather than fall back to b = 0.

• re = :zero draws with every random effect set to zero (population level), which was the only behaviour before #1092. It needs no data. drmTMB has no such mode; its re.form = NULL (conditional on the fitted modes) has no counterpart here. Do not use re = :zero draws for a parametric bootstrap or a posterior-predictive check of random-effect variation: they carry no between-group signal.

A fit without random effects ignores re and data. A temporal() fit always draws the full marginal model (a fresh chain per series), as before.

Return value (univariate / random-effect / meta models):

• nsim == 1 → a length-nobs response Vector (back-compatible).

• nsim > 1 → a nobs × nsim Matrix, one independent replicate per column.

Bivariate Gaussian models return Dict(:mu1=>…, :mu2=>…) for nsim == 1, or a length-nsim Vector of such Dicts for nsim > 1 (a matrix of paired responses is not well defined).

Supported families: Gaussian (univariate & bivariate), Student-t, SkewNormal, Poisson (+ zero-inflated / hurdle), NegBinomial2 (+ zero-inflated / hurdle / truncated), Beta, BetaBinomial, Binomial, Gamma, LogNormal, ZeroOneBeta, Tweedie, and CumulativeLogit.

Example

fit = drm(bf(@formula(y ~ x), @formula(sigma ~ x)), Gaussian(); data)
y1  = simulate(fit)               # Vector, length nobs (a homtoep fit that dropped
                                  # missing responses: one per data row, NaN there)
Y   = simulate(fit; nsim = 100)   # 100 columns, one per draw
ys  = simulate(fit; newdata = (x = [0.0, 1.0],))  # one row per newdata row

newdata (optional) is a NamedTuple / column table. Draws are population-level (random / structured effects set to zero, whatever re says), one row per newdata row, using predict_parameters for the distributional parameters. Omit it to draw for the training design. Temporal fits refuse newdata.

source

StatsAPI.fitted Function

fitted(fit)

Fitted value(s) of the mean parameter: the μ vector for univariate models, Dict(:mu1=>…, :mu2=>…) for bivariate models. On a random-effect fit this is the population-level value with the random effects at zero, not the marginal mean averaged over them. For the GLM-type families (Gaussian, Poisson, NegBinomial2, Gamma, Beta, Binomial) it is g⁻¹(Xβ̂), the value predict gives in-sample. Which mean is stored depends on the route:

• Population mean Xβ̂: fixed-effects fits, ordinary (1 | g) and crossed random intercepts, relmat / animal, and the dense algorithm = :gls / :lbfgs Gaussian phylo(1 | g) fitters.

• Conditional mean Xβ̂ + Sû (the BLUPs included, as drmTMB's fitted-row predict()): the sparse L-BFGS (the default) and EM Gaussian phylo(1 | g) routes, and so residuals(fit) on them are conditional residuals. predict(fit, data) is the fixed-effect Xβ̂ on every route.

Other exceptions:

• ZeroOneBeta stores the zero/one-inflated mixture mean and CumulativeLogit the expected category score;

• LogNormal fits report the response-scale median exp(·) of the corresponding log-scale mean above, including the phylo / relmat and crossed routes (#1102), so on its sparse phylo(1 | g) route the median includes the fitted phylogenetic effects.

drmTMB's fitted() adds the fitted random-effect modes and, for some families (LogNormal, zero-inflated, hurdle and truncated counts), returns the unconditional response mean, so the two packages' fitted values are not directly comparable.

source

StatsAPI.residuals Function

residuals(fit; type = :response, rng = Random.default_rng())

Model residuals. type selects the kind:

• :response (default) — raw response residuals (observed − fitted mean), matching fitted's shape. residuals(fit) is unchanged.

• :quantile — standardised residuals. Outside temporal() fits these are randomized quantile residuals (Dunn & Smyth; DHARMa / glmmTMB style): for observation i with reference distribution F_i, r_i = Φ⁻¹(u_i), where u_i is the (randomized, for discrete families) probability-integral transform of y_i. Univariate only.

What a correct model implies depends on the route:

• no random effect: F_i is the fitted distribution, and the r_i are approximately i.i.d. standard normal (zero-inflated and hurdle count fits included: their PIT uses the whole mixture CDF, below);

• one random intercept: F_i integrates it out over its fitted SD (32-node Gauss–Hermite), so each r_i is approximately standard normal, but residuals of one group stay correlated. On the mean, this covers any lone mean random intercept, including a phylo() or relmat() term; the fitted SD σ̂_b is used as every observation's random-effect SD, which matches the model when the relatedness matrix has a unit diagonal (for example an ultrametric tree of height 1). The 32-node rule is accurate to about 1e-6 in r_i while σ̂_b / σ̂ ≤ 2, but its error grows beyond that (6e-3 at 4, 0.12 at 10). The default sparse / EM Gaussian phylo(1 | g) route is exact instead: F_i is the closed-form marginal N(x_i'β̂, σ̂² + σ̂_b² C_ii), centred on the fixed-effect mean, with the tip's own variance C_ii from the tree (#995, #996). On the scale, for a Gaussian sigma ~ 1 + (1 | g), F_i mixes N(μ_i, (σ_i e^b)²) over b ~ N(0, τ̂²), with σ_i the fixed-effect scale. When τ̂ is moderate or large these residuals can be under-dispersed (SD below 1) even under a correct model, because the default marginal = :LA fit integrates the log-σ intercept on a fixed prior-scale grid that loses accuracy as the SD grows; refitting with marginal = :AGHQ is the check. drmTMB's quantile residuals for these fits condition on the fitted modes instead;

• other random-effect shapes (crossed, correlated slopes, a mean and a sigma intercept together) and CumulativeLogit fits: judged with the random effects set to 0. The group-level variation then stays in the residuals, so under a correct model they are over-dispersed (SD above 1) and a QQ plot that bends away from the line is not by itself evidence of misfit;

• temporal() AR1 / OU, and the paired phylo() + OU fit: residuals conditional on the fitted modes (below). They are not PIT residuals, and their variance is below 1 even when the model is true;

• temporal() homogeneous Toeplitz: whitened residuals (below), approximately i.i.d. standard normal.

The far tail differs too. The PIT routes (every fit without temporal(), Gaussian included) clamp u_i to [eps, 1 − eps], so |r_i| ≤ 8.126 however extreme y_i is: a Gaussian residual of 13.6 σ̂ is reported as 8.126. The temporal() routes do not clamp.

Quantile residuals are implemented for every DRModels.jl response family except Tweedie (no closed-form CDF in Distributions.jl) and SkewNormal:

• continuous (PIT u_i = F(y_i), no RNG): Gaussian, Student-t, LogNormal, Gamma, Beta;

• discrete, randomized (u_i = F(y_i−1) + (F(y_i) − F(y_i−1))·U, U ~ Uniform(0,1) drawn from rng): Poisson, TruncatedPoisson, NegBinomial2, TruncatedNegBinomial2, Binomial, BetaBinomial, CumulativeLogit (ordinal);

• atomic (point-mass mixture; the mass is randomized across): ZeroOneBeta.

Zero-inflated and hurdle Poisson / NegBinomial2 fits (zi ~ …, hu ~ …) are randomized within the jumps of the whole mixture CDF, zero part included: F(y) = π + (1 − π)·F_c(y) with zero-inflation probability π, and F(0) = p₀, F(y) = p₀ + (1 − p₀)·F_t(y) for y ≥ 1 with hurdle probability p₀, where F_c is the count CDF at the fitted count mean and F_t its zero-truncated CDF.

The per-family parameter → distribution map lives in _conditional_dist (reused by future simulate/PIT checks). Tweedie and SkewNormal throw an ArgumentError.

temporal() fits (as drmTMB's Pearson residuals; for AR1 and OU, drmTMB's quantile residuals are the same numbers):

• AR1 and OU, with or without (1 | id) or a paired phylo() term, are conditional on the fitted modes: (y_i − x_iᵀβ̂ − b̂_id − ŝ_i) / σ̂, where ŝ_i is the temporal mode (ranef(fit)[id]) and b̂_id the (1 | id) or phylogenetic mode. This is the estimated residual noise E[ε_i | y] / σ̂, after the fitted temporal path: it is not whitened against the marginal covariance, and its variance is below 1 even when the model is true. At a residual-SD boundary (sigma_ratio among check_drm(fit).temporal_boundary.findings) the temporal path absorbs the data: σ̂ → 0 drives these residuals toward 0, so they cannot reveal outliers there. Check check_drm(fit).temporal_boundary before reading them.

• Homogeneous Toeplitz is the whitened L⁻¹(y − Xβ̂), with L the Cholesky factor of each series' σ̂²R̂: there are no modes to condition on. drmTMB's type = "quantile" residuals for this structure are not whitened ((y − Xβ̂)/σ̂), so compare with its type = "pearson".

type = :response stays y − Xβ̂ (population level, matching fitted) for these fits; drmTMB's residuals(fit, type = "response") subtracts its conditional fitted() instead.

source

Predicting distributional parameters

Prediction interpretation

predict_parameters sets random and structured effects to zero; it does not integrate over their fitted distributions. With a nonlinear link, these are different predictions: for example, exp(η) differs from averaging exp(η + b) over a non-degenerate random effect b. Interpret the result as a fixed-effect prediction.

predict with type = :response returns invlink(Xβ) for the mean or location parameter, with random and structured effects set to zero. Its interpretation depends on the family: LogNormal returns a response-scale median; truncated Poisson and NB2 return the untruncated base rate or mean, rather than the expected positive count. ZeroOneBeta returns the interior beta-component mean, and CumulativeLogit returns the location predictor rather than an expected category score. A distributional regression also models scale and, for bivariate responses, correlation. Use predict_parameters to obtain the population-level value of every distributional parameter the model carries (:mu, :sigma, plus any family extras) at new covariate values, also with random / structured effects set to zero. These values apply each parameter's inverse link to its fixed-effect linear predictor; they do not integrate over the random-effect distribution. marginal_parameters is the cheap in-sample accessor that reads the fitted per-observation parameters straight off the fit. prediction_grid builds the new-data table to sweep over (varying chosen predictors, holding the rest at a reference value).

using DRModels, Random
Random.seed!(20260603)

x = randn(500)
y = 0.5 .- 0.8 .* x .+ exp.(-0.3 .+ 0.4 .* x) .* randn(500)
data = (; y, x)

# Gaussian location–scale fit: both μ and σ depend on x.
fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data = data)

# Per-distributional-parameter prediction over a covariate sweep.
grid = prediction_grid((; x = data.x), x = range(-2, 2; length = 11))
p = predict_parameters(fit, grid)   # Dict(:mu => …, :sigma => …) on the response scale
p[:mu]                              # predicted mean across the sweep
p[:sigma]                          # predicted scale across the sweep

# Working (link) scale instead: returns Xβ̂ per parameter.
predict_parameters(fit, grid; type = :link)[:mu]

# In-sample fitted parameters, straight off the fit (no recomputation).
marginal_parameters(fit)            # == predict_parameters(fit, data) in-sample

Formula-fitted finite-state missing predictors

Experimental

Exported for evaluation, not yet stable. API and numerics may change. The R bridge accepts the bounded inputs shown in the joint missing-predictor guide.

For one ordinal or categorical missing predictor in the bounded Gaussian joint route, construct the predictor model with impute_model; use CategoricalLogit() for nominal states. JointFiniteDrmFit retains the raw kernel coefficients and covariance. For ordinal predictors, cutpoints(fit) provides the constrained cutpoints separately. Modelled missing predictors shows the public formula route, predictions, and current limits.

DRModels.CategoricalLogit Type

Predictor-only multinomial-logit marker for finite categorical imputation.

source

DRModels.JointFiniteDrmFit Type

JointFiniteDrmFit

Formula-facing wrapper for one ordinal or categorical missing predictor. Its raw coefficient vector and covariance keep the kernel order beta, delta, alpha, cutraw; ordinal raw cutpoints are deliberately separate from the predictor-logit coefficients and cutpoints(fit) returns their natural scale. Direct formula fits retain applied mean and sigma schemas for known-state new-data prediction; the three-argument constructor remains available for manually wrapped prepared fits but those fits cannot predict new data.

source

DRModels.cutpoints Function

cutpoints(fit::JointFiniteDrmFit)

Natural cumulative-logit cutpoints for an ordinal finite predictor. This is a separate transformed view of raw coef(fit, :rawcut_<variable>); vcov(fit) continues to use the raw kernel coordinates.

source

Prepared-model fitting

The low-level Experimental fit_prepared_joint overloads accept PreparedJointModel, PreparedFiniteJointModel, and PreparedTwoJointGaussianModel. Each uses initial = prepared_joint_initial(model) and g_tol = 1e-8 by default. The initial vector must have the model's parameter length and finite entries; g_tol must be finite and positive. These overloads use L-BFGS with forward automatic differentiation. They return PreparedJointFit, PreparedFiniteJointFit, and PreparedTwoJointGaussianFit, respectively, with convergence and covariance status retained in their metadata. See the fit_prepared_joint for the lower-level interface; these overloads do not expand the formula or R-bridge admission boundaries above.

DRModels.PreparedJointModel Type

PreparedJointModel

Validated, formula-free input for a Gaussian response with one potentially missing Gaussian or Bernoulli predictor. Xmu excludes that predictor; its coefficient is the separate b coordinate in the fixed parameter order beta, b, delta, alpha, kappa (with kappa only for a Gaussian predictor).

source

DRModels.PreparedJointFit Type

PreparedJointFit

Result wrapper for the prepared prototype. fit is a normal DrmFit; the separate metadata field prevents missing-data information from being confused with random effects. Generic confint, simulate, and AIC forwarding is intentionally absent.

source

DRModels.PreparedFiniteJointModel Type

Validated Gaussian-response model with one ordinal or categorical predictor.

source

DRModels.PreparedFiniteJointFit Type

Fitted exact prepared ordinal or categorical missing-predictor model.

source

DRModels.prepared_joint_model Function

prepared_joint_model(y, x, Xmu, Xsigma, Xpredictor;
                     predictor, mu_names, sigma_names, predictor_names,
                     original_row)

Construct a complete-design, original-row-preserving prepared model. Missingness is permitted only in y and x; masks are data properties, never parameters.

source

prepared_joint_model(y, x::AbstractMatrix, Xmu, Xsigma, Xpredictor;
                     predictor = :gaussian, predictor_variables = (:x1, :x2), ...)

Construct the exact prepared model for two conditionally independent Gaussian predictors. Xmu excludes both marked predictor columns. The fixed parameter order is beta, b1, b2, delta, alpha1, logtau1, alpha2, logtau2.

source

prepared_joint_model(y, x, X_mu_state, Xsigma, Xpredictor;
                     predictor, levels, variable, ...)

Build the prepared finite-state Gaussian-response joint model. X_mu_state contains the complete mean design for every observation-state combination in state order; it must not be reconstructed from an expected score or a modal state. The prepared raw parameter order is beta, delta, alpha, cutraw. For ordinal predictors alpha has q entries and cutraw has K - 1; for categorical predictors alpha is level-major over the K - 1 non-baseline logits and no cutpoint block is present.

source

DRModels.prepared_joint_rowloglik Function

prepared_joint_rowloglik(model, theta)

Return exact marginal log-likelihood contributions in original prepared-row order. Missing x is integrated analytically for a Gaussian predictor and by finite logit-stable enumeration for a Bernoulli predictor.

source

Exact marginal row log likelihood in prepared-row order for the two-Gaussian route.

source

Exact prepared-row likelihood for an ordinal or categorical missing predictor.

source

DRModels.prepared_joint_conditional_moments Function

prepared_joint_conditional_moments(model, theta)

Return named vectors (mean, variance, status) for x conditional on the available row data. The returned values retain the original prepared-row order.

source

prepared_joint_conditional_moments(model::PreparedTwoJointGaussianModel, theta)

Return marginal means, full 2-by-2 row posterior covariance matrices, and per-predictor statuses. Off-diagonal covariance is retained when both predictors are missing and the response is observed.

source

prepared_joint_conditional_moments(model, theta)

Return posterior state probabilities in original prepared-row order. Ordinal mean and variance are the conditional score mean and variance. Categorical mean is the first modal state code and its variance is unavailable (NaN).

source

DRModels.fit_prepared_joint Function

fit_prepared_joint(model; initial = prepared_joint_initial(model), g_tol = 1e-8)

Fit the exact prepared prototype by forward-mode AD and LBFGS. Optimizer and covariance status are recorded separately; no failed or non-finite calculation is converted into a finite likelihood.

source

Fit the finite-state joint likelihood by forward-mode AD and LBFGS.

source

Prepared negative log-likelihood

prepared_joint_nll(model, theta) returns -sum(prepared_joint_rowloglik(model, theta)) for the prepared Gaussian, finite-state, and two-Gaussian-predictor routes. It evaluates the supplied parameter vector without fitting the model. Predictor-only rows retain their predictor likelihood; rows with the response and all modelled predictors missing contribute zero. These are Experimental array interfaces with the same admission and uncertainty limits described above.

Prepared-model summaries and finite-state imputations

joint_missing_summary(fit) returns a NamedTuple whose fields depend on the prepared route. Finite-state fits return predictor, variable, levels, original_row, observed_y, observed_x, conditional_probabilities, conditional_mean, conditional_variance, conditional_status, all_rows, predictor_only_rows, uncertainty_status, optimizer_status, and covariance_status. A JointFiniteDrmFit forwards to this same summary. Two-predictor Gaussian fits return predictor_variables, original_row, observed_y, observed_x (a pair of masks), conditional_mean, conditional_covariance, conditional_status, all_rows, predictor_only_rows, uncertainty_status, optimizer_status, and covariance_status. Returned arrays are copies of the retained metadata.

For a missing finite-state predictor, imputed(fit) reports an ordinal conditional expected score or a categorical modal state index in estimate. The states are numbered 1:length(levels); levels supplies their labels. The source column distinguishes conditional_expected_score from conditional_modal_category. With se = true and usable retained covariance, an ordinal std_error is the conditional state-score SD; it does not include parameter-estimation uncertainty. Categorical std_error is unavailable. Check uncertainty_status, and use the summary's conditional_probabilities for state probabilities. Observed rows keep their observed state index when rows = :all; the default rows = :missing selects missing predictor rows.

Inference

Check the status as well as the bounds of a profile interval. A signed infinite bound can mean that no crossing was found within the searched range, or that an endpoint solve failed. profile_result distinguishes these outcomes in stats and failed; confint warns about failed endpoint solves.

For coupled non-Gaussian location–scale fits, endpoint_diagnostics also records why each endpoint search stopped, its last evaluated candidate and its residual. A candidate from a failed search is diagnostic information, not a confidence limit. Through the experimental R bridge (engine = "julia"), inspect conf.status and profile.message in the returned interval table. A failed result must not be read as evidence of an unbounded interval.

StatsAPI.confint Function

confint(fit; level = 0.95, method = :wald, threads = false, parm = nothing)

Confidence intervals for every coefficient, as a vector of (param, coef, estimate, lower, upper) rows on each parameter's working scale (μ on the response scale; σ on log σ; ρ12 on atanh ρ12; random-effect SDs on log σ_b).

• method = :wald (default) — estimate ± z·se from the stored covariance.

• method = :profile — profile-likelihood interval: the endpoints where 2(ℓ̂ − ℓ_profile) = χ²₁(level), re-optimising the nuisance parameters at each fixed value (asymmetric, and exact under the LR statistic where Wald is only quadratic-approximate). Works on any fitted Gaussian model, and on the non-Gaussian canonical location–scale fit ((1 | tag | group) coupled RE), where it routes to a constrained L-BFGS nuisance solve with fresh objective and gradient acceptance checks on the variance boundary. The endpoint search uses warm-start continuation (each profiled solve starts from the previous point's optimum) and a guarded-Newton root-find driven by the envelope-theorem slope ∂nll/∂θ_k, falling back to bisection. This guarded local search records explicit failed endpoint arms when it cannot certify an evaluated root. Endpoint validity assumes an accurate inner nuisance solve: the :finite autodiff path can leave the profiled NLL slightly non-monotone. When that is detected the bracket is reset to [0, t] and pure bisection is used. This is a guarded local search, not a guarantee of the globally first LR crossing; a nonmonotone flag is set on the profile stats row so callers can inspect that limitation. An endpoint arm that the search cannot certify is REFUSED, not returned: this method throws an ArgumentError naming the coefficient, the arm, and the nuisance-solve reason rather than reporting the failed side as a signed Inf. Use profile_result when you want the same rows plus the per-endpoint diagnostics instead of an exception. Pass threads = true to profile coefficients in parallel when the fitted objective is thread-safe; if only one coefficient is profiled, its lower and upper endpoint searches are run in parallel instead. Coupled location–scale random-effect fits (a shared (1 | p | group) term on mu and sigma) also profile coefficients in parallel, but each coefficient's lower and upper endpoint searches stay serial, so profiling a single coefficient there gains nothing from threads = true. A fit that explicitly selected DRMControl(optimizer = :lbfgs | :bfgs) is refused before profiling, because nuisance refits cannot yet forward it.

• parm = :resd or parm = [:mu, :resd] — restrict intervals to one or more parameter blocks. This is especially useful for profiling random-effect SDs without also profiling fixed-effect coefficients. A random-effect SD (:resd) whose lower arm reaches the numerical floor without an LR crossing is reported on the SD scale as [0, u] (drmTMB profile.boundary), not as a thrown error and not as (-Inf, Inf).

Mirrors drmTMB's confint(fit, method = "wald" | "profile").

source

StatsAPI.stderror Function

stderror(fit) -> Vector{Float64}

Wald standard errors (√diag of the covariance), in the fit's coefficient order.

A coefficient's Wald SE is defined only where its estimated variance is finite and positive. At a singular boundary the observed information is not positive-definite, so the stored covariance carries a non-positive (or non-finite) variance for the unidentified direction. Rather than return a silent NaN there, that coefficient reports Inf — an undefined, infinitely wide standard error — which propagates to an unbounded (-Inf, Inf) Wald interval. check_drm flags the same situation via vcov_posdef; drmTMB returns all-NaN from sdreport in this case.

source

DRModels.profile_result Function

profile_result(fit; level = 0.95, threads = false, parm = nothing)

Auditable profile-likelihood confidence intervals. Returns a NamedTuple with:

• ci — the same rows confint(fit; method = :profile) returns, except that a FAILED endpoint arm is kept here as a signed Inf alongside its lower_endpoint_failed / upper_endpoint_failed flag. This is the auditable surface; confint refuses such a row rather than returning it;

• stats — per-coefficient endpoint work counts;

• endpoint_diagnostics — canonical location–scale endpoint reason, last evaluated candidate, and residual for each arm; other profile backends omit this additive diagnostic;

• attempted, used, failed — coefficient counts;

• threaded, worker_threads, julia_threads, blas_threads, blas_oversubscribed, elapsed — CPU context and wall-clock timing;

• autodiff — profile nuisance-gradient backend (:stored, :forward, :finite, or canonical location–scale :locscale).

For generic thread-safe objectives, set threads = true to parallelise independent profile coefficients; when one coefficient is profiled, its endpoint arms may run in parallel. Canonical location–scale profile jobs use only the coefficient-level policy: each job owns its nuisance state, while its lower and upper endpoint chains remain serial. A fit with an explicit DRMControl.optimizer is refused before profiling. Threaded profiling is refused for structured Binomial random-effect SD fits because their Laplace objectives use mutable shared warm-start state; ordinary independent Binomial quadrature routes are not subject to that restriction.

source

DRModels.bootstrap_ci Function

bootstrap_ci(formula, family::Gaussian; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8)
bootstrap_ci(formula, family; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true)
bootstrap_ci(fit::DrmFit{<:Gaussian}; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8)
bootstrap_ci(fit::DrmFit; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true)

Parametric bootstrap confidence intervals: fit the model, then simulate B replicate responses, refit each, and take percentile intervals per coefficient. Univariate-response models (fixed / random-effect / meta / structured). Same row shape as confint. Set threads = true to refit bootstrap replicates in parallel. Pass through the structured-provider keywords (K / A / tree / tree1 / tree2 / coords) and, for Gaussian fits, the solver controls (algorithm / g_tol) exactly as to drm. Use bootstrap_result when you need attempted/used/failed counts and per-replicate failure messages. A thrown refit is recorded; the default failures = :error then aborts. Pass failures = :skip to keep the successful draws. If you already have fit = drm(...), pass the fit directly to avoid refitting the base model before the bootstrap replicates. A homogeneous Toeplitz temporal() fit returns the mean-coefficient rows only (drmTMB refuses the temporal bootstrap and defers σ and lag-correlation intervals). tree1 and tree2 are used by the direct Julia Binomial phylo_interaction() SD route. The fit-based form refuses a fit with an explicit DRMControl.optimizer before any bootstrap refit; the formula-based form keeps its default solver.

The intervals are plain percentile intervals. For a sigma coefficient in a mixed model (a model with random effects) they are too narrow, and bootstrap_ci warns and recommends confint(fit, method = :profile) (percentile covered 83 to 90% at a nominal 95% in a simulated Gaussian location-scale model with a mean random intercept; profile covered 92 to 94%).

source

DRModels.bootstrap_summary Function

bootstrap_summary(formula, family::Gaussian; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8)
bootstrap_summary(formula, family; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true)
bootstrap_summary(fit::DrmFit{<:Gaussian}; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8)
bootstrap_summary(fit::DrmFit; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true)

Parametric bootstrap coefficient summaries in one pass: point estimate, bootstrap standard error, and percentile confidence interval. This is the bootstrap analogue of using stderror(fit) plus confint(fit), avoiding a second bootstrap run when both SEs and intervals are needed. Row fields are (param, coef, estimate, std_error, lower, upper). By default, a thrown refit is recorded and then aborts (failures = :error). Set failures = :skip to compute summaries from successful replicates; call bootstrap_result to inspect the skipped failures. A homogeneous Toeplitz temporal() fit returns the mean-coefficient rows only, as bootstrap_ci. For the direct Julia Binomial phylo_interaction() SD route, supply tree1 and tree2 to preserve the fitted pair field during refits. The fit-based form refuses a fit with an explicit DRMControl.optimizer before any bootstrap refit.

source

DRModels.bootstrap_result Function

bootstrap_result(formula, family::Gaussian; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8, method = :ML, marginal = :LA, penalty = nothing)
bootstrap_result(formula, family; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, method = :ML, marginal = :LA, penalty = nothing)
bootstrap_result(fit::DrmFit{<:Gaussian}; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true, algorithm = :auto, g_tol = 1e-8)
bootstrap_result(fit; data, B = 300, level = 0.95, rng = default_rng(), K = nothing, A = nothing, tree = nothing, coords = nothing, threads = false, failures = :error, check_converged = true)

Auditable parametric bootstrap. Returns a NamedTuple with:

• summary: the same rows returned by bootstrap_summary;

• failures: rows (replicate, seed, message) for failed refits;

• attempted, used, failed: replicate counts;

• status: "bootstrap" when every replicate was kept, "bootstrap_incomplete" when any replicate was dropped and at least two remain, "bootstrap_at_boundary" when those retained draws pile up on a variance-component or correlation bound, or "bootstrap_unavailable" when exactly one replicate remains. If none remain, the call raises an error;

• boundary_params: "param:coef" names that triggered the boundary status;

• seeds: the per-replicate seeds used for reproducibility;

• estimation_method, marginal: the estimator of the seed fit (:ML, :REML or :MAP; :LA, :Laplace, :VA or :AGHQ);

• refit_estimation_method, refit_marginal: the estimator the retained replicates used (:mixed if their estimators differ). A warning marks a mismatch with a known seed estimator;

• threaded: whether threaded refits were actually used.

• worker_threads, julia_threads, blas_threads, blas_oversubscribed, elapsed: CPU context and wall-clock time for the simulated-refit phase.

failures = :error (default) records a thrown refit and then aborts. Pass failures = :skip to keep the successful draws and return status = "bootstrap_incomplete" with a warning. A refit that does not converge, or that lands on a degenerate optimum, is not that kind of error: with check_converged = true (the default) it is dropped, status is "bootstrap_incomplete" (or "bootstrap_at_boundary" when the retained draws also sit on a bound), and a warning names that status. A refit at a variance boundary is kept: convergence_status calls it :boundary, but it is a real draw. Set check_converged = false only to reproduce the old behaviour that counted unconverged refits in the interval. Passing an existing DrmFit reuses that point estimate as the bootstrap seed fit and starts directly with the B simulated refits. g_tol controls the optimizer tolerance for the seed fit and refits when a formula is supplied, and for refits when an existing fit is supplied. The q4 Gaussian phylogenetic shortcut delegates to bootstrap_sigma_a and does not use this g_tol argument; call that function directly with q4_g_tol to change its optimizer tolerance. Other Gaussian bootstrap refits also pass algorithm through to drm(...); this is useful for large structured models where :auto selects a sparse route. For a homogeneous Toeplitz temporal() fit, summary holds the mean-coefficient rows only, as bootstrap_ci. For the direct Julia Binomial phylo_interaction() SD route, pass the same tree1 and tree2 to preserve the fitted product field in simulation and refits. The fit-based form refuses a fit with an explicit DRMControl.optimizer before any bootstrap refit. The fit-based form repeats the seed fit's REML, marginal (:Laplace, :VA, :AGHQ) and MAP penalty in every replicate. The formula-based form takes method (:ML or :REML), marginal and penalty and uses them for the seed fit and the replicates alike (DRModels.jl#1025). The selected family route must support any requested non-default estimator or penalty. algorithm and g_tol are Gaussian-overload controls. A joint mi(x) fit (impute/missing controls) is not bootstrapped: there is no bootstrap_result method for JointDrmFit, and the call raises a MethodError rather than refitting with default settings.

Silent conditional fallback (#997)

When the marginal simulator cannot be built, bootstrap_result draws replicates from the CONDITIONAL simulate(fit) (random effects held at b = 0), which collapses variance-component intervals (#459), and it does so without a warning. Known cases, kept as @test_broken in test/test_marginal_simulator_fallback_997.jl: crossed ordinary bars (1 | g1) + (1 | g2) on the mean, a mean grouping column missing from data, and a correlated slope (1 + x | g). Only a meta_V(...) + random-effect fit throws. Pass the same data (and tree / K / A) used to fit, and treat variance-component intervals from these fits as unreliable.

source

Information criteria

DRModels.loglik Function

loglik(fit) -> Float64

Maximised log-likelihood of the fitted model.

For a REML fit (drm(...; method = :REML)) this returns the restricted log-likelihood (reml_loglik(fit)). REML log-likelihoods are not comparable across different fixed-effect (mean) structures — the error-contrast basis differs — so do not use them for model selection across mean structures; the aic/bic/lrtest guard enforces this. Use ml_loglik (the plain ML log-likelihood at the REML estimate) when an ML-comparable value is needed.

source

StatsAPI.loglikelihood Function

loglikelihood(fit) -> Float64

The StatsAPI/StatsBase-facing alias for loglik — anything that dispatches on the generic loglikelihood (from the StatsAPI.StatisticalModel interface) can call it on a DrmFit. Returns exactly the same value as loglik(fit), with the same REML caveat.

source

DRModels.ml_loglik Function

ml_loglik(fit) -> Float64

The plain (unrestricted) maximum-likelihood log-likelihood. For an ML fit this equals loglik; for a REML fit it is the ML log-likelihood evaluated at the REML parameter estimate — the value to use when an ML-comparable log-likelihood is needed (e.g. across different mean structures).

source

DRModels.reml_loglik Function

reml_loglik(fit) -> Float64

The restricted (REML) log-likelihood. Returns NaN for an ML fit (REML was not used). See loglik for the cross-structure-comparison caveat.

Normalisation convention

For the univariate fixed-effect Gaussian location–scale REML and the Gaussian mean (1 | g) REML, this value is the normalised Patterson– Thompson restricted log-likelihood — the same convention lme4, glmmTMB and TMB report, so it is directly comparable to logLik() from those packages.

The bivariate q=2 and q=4 Laplace REML routes (src/reml_q2.jl, src/reml_q4.jl — reached via structured/phylo bivariate fits with method = :REML) report the same normalised scale.

The normalised form includes (n_β/2)·log(2π), where n_β counts the marginalised fixed effects. For a q=4 phylogenetic model with n_β = 6, this constant is 3·log(2π) ≈ 5.51. Omitting it changes the reported likelihood, but not the parameter estimates.

Every REML route in DRModels.jl uses this normalised form, matching lme4, glmmTMB, TMB and drmTMB. See the fit_q4_reml docstring for the derivation and the limits of the cross-engine comparison.

source

DRModels.estimation_method Function

estimation_method(fit) -> Symbol

The estimator used to fit the model: :ML (default) or :REML (drm(...; method = :REML)).

source

DRModels.reml_objective_at Function

reml_objective_at(prob, Q_cond, phi; beta0=nothing, u0=nothing, n_newton=60)
    -> NamedTuple(reml_loglik, raw_reml_ll, beta, u, converged)

Evaluate the q=4 REML objective at a SUPPLIED phi = (beta_rho, lc) — the variance-component / correlation half of fit_q4_reml's parameter space — reprofiling beta_mu1/mu2/s1/s2 (the marginalised mean and scale fixed effects) by the SAME conditional-Newton alternation fit_q4_reml runs internally (reml_ll_and_mode), rather than re-running the outer LBFGS.

This is the diagnostic primitive behind cross-engine mode-finder-vs- objective-translation checks: given another engine's fitted point (mapped into DRModels.jl's phi/beta scale — see pack_phi, Λ_to_lc), call this to ask "is DRModels.jl's OWN objective, evaluated AT that point, better or worse than what DRModels.jl's own solver returned?" without hand-writing the inner E-step/Newton alternation each time.

beta0/u0 seed the conditional-Newton warm start (pass the other engine's fixed effects, or fit_q4_reml's returned .beta/.u_hat, as beta0/u0 respectively) so the inner profile has the best chance of reaching its true conditional optimum at phi rather than stalling from a cold start.

reml_loglik is on the SAME normalised (Patterson–Thompson, + (n_β/2)·log(2π)) scale as fit_q4_reml's .reml_loglik and as drmTMB/TMB/lme4/glmmTMB; raw_reml_ll is the pre-normalisation value, exposed for constant-offset sanity checks. n_β is prob's combined X1+X2+Xs1+Xs2 design width — the same marginalised-fixed-effect count fit_q4_reml uses.

Returns (reml_loglik, raw_reml_ll, beta, u, converged). converged is the inner alternation's own convergence flag (inner_converged from reml_ll_and_mode) — a barrier hit (non-PD Schur complement) surfaces as raw_reml_ll = reml_loglik = -Inf, converged = false, rather than an error.

source

StatsAPI.dof Function

dof(fit) -> Int

Degrees of freedom — the number of estimated parameters (length of θ).

source

StatsAPI.aic Function

aic(fit) -> Float64

Akaike information criterion, -2·loglik + 2·dof. Lower is better; compares models fit by ML (not REML) on the same data.

On a REML fit this uses the restricted log-likelihood and is only valid for comparing models that differ in variance structure only (same mean structure); a one-time warning is emitted. Use ML for cross-mean-structure selection.

source

StatsAPI.bic Function

bic(fit) -> Float64

Bayesian (Schwarz) information criterion, -2·loglik + dof·log(nobs).

On a REML fit this carries the same variance-only-comparison caveat as aic and emits a one-time warning.

source

DRModels.aicc Function

aicc(fit::DrmFit) -> Float64

Corrected Akaike information criterion (AICc) — the small-sample, second-order correction to aic:

AICc = AIC + 2k(k + 1) / (n − k − 1)

with k = dof(fit) estimated parameters and n = nobs(fit) observations. The correction is always positive, so aicc(fit) ≥ aic(fit), and it converges to aic(fit) as n → ∞. Prefer AICc over AIC when n / k is small (a common rule of thumb is n / k < 40). Like AIC, AICc compares models fit by ML on the same data; lower is better.

If n - k - 1 <= 0 (too few observations for the correction to be defined), returns Inf. On a VA fit this errors before that short-circuit: loglik carries an ELBO, not a marginal log-likelihood.

Example

fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data)
aicc(fit) > aic(fit)        # the correction is strictly positive
isfinite(aicc(fit))         # finite whenever n - k - 1 > 0

source

Model comparison

DRModels.lrtest Function

lrtest(reduced::DrmFit, full::DrmFit) -> NamedTuple

Likelihood-ratio test for two nested, ML-fitted models, mirroring drmTMB's anova(reduced, full). reduced must be a special case of full (fewer parameters); both must be fit by maximum likelihood (DRModels.jl's default — REML likelihoods are not comparable across fixed-effect structures).

Returns a NamedTuple (; statistic, dof, pvalue):

• statistic = 2 * (loglik(full) - loglik(reduced)) — the LR statistic, which is asymptotically χ² with dof degrees of freedom under the null that the reduced model is adequate.

• dof = dof(full) - dof(reduced) — the number of extra parameters in full.

• pvalue = ccdf(Chisq(dof), max(statistic, 0)) — the upper-tail χ² p-value.

dof must be positive (full must have more parameters than reduced), otherwise an ArgumentError is thrown. Both fits must have been estimated on the same observation rows, not only the same nobs: two fits that drop different rows and happen to keep the same count are refused (#1053). A negative statistic (the reduced model fits better, a sign the models are not actually nested, or one did not converge) is still returned as-is, but the p-value clamps the statistic at zero (so pvalue stays in [0, 1]); inspect statistic directly in that case.

Variance components and Poisson ⊂ NB2 are a boundary null

The χ²(dof) reference is only valid when the extra parameters in full are interior (regular). When full adds a variance component that reduced lacks (a random-effect SD, :resd/:resid; a Cholesky covariance entry, :recov; a group-level covariance, :phylocov), or when full is NegBinomial2() and reduced is Poisson() (the extra sigma tested at 0), the tested null sits on the boundary of the parameter space and the statistic follows a chi-bar-square mixture, not χ²(dof) (Self & Liang 1987; Stram & Lee 1994). The naive χ² p-value is then conservative (too large — the test loses power). lrtest detects these cases and emits a one-time warning; use lrt_boundary (or a parametric bootstrap) for a boundary-correct p-value. This function does not swap the reported p-value (#304). drmTMB's anova() does swap in the mixture for the same pairs; call lrt_boundary here to match that p-value.

Example

x = randn(400)
y = 0.5 .- 0.8 .* x .+ exp.(-0.3 .+ 0.4 .* x) .* randn(400)
data = (; y, x)

full    = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data)
reduced = drm(bf(@formula(y ~ 1),     @formula(sigma ~ 1)),     Gaussian(); data)

t = lrtest(reduced, full)
t.statistic    # 2·(logLik_full − logLik_reduced), > 0 when x helps
t.dof          # 2 extra parameters (x in μ and in log σ)
t.pvalue       # < 0.05 when x is truly predictive

source

DRModels.anova Function

anova(reduced::DrmFit, full::DrmFit) -> NamedTuple

Alias for lrtest, matching drmTMB's anova(reduced, full) spelling for a nested likelihood-ratio test. Returns the same (; statistic, dof, pvalue) NamedTuple.

Example

anova(reduced, full) == lrtest(reduced, full)   # true

source

StatsAPI.weights Function

weights(fit::DrmFit) -> Union{Vector{Float64}, Nothing}

Prior (per-observation) weights the fit was made with, from drm(...; weights = w), or nothing when no weights were given. Extends StatsAPI.weights.

The weights multiply each row's log-likelihood (drmTMB semantics), so loglik, aic, bic, profiles and likelihood-ratio tests use the weighted likelihood. nobs(fit) still counts rows, zero-weight rows included.

Prior weights are implemented on the fixed-effect ML routes of the univariate families. Any other route refuses weights = with an ArgumentError (#931; the family-by-family rollout to random-effect routes is tracked there).

Before #931 step (b), this accessor returned ones(nobs(fit)) for every fit, which looked like weights had been used.

Example

fit  = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1)), Gaussian(); data)
weights(fit) === nothing          # no prior weights given
fitw = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1)), Gaussian(); data, weights = w)
weights(fitw) == w

source

DRModels.update Function

update(fit::DrmFit, formula, args...; data, kwargs...) -> DrmFit

Refit fit with a new formula (a bf bundle), reusing the fitted family and every estimation option the fit stores:

• method — :REML when fit.estim_method === :REML (the default :ML is drm's own default and is not repeated);

• marginal — when it is not the default :LA;

• penalty — the stored drm_phylo_penalty specification of a :MAP fit.

Explicit keywords override those stored options (update(fit, formula; data, method = :ML) refits a REML seed by maximum likelihood). Changing method on a penalized MAP fit drops the stored penalty unless penalty is passed again; otherwise method = :ML would still be a MAP refit. family is not a refit keyword: the fitted family is kept, and passing one is an error. An unnamed extra positional argument is an error: a keyword-only refit would otherwise drop it and silently ignore the value.

data must be supplied. A DrmFit does not retain its observations, nor the options that were only arguments to the original drm call: K, A, tree, coords, algorithm, g_tol, profile_ci, phylo_coupled, sparse, impute, and missing. Pass those again as keywords when the refit needs them.

Example

full    = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1 + x)), Gaussian(); data)
# Drop x everywhere, keeping the same Gaussian family and the seed estimator:
reduced = update(full, bf(@formula(y ~ 1), @formula(sigma ~ 1)); data = data)
length(coef(reduced)) < length(coef(full))   # fewer parameters

source

Diagnostics and accessors

DRModels.check_drm Function

check_drm(fit) -> NamedTuple

Post-fit convergence / identifiability diagnostics — drmTMB's check_drm(). Returns a NamedTuple and logs a short report:

• convergence_status — convergence_status(fit): :converged, :boundary or :degenerate, by the same rules as drmTMB (FX3b).

• converged — is_converged(fit), i.e. convergence_status === :converged. A Gaussian fit whose scale collapsed is reported false here even when fit.converged is true (#949), and so is a :boundary fit (FX3b).

• max_abs_grad — max|∇nll| at the optimum (≈ 0 at a clean interior optimum; NaN when no gradient could be produced: see grad_source for which of the two reasons applies).

• grad_source names which producer supplied max_abs_grad, so a caller can tell an exact gradient from an approximate one and either from no gradient at all. Reuses the autodiff vocabulary of profile_result: :locscale (the canonical location-scale objective's exact analytic outer gradient), :stored (the fit's own gradient callback), :forward (ForwardDiff through the stored objective), :finite (a central finite difference of the stored objective, used when the objective is exact on Float64 but not dual-number safe, and accurate to roughly 1e-6 relative, NOT to machine precision), :none (the fit stores no objective, so there is nothing to differentiate), and :unavailable (an objective IS stored but neither its derivative nor a finite difference of it produced a finite value). :none and :unavailable are the two NaN cases and they mean different things; both leave ok unscored against the gradient criterion, and :unavailable also warns.

• vcov_complete — whether the stored covariance is finite throughout. Some routes report a partial covariance by design: the sparse phylo fitter computes the fixed-effect block and leaves the variance-component block NaN. When this is false the three fields below cannot be computed and are reported as false / NaN / Inf rather than raising.

• vcov_posdef — whether the stored covariance is positive-definite (drmTMB's sdreport is all-NaN exactly when this fails).

• min_eigval / cond — smallest eigenvalue and condition number of the covariance; a near-zero min_eigval flags a singular / weakly identified direction (e.g. a variance pinned at the boundary).

• penalized_map — whether this is a penalized (MAP) fit, using a phylogenetic penalty or fixed-effect priors. Standard errors describe local penalized curvature, and loglik is the unpenalized data log-likelihood.

• variance_boundary — nothing unless the fit is Gaussian with a grouped / structured random effect and a homoscedastic residual (sigma ~ 1); then the NamedTuple of _variance_boundary(fit) (residual_at_boundary, residual_ratio, structured_at_boundary, structured_ratios, one_obs_per_group). Flags a variance component at its lower boundary (the residual σ̂ of a phylogenetic fit whose structured term absorbs all the variance is an optimiser-stopping artefact, not an estimate — #724) and says when σ_a and σ_e are separated only by the covariance structure (one observation per group — #697). Does not affect ok.

• temporal_boundary — nothing unless the fit has a temporal() term; then (at_boundary, findings), drmTMB's temporal_boundary rules: residual σ̂ below 1e-3 · sd(y), OU decay × longest series span below 1e-4 or × shortest gap above 30, AR1 |φ̂| above 0.999. drmTMB reports such a fit as convergence_status() == "boundary" (so its is_converged is false since FX3b). Does not affect ok.

• skew_normal_nu_wald — nothing unless the fit is SkewNormal; then (near_zero, min_abs_nu, cut). near_zero is true when any fitted slant |ν̂| < 0.5 (cut): the information matrix is singular at ν = 0, so the Wald SE for nu is unreliable there (#1060), and a warning is emitted. The cutoff is an uncalibrated heuristic of this package; no Wald test is run. Does not affect ok.

• fixed_effect_rank_deficiency — nothing unless a fixed-effect design (mu or sigma) had an aliased column (sigma ~ 1 + x + z with z = 2x + 1, a duplicated column); then (deficient = true, blocks), one (param, rank, ncol, aliased) per rank-deficient block, naming the aliased coefficients (#1110; drmTMB's fixed_effect_rank_deficiency). Such a fit is not converged (is_converged is false), the aliased coordinates' SEs are NaN, and ok is false. Detected on the Gaussian location-scale route.

• ok — true when the optimiser converged at a non-degenerate optimum (the test is_converged made before FX3b), the gradient is small, the covariance is PD, and no fixed-effect design is rank-deficient. It does not read convergence_status, which is reported in its own field: a :boundary fit, or one the twin rules call :degenerate for an inflated standard error, can still be ok. On a phylogenetic-penalty fit the gradient criterion is dropped because its stored objective is the unpenalized likelihood. Fixed-prior fits retain the penalized objective and its gradient is checked. The criterion is also dropped whenever max_abs_grad is NaN, so read grad_source before reading ok: ok = true alongside :none or :unavailable means converged and PD only, with stationarity untested.

A non-ok result is informative, not an error: a model sitting on a variance boundary (Watanabe-singular) can be the data's MLE, with valid Wald SEs on the remaining directions — see confint.

source

DRModels.bridge_diagnostics Function

bridge_diagnostics(fit::DrmFit) -> NamedTuple

Route-aware convergence diagnostics for the engine = "julia" R bridge (#569) — the Julia twin of check_drm, reshaped for drm_bridge's payload plus the two quantities check_drm does not itself report: which internal route produced this fit and which integrator/optimiser it used.

Every field is read straight off fit, or computed by the SAME logic check_drm uses (_check_max_abs_grad, the covariance finiteness/positive-definiteness check) — deliberately NOT by calling check_drm(fit) itself, which additionally @info/@warn-logs a report on every call. drm_bridge calls this for every bridged fit, and a bridge boundary that writes to stderr on every ordinary fit is a regression in its own right (an R caller doing @test_nowarn drm_bridge(...)-equivalent checking, or just tailing its own logs, would see one check_drm report per engine = "julia" fit that nobody asked to be told about). Nothing here is fabricated. A quantity a route does not record is missing, never a fabricated zero or NaN standing in for it:

• route — the fitted objective's Julia type name (fit.nll === nothing reports "none"); the honest, ungeneralised answer to "which internal objective fitted this model".

• integrator — fit.marginal (:LA, :Laplace, :VA, :AGHQ).

• optimizer — "Optim.LBFGS" when niterations recorded an achieved iteration count (see its docstring for exactly which routes that covers); missing on a route with no single outer optimiser call to attribute one to.

• converged — the optimiser flag AND the degeneracy bar, the same value as the top-level bridge converged (#949). Since FX3b this is not is_converged, which is also false at a :boundary optimum; the bridge sends convergence_status as its own key.

• iterations — niterations(fit); missing when unrecorded (niterations returns -1).

• max_abs_grad, grad_source — the same (magnitude, source) pair check_drm reports as max_abs_grad/grad_source (_check_max_abs_grad); max_abs_grad is missing (not NaN) when grad_source is :none or :unavailable, i.e. no gradient was actually produced.

• vcov_complete — whether fit.vcov is finite throughout (check_drm's vcov_complete).

• vcov_posdef, min_eigval, cond — check_drm's positive-definiteness / eigenvalue / condition-number trio, but missing (not check_drm's documented false/NaN/Inf placeholders) whenever vcov_complete is false, since those three cannot actually be computed then.

• penalized_map — fit.estim_method === :MAP (check_drm's penalized_map).

• boundary — 1-based indices into fit.theta/diag(fit.vcov) whose stored variance is non-finite or negative — the same condition stderror already reports as an infinite standard error. Empty (not missing) when none are.

Example

d = bridge_diagnostics(fit)
d.route        # e.g. "LocScaleObjective"
d.grad_source  # e.g. :stored

source

DRModels.family Function

family(fit::DrmFit)

Return the response family object the model was fitted with, e.g. Gaussian(), Poisson(), Student(). This is the post-fit accessor for the family slot passed to drm; family(fit) === fit.family.

source

DRModels.convergence_status Function

convergence_status(fit) -> Symbol

Classify a fit as :converged, :boundary or :degenerate, by the same numeric rules as drmTMB's convergence_status() (FX3b; drmTMB docs/design/278-twin-convergence-status.md). The first matching rule wins:

1. :degenerate unless the optimiser stopped normally (fit.converged) at a finite log-likelihood strictly above the -1e15 failed-objective sentinel.

2. :degenerate when a geometry screen fired: a Poisson or NegBinomial2 fixed-effect MLE that does not exist (#1109), a fixed-effect-only Binomial fit with (quasi-)complete separation (#1186; drmTMB binomial_separation), or a Gaussian residual scale collapsed onto the data. The geometry flag is recorded at fit time, so it applies with se = false too.

3. :degenerate unless vcov(fit) is square, finite and positive definite.

4. :degenerate unless every standard error is finite.

5. :degenerate when a standard error is at least 50 and at least 1000 times the median of the standard errors below 50 (the smallest standard error when none is below 50): a near-flat, weakly identified direction.

6. :boundary when a fitted residual correlation has |rho12| > 0.98, or a fitted random-effect SD (re_sd) is below 1e-4. A random-effect log SD within 1e-3 of the sparse-Laplace floor of -8 also counts (the floor e^-8 = 3.4e-4 never falls under 1e-4). Its standard error is NaN (the Hessian of the other coordinates is inverted, #1190), and rules 3 to 5 then read the covariance of the remaining coordinates.

7. :converged otherwise.

Rules 3 to 5 are skipped when standard errors were turned off (se = false, or q4_vcov = false), as drmTMB skips them for drm_control(se = FALSE): such a fit is :converged unless rules 1 or 2 apply. A fit with an aliased fixed-effect column (#1110) is :degenerate through rule 2, as the design-rank rule of drmTMB.

is_converged and converged are true only for :converged. The raw optimiser flag is fit.converged.

source

DRModels.is_converged Function

is_converged(fit) -> Bool

true exactly when convergence_status(fit) === :converged, as in drmTMB. A false here means the reported estimates or standard errors should not be taken as an interior, identified optimum.

This is stricter than the raw fit.converged flag, which only says that the optimiser's own test was met. The Gaussian log-likelihood is unbounded as the residual scale goes to zero: with one row per group a structured random effect can interpolate the data, so sigma collapses and the objective runs away. A one-row-per-species phylogenetic fit produced sd_phylo = 22980, sigma = 7.5e-15, loglik = 6.8e13 with fit.converged = true, and 25% of parametric-bootstrap replicates in that example reached such a point. Since FX3b, a fit with a variance component or residual correlation at its boundary (:boundary) is also false here.

source

DRModels.converged Function

converged(fit) -> Bool

The same predicate as is_converged, under the name drmTMB also exports. The raw optimiser flag is fit.converged.

source

StatsAPI.deviance Function

deviance(fit::DrmFit) -> Float64

Deviance of the fitted model, -2 · loglik(fit) — drmTMB's deviance(). NaN (with a warning) on a fit stranded on the failed-objective sentinel. Extends StatsAPI.deviance.

source

StatsAPI.dof_residual Function

dof_residual(fit::DrmFit) -> Int

Residual degrees of freedom, nobs(fit) - dof(fit) (R's df.residual). Extends StatsAPI.dof_residual.

source

DRModels.r2_constant_sigma Function

r2_constant_sigma(fit::DrmFit) -> Float64

Coefficient of determination for a Gaussian mean model with constant σ:

R² = 1 - Σ(y - ŷ)² / Σ(y - ȳ)²

This is the ordinary lm() R², and on a constant-σ Gaussian fit it equals what lm() reports for the same mean model.

It refuses, rather than returning a number, whenever that formula would be ambiguous or misleading:

• a modelled σ (any non-intercept term in the sigma formula) — the total variance then depends on the covariates, so "the" denominator does not exist;

• random effects — the quantity wanted there is a marginal or conditional R² (Nakagawa & Schielzeth), which are two different numbers, neither of them this one;

• a non-Gaussian family, or a bivariate fit — a different definition applies, not this one;

• a mean model with no intercept — Σ(y - ȳ)² is then not the baseline the model is being compared against.

Examples

fit = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Gaussian(); data = d)
r2_constant_sigma(fit)      # e.g. 0.397

bad = drm(bf(@formula(y ~ x), @formula(sigma ~ sex)), Gaussian(); data = d)
r2_constant_sigma(bad)      # ArgumentError naming the modelled σ

See also coeftable, summary.

source

Fit summaries and route inventories

Base.summary Method

summary(fit::DrmFit)

Coefficient table for a fitted model — the DRModels.jl analogue of drmTMB's summary(). Returns the same CoefTable as coeftable. Fixed-prior MAP fits withhold z, p and interval columns; their SEs describe local penalized curvature.

source

DRModels.niterations Function

niterations(fit) -> Int

Optimiser iterations actually taken, or -1 when the fitter does not record it.

Deliberately NOT named iterations: Optim.iterations already means this, and DRModels.jl also uses iterations as a fitting OPTION (the cap). Keeping the accessor distinct stops "max allowed" and "actually taken" being confused for each other.

-1 is not a placeholder to be filled in later on every route — it is the honest answer for a fit that has no single outer optimiser call to count, and it is preferred over an approximated count.

Coverage by family

Reports Optim.iterations(res) from the LBFGS run that produced θ̂ for: Gaussian (both the plain ML fixed-effects fit and the Cox–Reid REML fixed-effects fit), Student, SkewNormal, Poisson, NegBinomial2, TruncatedNegBinomial2, Beta, BetaBinomial, Binomial, Gamma, LogNormal, ZeroOneBeta, Tweedie, CumulativeLogit — for their fixed-effects fit and, where the family has one, its scalar (1 | g) random-intercept, correlated (1 + x | g), zero-inflated (zi), hurdle (hu), and (Poisson only) AGHQ / coordinate-spatial-range variants. The bivariate residual routes (Gaussian/Gaussian, Student/Student, LogNormal/LogNormal) report counts the same way; LogNormal's bivariate fit borrows the Gaussian-on-log-y optimiser run wholesale (only the reported likelihood is Jacobian-shifted), so it carries that run's iteration count rather than re-deriving one. fit_mixed_family's cross-family latent-rho route reports its own optimiser run too, but through the returned NamedTuple's iterations field — that route does not produce a DrmFit, so niterations does not apply to it.

Returns -1 (no single outer LBFGS call to attribute the count to, or no recorded count): Gaussian's meta_V, phylo/relmat/animal/spatial, and multi-random-effect routes (the Cox–Reid REML random-intercept route, e.g. Poisson (1 | g) with method = :REML, also stays -1 — its reported θ̂ comes from a secondary restricted refit, not the counted LBFGS run, so attributing that run's count to it would be a mismatch, not a full count); the bivariate Gaussian phylo/structured (q2/q4) sparse-Laplace routes; and every family's phylo/relmat/animal/coordinate-spatial random-effect routes other than Poisson's spatial-range fit above. These share the sparse augmented-state Laplace engine (src/sparse_*.jl, src/*_phylo.jl) rather than a single top-level Optim.optimize call, so there is no one iteration count to report honestly; do not infer non-iteration (e.g. "closed form") from -1 on these routes — check the family/route, not just the flag. The ordinary (1 | g) marginal = :Laplace route (Poisson, Binomial, NegBinomial2, Gamma, Beta) is also -1: its θ̂ comes from a chain of optimiser stages (LBFGS, a short polish, a boundary polish and, when needed, Newton steps), not one counted run.

source

DRModels.profile_targets Function

profile_targets(fit::DrmFit; ready_only = false) -> Vector{NamedTuple}

Every parameter profile_result / confint(..., method = :profile) can be asked for on this fit, with an honest readiness flag — drmTMB's profile_targets().

Runs no optimisation: it walks the fitted object. One row per coefficient with

• parm — the coefficient name;

• param — its block (:mu, :sigma, :resd, …);

• index — its position in fit.theta;

• estimate — the fitted value on the estimation scale;

• scale — :log for a variance-component / scale coefficient (and the temporal OU decay), :atanh for the temporal AR1 persistence and the Toeplitz partial autocorrelations, :identity otherwise;

• profile_ready — whether a profile interval can actually be computed here;

• profile_note — why, when it cannot.

Pass ready_only = true to drop the unavailable rows.

Example

tg = profile_targets(fit)
filter(r -> !r.profile_ready, tg)        # what will refuse, and why

source

DRModels.structured_effects Function

structured_effects(fit::DrmFit) -> Vector{NamedTuple}

One row per structured marker in the fitted formula — drmTMB's structured_effects(). Fields dpar, kind, grouping.

kind is the marker (:phylo, :relmat, :animal, :spatial, :temporal), grouping the factor it wraps, and dpar the distributional parameter whose formula carried it. Exists so downstream code never has to grep or re-parse formula text.

Returns an empty vector for a model with no structured markers, and for a fit whose formula was not retained.

Example

structured_effects(fit)
# 2-element Vector{NamedTuple}:
#  (dpar = :mu,    kind = :phylo, grouping = :species)
#  (dpar = :sigma, kind = :phylo, grouping = :species)

source

Derived phylogenetic quantities and boundary comparisons

DRModels.gaussian_locscale_phylo_sds Function

gaussian_locscale_phylo_sds(fit::DrmFit) -> NamedTuple

Extract sd(μ-phylo) and sd(σ-phylo) from a SEPARATE-block Gaussian phylo location-scale fit. Reads the :resd_mu / :resd_sigma blocks (exp of their stored log values). For a COUPLED fit use fit.scales[:lambda_sd_mu] etc.

source

DRModels.profile_sigma_a Function

profile_sigma_a(fit; level=0.95, axes=:all, chibar=false, n_newton=nothing,
                g_tol=1e-3, max_bisect=14) -> NamedTuple

Profile-likelihood CIs for the among-axis SDs sqrt.(diag(Σ_a)) of a bivariate q=4 phylogenetic location–scale fit. Returns a NamedTuple with summary rows (param, coef, estimate, lower, upper, deviance_floor, bounded) for param ∈ (:sd_mu1, :sd_mu2, :sd_sigma1, :sd_sigma2). lower=0 ⇒ the axis is at / near the boundary; upper=Inf with bounded=false ⇒ no upper crossing was found (the profile is too flat — report it honestly, do not invent a bound). Also returns level, chibar, threshold, nll_hat, axes.

Threshold: chibar=false (default) uses the standard profile-LR cutoff χ²₁(level). For an identified (interior) SD this gives ~nominal coverage; for a collapsed (boundary) axis it is conservative (over-covers — the LR statistic is then ≈0.5·χ²₀+0.5·χ²₁, so the exclusion rate is ≈0.5·(1−level), e.g. a 90% naive CI covers ~95% on the boundary) and the lower bound is 0. This conservatism is the safe failure mode for an axis of unknown status, so naive is the publication default. chibar=true uses the boundary mixture cutoff χ²₁(2·level−1) — nominal only when the axis is known to be at the boundary AND the other variance/correlation nuisances are interior (a single boundary component); with several axes near their boundaries the 0.5:0.5 mixture is itself approximate. Needs no Hessian; the fit must carry fit.ranef.prob. The calibrated complement to [bootstrap_sigma_a].

source

DRModels.bootstrap_sigma_a Function

bootstrap_sigma_a(fit; data, B = 300, level = 0.95, rng = default_rng(),
                  failures = :warn, check_converged = true,
                  q4_g_tol = 1e-3, q4_iterations = 300) -> NamedTuple

Parametric bootstrap of the among-axis standard deviations sqrt.(diag(Σ_a)) for a bivariate q=4 phylogenetic location–scale fit (drm(bf(mu1=…, mu2=…, sigma1=…, sigma2=…, rho12=…), Gaussian(); tree = …)).

Returns a NamedTuple with

• summary — a vector of rows (param, coef, estimate, std_error, lower, upper) for param ∈ (:sd_mu1, :sd_mu2, :sd_sigma1, :sd_sigma2). estimate is the fitted SD; lower/upper are the level percentile interval over the successful replicates.

• cor_summary — the same rows for the 6 among-axis correlations (:cor_mu1_mu2, :cor_mu1_sigma1, :cor_mu1_sigma2, :cor_mu2_sigma1, :cor_mu2_sigma2, :cor_sigma1_sigma2) — the coevolutionary correlations of coevolution_cor, now with CIs. A correlation whose axis collapses is unidentified and comes back with a wide interval (e.g. ρ_a(μ1,μ2) spanning the sign when a σ-axis pins).

• status — "bootstrap" when every replicate is kept, "bootstrap_incomplete" when any is dropped and at least two remain, "bootstrap_at_boundary" when the retained draws also sit on an SD or correlation bound, or "bootstrap_unavailable" when fewer than two remain. A dropped replicate warns.

• attempted, used, failed, failures, level, draws (the used × 4 matrix of replicate SDs), cor_draws (used × 6), axes, cor_pairs, elapsed.

The interval is the honest report at a boundary: when a scale axis carries no phylogenetic signal its SD collapses and the percentile interval includes ~0, where the Hessian-based SE / profile is undefined. data must be the table the fit was estimated on (the covariate / grouping columns are reused; the two response columns are overwritten per replicate).

With failures = :warn (default) failed refits are dropped and reported in failures; failures = :error rethrows on the first failure.

source

DRModels.repeatability Function

repeatability(fit; component = nothing, level = 0.95, method = :delta) -> NamedTuple
repeatability(fit, newdata; level = 0.95) -> NamedTuple

Alias for icc: the adjusted repeatability R = σ²_g / (σ²_g + σ²_resid) for the chosen grouping factor. With a single structured component and no other components, repeatability and heritability coincide.

A random intercept on the scale

A random intercept on sigma (sigma ~ (1 | g), SD reported as re_sd(fit)[:g_logsigma]) lives on the log σ scale. It is not a variance component of the response, so it is never a component and never enters the numerator. It makes the residual variance vary by group, so the residual entry of the denominator becomes the marginal residual variance

E[σ²] = exp(2 b₀ + 2 Σ_k ω_k²),

where b₀ is the sigma intercept and ω_k the log-σ random-intercept SDs (not the squared median exp(2 b₀)); the delta-method gradient runs through ω_k too. This is drmTMB's definition. A fit whose only random effect is on sigma has no mean component and is refused.

Location–scale–scale fits: the estimand is conditional on covariates

When the between-individual SD and the residual SD both depend on covariates (sd(id) ~ z, sigma ~ z),

y_ij ~ N(μ_ij, σ_e(z_i)²),   b_i ~ N(0, σ_b(z_i)²),
log σ_b(z) = α'z,            log σ_e(z) = γ'z,

repeatability is a FUNCTION of the covariates, not a number:

R(z) = σ_b(z)² / (σ_b(z)² + σ_e(z)²) = logistic( 2 (α'z − γ'z) ),

the correlation between two observations of an individual whose covariates are z. No single scalar is the repeatability of such a model — any one number is a choice of covariate distribution (and the ratio of average variances is NOT the average of the ratios) — so repeatability(fit) / icc(fit) / heritability(fit) refuse these fits rather than pick one silently. Use the two-argument form: newdata is a column table holding every predictor of the sd(g) and sigma formulas (categorical predictors must contain all levels, in the training coding, so the design matches). It returns, per row of newdata,

• estimate — R(z);

• se_logit — the Wald SE of logit R(z) = 2(α'z − γ'z), which is LINEAR in the coefficients, so it is exact up to the usual Wald approximation (uses the joint vcov of the sd and sigma blocks, so the covariance between them counts);

• lower, upper — the Wald interval on the logit scale mapped back to (0, 1);

• level, method = :wald_logit.

fit = drm(bf(@formula(y ~ sex + (1 | id)), @formula(sigma ~ sex),
             @formula(sd(id) ~ sex)), Gaussian(); data)
repeatability(fit, (; sex = [0.0, 1.0]))   # R for females, R for males

sd(g, phylogenetic) ~ z fits are handled by the same call to heritability: heritability(fit, newdata) returns h²(z) = σ_a(z)² / (σ_a(z)² + σ_e(z)²), the tip-level share of variance that is phylogenetic.

source

DRModels.lrt_boundary Function

lrt_boundary(fit_full::DrmFit, fit_reduced::DrmFit; q::Integer = 1) -> NamedTuple

Boundary-corrected likelihood-ratio test for q variance components = 0, comparing two nested, ML-fitted models. fit_reduced must be fit_full with the q variance component(s) removed (e.g. dropping a random effect (1 | g)). Unlike lrtest — which uses the naive χ²(q) reference and is conservative for variance-component tests — this uses the chi-bar-square mixture appropriate to a parameter on the boundary of its space.

Returns a NamedTuple (; statistic, q, pvalue, pvalue_naive):

• statistic = 2 · (loglik(fit_full) − loglik(fit_reduced)) — the LR statistic.

• q — the number of boundary variance components tested (the mixture order).

• pvalue = chibar_pvalue(statistic, q) — the χ̄² boundary p-value (see chibar_pvalue).

• pvalue_naive = ccdf(Chisq(q), max(statistic, 0)) — the naive χ²(q) p-value, for comparison. Always pvalue ≤ pvalue_naive (the correction makes the test less conservative, i.e. more powerful), so reporting both is informative.

Assumptions

Same as chibar_pvalue: the q dropped parameters are variances tested at 0; for q = 2 they are independent; all other parameters are interior; ML fits. A negative statistic clamps to the boundary p-value.

Finite-sample behaviour: Gaussian (1 | g) with few groups

The 50:50 mixture is the asymptotic null law, and it is accepted as is for the Gaussian random-intercept test lrt_boundary(full, reduced; q = 1) (dropping (1 | g)). With few groups it is conservative: in 500 null replicates (sd_b = 0, 25 groups of 12, y = 0.4 - 0.5 x + 0.7 ε, nominal α = 0.05) the chi-bar p-value rejected in about 3.4% of replicates (17/500), against 1.2% (6/500) for the naive lrtest χ²₁ p-value; 58.6% of the statistics sat on the χ²(0) atom. The statistics match drmTMB on the same data to about 1e-8, so this is the behaviour of the reference distribution, not of the fit. A pvalue just above 0.05 is therefore weak evidence for dropping the random effect, and no mixture weights are changed.

The same refusals as lrtest apply before the statistic is formed: REML fits whose fixed structure differs, a sentinel log-likelihood, different marginal approximations, a penalized MAP fit, and fits on different observation rows (not only a different nobs).

These guards are the same set lrtest runs (#1002; pinned in test/test_lrt_boundary_guards_1002.jl). They do not read converged, and lrt_boundary does not check the family or that fit_reduced is truly nested in fit_full; a non-converged fit that is not a sentinel still gets a p-value.

Finite-sample behaviour: Poisson vs NegBinomial2

The 50:50 mixture is the asymptotic null law, and it is accepted as is for the family nest lrt_boundary(nb, pois; q = 1) (Poisson ⊂ NegBinomial2, where the extra dispersion sigma sits on the boundary). In finite samples it is conservative for this pair: in 500 Poisson-null replicates (n = 250, y ~ Poisson(exp(0.4 + 0.3 x)), nominal α = 0.05) the 50:50 p-value rejected in about 3% of replicates (15/500), against 1.8% (9/500) for the naive lrtest χ²₁ p-value. So a pvalue just above 0.05 here is not evidence for the Poisson model, and the test is somewhat low-powered near the boundary; no mixture weights are changed.

Example

# Random intercept vs no random effect: dropping (1 | g) removes ONE variance.
full    = drm(bf(@formula(y ~ x + (1 | g)), @formula(sigma ~ 1)), Gaussian(); data)
reduced = drm(bf(@formula(y ~ x),           @formula(sigma ~ 1)), Gaussian(); data)

t = lrt_boundary(full, reduced; q = 1)
t.statistic        # 2·(ℓ_full − ℓ_reduced)
t.pvalue           # χ̄²: 0.5·P(χ²₁ > stat)   — correct for the boundary
t.pvalue_naive     # χ²(1): P(χ²₁ > stat)     — conservative (≈ 2× too large)

source

R bridge and preprocessing

These entries document the Julia-side interface used by the optional R bridge. Coming from R explains the optional route and refusals; these docstrings do not expand that contract.

DRModels.drm_bridge Function

drm_bridge(; formula, family, data, tree = nothing, K = nothing,
           A = nothing, Ainv = nothing, coords = nothing, newdata = nothing,
           options = Dict(), K_labels = nothing, A_labels = nothing,
           structured_group = nothing)

Fit a DRModels.jl model through a marshalling-friendly boundary for R callers. formula may be a semicolon-separated string such as "y ~ x; sigma ~ x" or a dictionary / named tuple whose values are formula strings. family is a string such as "gaussian", "student", "nbinom2", or "biv_gaussian". data is a column table, dictionary, or named tuple.

formula accepts R syntax beyond plain @formula: : interactions, * crossing, - 1/general - term removal, (...)^k crossing, and scale(x)/I(expr)/factor(x)/poly(x, k) (materialised into real columns — I(...) only over a safe + - * / ^ grammar, never arbitrary code). poly(x, k) is R's ORTHOGONAL basis (raw = FALSE, the default) and expands to k columns; only poly(x, k) itself is accepted — raw = TRUE is spelled I(x^k), and an explicit coefs = or the multivariate poly(x, y, degree) are rejected rather than approximated. These materialised columns are NOT (yet) reconstructed for newdata — a formula using them combined with newdata fails with a missing-column error.

K / A are matched to the levels of the grouping variable by position (first appearance in data). To match by name instead, pass K_labels / A_labels (the matrix's row/column names) together with structured_group (the grouping column, as a string): the matrix is reordered to first-appearance order, and a duplicated label or a missing group level throws ArgumentError rather than fitting against the wrong rows.

The return value is a Dict{String,Any} made of primitive R-reconstructable pieces: named coefficients, covariance matrix, likelihood summaries, fitted values, residuals, scales, and residual correlations when present.

It also carries what drmTMB's post-fit surface consumes: "dpars", the per-observation distributional parameters on the response scale keyed by dpar name, plus "trials" for binomial-type families. Pass newdata (a column table) to add "dpars_newdata", the same parameters evaluated on fresh rows — what fitted_distribution(object, newdata = ...) needs.

For animal models, Ainv supplies a precision matrix directly. It is passed to the family fitter without covariance normalization; only routes that explicitly accept animal precision matrices can use it.

source

DRModels.drm_bridge_inference Function

drm_bridge_inference(; formula, family, data, tree = nothing,
                     K = nothing, A = nothing, coords = nothing,
                     options = Dict(), method = "profile",
                     level = 0.95, B = 199, seed = nothing,
                     threads = false, parm = nothing)

Run a narrow inference primitive for the R bridge.

With parm = nothing (the default) this targets the Gaussian phylogenetic SD block (param = :resd / :resd_mu / :resd_sigma). Pass parm = "sd:mu" or "sd:sigma" to select the coupled location- or scale-axis row explicitly, or "sd:resd" / "sd:resd_mu" / "sd:resd_sigma" when an R bridge caller must preserve the fitted parameter block exactly.

Pass parm = "fixef:<dpar>:<coef>" (e.g. "fixef:mu:x") to instead profile or bootstrap a single ordinary fixed-effect coefficient, on its link scale, using DRModels.profile_result or DRModels.bootstrap_result. Returns the same payload shape either way. For an explicit structured fixed-effect target, the supplied covariance provider (tree, K, A, or coords) is reused for the initial fit, marginal simulation, and every bootstrap refit.

source

DRModels.drm_bridge_objective_at Function

drm_bridge_objective_at(formula, family, data, tree, options = Dict();
                        beta, Lambda, rho12)

The bivariate q=4 phylogenetic REML counterpart to reml_objective_at, using the same R-compatible data interface as drm_bridge. The drmTMB R wrapper is drm_julia_reml_objective_at(). formula, family, data, tree, options are exactly the payload drm_bridge takes for a bivariate q=4 phylogenetic model (a formula with phylo(...) markers shared by mu1, mu2, sigma1, sigma2); options is accepted for positional parity with that payload but is not otherwise used — this is a read-only diagnostic, not a fitting route.

beta, Lambda, rho12 are the outer evaluation POINT, in the SAME (fixed-effect, among-axis-covariance, residual-correlation) parameterisation drm_bridge's q4 phylo REML fit reports: beta is a NamedTuple or Dict with mu1/mu2/sigma1/sigma2 numeric vectors (inner-Newton warm starts for the profiled-out fixed effects — reml_objective_at reprofiles them at phi regardless of the warm start supplied), Lambda is the 4×4 symmetric among-axis covariance matrix (axis order mu1, mu2, sigma1, sigma2), and rho12 is the residual correlation. Internally: Lambda/rho12 are packed into DRModels.jl's own phi = (beta_rho, lc) via pack_phi and reml_objective_at evaluates the q=4 REML objective there.

Returns a Dict{String,Any} with "objective" and "reml_loglik" (the normalised Patterson–Thompson restricted log-likelihood reml_objective_at reports — the two keys carry the same value; "objective" is the route-agnostic name, "reml_loglik" names the DRModels.jl convention explicitly), "raw_reml_ll" (the pre-normalisation value), "converged_inner" (the inner conditional-Newton alternation's own convergence flag — a barrier hit surfaces as -Inf/false rather than an error), and "contract" => "bridge_objective_at_v1" so R callers can assert the return shape.

This diagnostic evaluates a supplied point; it does not fit or select a model. It helps distinguish differences in objective values from differences in the optima found by separate engines. Agreement at one point does not establish general cross-engine equivalence.

source

DRModels.drm_listwise Function

drm_listwise(formula, data; verbose = true) -> NamedTuple

Drop every row of data that has a missing or non-finite (NaN/Inf) value in any response or predictor column referenced by formula (a bf bundle or the bivariate keyword form), returning a column NamedTuple ready to pass straight to drm. This is listwise deletion (a.k.a. complete-case analysis).

clean = drm_listwise(bf(@formula(y ~ x), @formula(sigma ~ x)), data)        # @warn names the dropped rows
fit   = drm(bf(@formula(y ~ x), @formula(sigma ~ x)), Gaussian(); data = clean)

Why this exists: drm itself drops rows with a missing predictor (and a missing response where the route has no observed-rows path), with a warning that gives the count. drm_listwise does the same drop explicitly, also drops NaN/Inf cells, and emits one @warn reporting how many rows (and which columns) were removed, so the data you pass to drm are exactly the rows it fits.

Listwise deletion discards information

For the bivariate / q=4 coevolution models a row is dropped if either trait is missing — exactly the information that full-information maximum

likelihood (FIML) is designed to keep. Listwise deletion is unbiased only under missing-completely-at-random (MCAR) and is generally less efficient than FIML under missing-at-random (MAR). Native FIML for missing responses and multiple imputation for missing predictors are not currently available; this helper is an interim complete-case path, not a substitute.

Only the columns the model uses are checked, so unrelated columns with missing values do not cause rows to be dropped. Set verbose = false to suppress the warning (the drop still happens). If no rows are dropped the columns are returned as they came, except that a Union{Missing,T} column that holds no missing is narrowed to T (as it is on the dropping branch), so it fits as numbers.

source

Staged pair association diagnostics

DRModels.PairAssociation Type

PairAssociation

Result of associate_pairs: a staged, frozen-margin association fit.

Fields: pair_class, coefficients (association scale, alpha), eta (0.999999·tanh(alpha)), loglik, score and curvature (finite-difference diagnostics at the optimum), near_boundary, multistart_disagreement, and nobs.

The uncertainty is conditional on the frozen margins. See association.

source

DRModels.integration_diagnostics Function

integration_diagnostics(a::PairAssociation)

Per-row quadrature quality for a both-censored staged fit, or nothing for the closed-form classes.

Returns the rectangle probability and QuadGK's absolute error estimate per row, plus the worst relative error. drmTMB retains abs.error from its own adaptive integration for exactly this reason: a rectangle probability that silently lost precision would corrupt the association without any visible failure.

source

Phylogenetic penalty (MAP)

drm_phylo_penalty is the Julia twin of drmTMB's drm_phylo_penalty(): pass penalty = drm_phylo_penalty(...) to drm to turn a phylogenetic variance-component fit into a MAP estimate. drm_phylo_penalty_sweep refits across a grid of correlation-penalty scales so you can see whether a conclusion moves with the prior. PhyloCorPenaltyNeedsTwoSD is the typed refusal when cor_sd is asked of a model that has no phylogenetic correlation to penalize.

DRModels.drm_phylo_penalty Function

drm_phylo_penalty(; sd_u = 1.0, sd_alpha = 0.05, cor_sd = nothing)

A PC-prior-style penalty on phylogenetic variance components — the Julia twin of drmTMB's drm_phylo_penalty(). Passing it to drm as penalty = … turns the fit into a MAP estimate.

• sd_u, sd_alpha: an exponential prior on each phylogenetic SD, calibrated so that a priori P(sd > sd_u) = sd_alpha. Smaller sd_u shrinks harder.

• cor_sd: optional SD of a mean-zero normal prior on the unconstrained (Fisher-z) phylogenetic correlation. Only meaningful for a route that actually estimates a correlation — the coupled mean↔σ block. Requesting it elsewhere raises PhyloCorPenaltyNeedsTwoSD rather than silently doing nothing.

The penalized fit reports an unpenalized loglik; the penalty value is on fit.phylo_penalty. Standard errors come from the penalized curvature and are credible-interval-shaped — do not read them as frequentist, and do not compare penalized fits with aic / lrtest (both refuse).

sd_u is measured on YOUR tree's scale

sd_u is a threshold on the phylogenetic SD, and that SD is only meaningful relative to the tree's scale. DRModels.jl builds its phylogenetic covariance from the branch lengths as supplied, so a tree of height h has tip variance h and the fitted SD carries a factor sqrt(h). drmTMB instead standardises via ape::vcv(tree, corr = TRUE), whose tips always have variance 1.

(The Gaussian phylo(1 | g) mean route with a constant sigma on an ultrametric tree is the exception: it fits at unit height and reports the correlation scale, so sd_u is already on drmTMB's scale there, #979.)

So the same sd_u is the same prior on both sides only when the tree has unit height. Copying sd_u = 1 from an R script onto a raw tree of height 2 gives a prior roughly sqrt(2) tighter than intended. Rescale first:

# ape: tree$edge.length <- tree$edge.length / max(diag(vcv(tree)))

Measured: on an ape::rcoal tree of height 1.5285 the two SDs differed by exactly sqrt(1.5285) = 1.2363 while the log-likelihoods agreed to five decimals — the same fit, two scales. See tools/parity_phylo_penalty.R.

Example

pen = drm_phylo_penalty(sd_u = 0.5, sd_alpha = 0.05)
fit = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
          Gaussian(); data = dat, tree = tree, penalty = pen)
fit.phylo_penalty        # the penalty at the optimum
loglik(fit)              # UNPENALIZED data log-likelihood

source

DRModels.drm_phylo_penalty_sweep Function

drm_phylo_penalty_sweep(f, fam; data, cor_sd = [0.25, 0.5, 1.0],
                        sd_u = 1.0, sd_alpha = 0.05,
                        phylo_coupled = true, kwargs...)

Prior-sensitivity sweep over the phylogenetic-correlation penalty — the Julia twin of drmTMB's drm_phylo_penalty_sweep(). Refits the model once per cor_sd value and reports whether the conclusion moves as the prior tightens.

Returns (; summary, fits). summary is a Vector{NamedTuple} with fields cor_sd, converged, vcov_posdef, loglik, cor, and error — a failed fit contributes a row of NaN/missing plus its message rather than aborting the sweep. fits is a Dict keyed "cor_sd=<value>".

phylo_coupled = true is the default because only the coupled mean↔σ block estimates a phylogenetic correlation at all; sweeping cor_sd over any other block would produce identical rows and look like a sensitivity check while being a no-op. That case raises PhyloCorPenaltyNeedsTwoSD from a probe fit before any of the sweep runs, so the no-op is reported rather than sold.

Extra kwargs... are forwarded to drm (tree, g_tol, …).

Example

sw = drm_phylo_penalty_sweep(
        bf(@formula(y ~ x + phylo(1 | species)),
            @formula(sigma ~ 1 + phylo(1 | species))),
        Gaussian(); data = dat, tree = tree, cor_sd = [0.25, 0.5, 1.0])
sw.summary          # one row per cor_sd

source

DRModels.PhyloPenalty Type

PhyloPenalty

A penalty specification for phylogenetic variance components — drmTMB's drm_phylo_penalty() object. Fields sd_u, sd_alpha, rate, cor_sd.

rate is derived, not supplied: rate = -log(sd_alpha) / sd_u, so that a priori P(sd > sd_u) = sd_alpha.

source

DRModels.PhyloCorPenaltyNeedsTwoSD Type

PhyloCorPenaltyNeedsTwoSD(msg)

Raised when cor_sd is requested on a model that has no phylogenetic correlation parameter to penalize. It is a distinct type, not a bare ArgumentError, because drm_phylo_penalty_sweep must catch precisely this condition — a cor_sd sweep over a model with no correlation would return identical rows and look like a prior-sensitivity check while being a no-op.

Mirrors drmTMB's classed condition drm_phylo_cor_penalty_needs_two_sd.

source

Staged pair association

associate_pairs estimates a latent-normal association between two already-fitted univariate models (drmTMB's staged, frozen-margin route). The kernel must be given explicitly as latent_normal; association is the post-fit summary.

DRModels.associate_pairs Function

associate_pairs(fit_1, fit_2; kernel, association = nothing, marginal = :LA)

Estimate a staged latent-normal association between two already-fitted univariate models — drmTMB's associate_pairs().

kernel must be given explicitly as latent_normal(). association defaults to an intercept (~ 1); only the intercept-only form is implemented in this slice, matching drmTMB's boundary for every pair class except its one Bernoulli × NB2 exception.

The margins are frozen: fit_1 and fit_2 are used as fitted and are never re-estimated. Reported uncertainty is therefore conditional on them.

Implemented pair classes are Gaussian × Bernoulli, Gaussian × NB2, Bernoulli × Bernoulli, Bernoulli × NB2, and NB2 × NB2. A Binomial margin must represent Bernoulli observations on these routes. The accepted classes have class-specific continuous/discrete likelihood calculations; other family pairs are refused. These implementation routes do not imply general interval calibration, and the margins remain frozen as described above.

fg = drm(bf(@formula(y ~ x), @formula(sigma ~ 1)), Gaussian(); data = d)
fb = drm(bf(@formula(z ~ x)), Binomial(); data = d)
a  = associate_pairs(fg, fb; kernel = latent_normal())
a.eta

source

DRModels.latent_normal Function

latent_normal()

The latent-normal association kernel. Must be passed explicitly: like drmTMB, there is no implicit association model.

source

DRModels.association Function

association(a::PairAssociation)

Association summary for a staged fit — drmTMB's association().

Returns a NamedTuple with eta (the latent correlation), alpha (the association-scale coefficient), score/curvature diagnostics, near_boundary, and multistart_disagreement.

The uncertainty is conditional on the frozen margins — it ignores margin estimation error, so it is not a joint standard error. Matching drmTMB, no simultaneous eta bands or profile intervals are offered.

source

Heritability, ICC, and derived quantities

DRModels.heritability Function

heritability(fit; component = nothing, level = 0.95, method = :delta) -> NamedTuple

Phylogenetic heritability / signal (a.k.a. λ / H²) of a structured-Gaussian fit: the share of the total variance carried by one structured component,

h² = σ²_component / ( Σ_k σ²_k + σ²_resid ),

where the sum runs over all structured variance components plus the residual. This is the comparative-biology "phylogenetic signal" — for a single phylo(1 | species) component it is Pagel/Lynch's phylogenetic heritability; with a second structured component (e.g. + animal(1 | id)) the denominator includes it too.

component selects which grouping factor is the numerator (a Symbol, e.g. :species); if omitted and the fit has exactly one structured component, that one is used. method is :delta (epsilon-method / generalized-delta via bias_correct, the default) or :profile (a true profile-likelihood CI on the ratio: at each fixed ratio the likelihood is re-maximised over ALL nuisance parameters — the other variance components, residual, and mean coefficients — so the co-components can absorb variance and the profiled deviance rises at the correct rate; it is NOT a substitution/ELR profile that freezes the nuisance variances at the MLE). For the dense phylogenetic correlation-scale parameterisation, fit with algorithm = :gls before using delta-method Wald intervals; the default sparse all-node phylogenetic route stores only partial covariance information in this slice, so profile intervals are the safer uncertainty path there.

For an active fixed-effect prior fit, variance-ratio intervals are not qualified; heritability, icc, and repeatability refuse both interval methods.

Returns a NamedTuple:

• estimate — the plug-in ratio g(θ̂) (exactly in [0, 1]);

• corrected — the bias-corrected estimate g(θ̂) + ½·tr(H_g·V) (delta only);

• se — the delta-method standard error (delta only);

• ci — the (lower, upper) CI, clamped to [0, 1] (when the delta-method SE is not finite the bounds are (NaN, NaN) (an all-NaN covariance or a finite but non-positive-definite one that gives a negative variance); read status, which is then :vcov_unavailable);

• level — the confidence level;

• method — the method used;

• status — :ok, or the first applicable caveat: :not_converged (the optimiser did not stop normally, or the optimum is degenerate, which includes a sentinel log-likelihood; a variance component on its boundary, convergence_status :boundary, is not flagged here), :penalized_map (a MAP fit, whose variance components are deliberately shrunk), :vcov_unavailable (delta SE not finite), :profile_edge (a profile bound is the 0/1 edge, i.e. one-sided, not a likelihood-ratio crossing) or :profile_unresolved (a profile bound is NaN). Every caveat except :ok also warns.

The gradient and Hessian of the ratio are threaded EXACTLY through the log σ → variance map by automatic differentiation. At a variance boundary (σ_component → 0 ⇒ h² ≈ 0, or σ_resid → 0 ⇒ h² ≈ 1) the Wald SE can be degenerate; the profile method gives a more honest (possibly one-sided) interval there.

Example

fit = drm(bf(@formula(y ~ x + phylo(1 | species)), @formula(sigma ~ 1)),
          Gaussian(); data, tree, algorithm = :gls)
h = heritability(fit)             # single component ⇒ no `component` needed
h.estimate, h.ci

Location–scale–scale fits

With sd(g) ~ z / sd(g, phylogenetic) ~ z and sigma ~ z the ratio depends on the covariates, so this form refuses. Use heritability(fit, newdata; level = 0.95) (phylogenetic sd) or repeatability(fit, newdata) (iid sd), which return the covariate-conditional R(z) = σ_b(z)² / (σ_b(z)² + σ_e(z)²) per row of newdata with a Wald-on-logit interval; see repeatability for the definition.

A random intercept on sigma is not a component: it turns σ²_resid into the marginal residual variance E[σ²]; see repeatability.

source

DRModels.icc Function

icc(fit; component = nothing, level = 0.95, method = :delta) -> NamedTuple

Intraclass correlation / repeatability for one grouping factor,

ICC = σ²_component / ( σ²_component + σ²_resid ),

the share of variance at the grouping level relative to that component plus the residual (the classic two-component repeatability). When the fit has more than one structured component this is the focal-vs-residual repeatability for the chosen component; use heritability for the full-variance share that also nets out the other components. Same return shape and method options as heritability; the CI is clamped to [0, 1].

For location–scale–scale fits (sd(g) ~ z) this form refuses; use icc(fit, newdata) — the covariate-conditional repeatability documented under repeatability.

A random intercept on sigma is not a component: it turns σ²_resid into the marginal residual variance E[σ²]; see repeatability.

source

DRModels.bias_correct Function

bias_correct(fit, g::Function; level = 0.95) -> NamedTuple

Epsilon-method (generalized-delta) bias correction for a smooth scalar derived quantity g(θ) of a fitted model, à la TMB sdreport(..., bias.correct = TRUE) (Thorson & Kristensen 2016).

g receives the full coefficient vector θ = coef(fit) (working scale: μ on the response scale, σ on log σ, ρ12 on atanh ρ12, random-effect SDs on log σ_b) and returns a real scalar.

Returns a NamedTuple with

• estimate — the raw plug-in value g(θ̂);

• corrected — the bias-corrected value g(θ̂) + ½·tr(H_g·V);

• bias — the correction itself, ½·tr(H_g·V) (so corrected = estimate + bias);

• se — the delta-method standard error √(∇gᵀ V ∇g) using the EXACT gradient and the stored covariance V = vcov(fit);

• ci — the (lower, upper) Wald interval corrected ± z·se at level;

• level — the confidence level used.

Method

With θ̂ ≈ N(θ, V), a second-order Taylor expansion of g about θ̂ gives the expectation E[g(θ̂)] ≈ g(θ̂) + ½·tr(H_g(θ̂)·V), and the first-order delta method gives the variance ∇g(θ̂)ᵀ V ∇g(θ̂). The gradient ∇g and Hessian H_g are obtained by automatic differentiation of g, so they are exact for any smooth g.

Anchors / what to trust

• Linear g (identity anchor): for g(θ) = θ_k (or any affine map) H_g = 0, so corrected == estimate and se, ci equal the plug-in Wald values exactly.

• Curved g (curvature anchor): for g = exp on a coordinate with Wald mean m and variance v, corrected = exp(m)(1 + v/2), the second-order expansion of the analytic E[exp(θ_k)] = exp(m + v/2) (they agree to O(v²)).

Caveats

This is a second-order correction under the asymptotic premise θ̂ ≈ N(θ, V); it does not fix non-normality of θ̂ or a non-PD V at a variance boundary (there, prefer bootstrap_ci / profile confint). It is kept distinct from the raw plug-in accessors on purpose.

Example

fit = drm(bf(@formula(y ~ x), @formula(sigma ~ x)), Gaussian(); data = (; y, x))
k   = only(findall(==(:sigma), first.(fit.blocks)))   # σ block
ki  = fit.blocks[k].second[1]                          # log σ intercept index
bc  = bias_correct(fit, θ -> exp(θ[ki]))               # back-transformed σ
bc.estimate, bc.corrected, bc.se, bc.ci

source

bias_correct(θ̂::AbstractVector, V::AbstractMatrix, g::Function; level = 0.95)

Lower-level form operating directly on a point estimate θ̂ and its covariance V (e.g. for testing against closed forms, or for derived quantities not tied to a DrmFit). Same return shape as the DrmFit method.

source

Boundary inference

DRModels.chibar_pvalue Function

chibar_pvalue(stat::Real, q::Integer = 1) -> Float64

Chi-bar-square (χ̄²) boundary-corrected upper-tail p-value for a likelihood-ratio statistic stat = 2·(ℓ_full − ℓ_reduced) that tests q variance components = 0. Because a variance is constrained to be non-negative, testing it at zero is a boundary problem and the naive χ²(q) reference distribution is wrong (and conservative); the correct null is a mixture of χ² distributions (Self & Liang 1987; Stram & Lee 1994).

Supported q:

• q = 1 (one boundary parameter): null is 0.5·χ²(0) + 0.5·χ²(1), so

  p = 0.5 · P(χ²₁ > stat).

  At stat = 0 this returns 0.5; for stat > 0 it is exactly half the naive χ²(1) p-value.

• q = 2 (two independent boundary parameters): null is 0.25·χ²(0) + 0.5·χ²(1) + 0.25·χ²(2), so

  p = 0.5 · P(χ²₁ > stat) + 0.25 · P(χ²₂ > stat).

  At stat = 0 this returns 0.75.

Assumptions

• The q tested parameters are variances tested at the boundary 0.

• For q = 2, the two components are independent (uncorrelated information); otherwise the 0.25/0.5/0.25 weights are only approximate.

• All other parameters are interior (regularity away from the boundary).

• ML fits (as with lrtest; REML likelihoods are not cross-comparable).

A negative stat (the reduced model fit better — non-nesting or non-convergence) is clamped to 0, returning the boundary value (0.5 for q = 1, 0.75 for q = 2); inspect stat directly in that case. Only q ∈ (1, 2) are supported; other values throw an ArgumentError.

Example

stat = 3.5
chibar_pvalue(stat, 1)              # 0.5 * P(χ²₁ > 3.5)
chibar_pvalue(stat, 1) ≈ 0.5 * ccdf(Distributions.Chisq(1), stat)   # true
chibar_pvalue(0.0, 1) == 0.5        # boundary point mass
chibar_pvalue(0.0, 2) == 0.75

source

Cross-family post-fit

Accessors for a cross-family drm(...) result. The cross-family bivariate route is experimental and not ready for routine use: narrow documented evidence, no interval coverage, and the dependence it reports is a latent-scale scalar correlation (fit.rho_latent), not a rho12 formula. mf_coef is the tidy coefficient table; the other mf_* helpers live beside it in the module.

Engine constructors (q=4 and coevolution)

AugProblem, make_problem, and fit_q4_sparse_tmb are low-level exported interfaces behind the q=4 PLSM path. Use drm(...) for ordinary model fitting; these bindings serve scripts that prepare engine inputs directly.

AugProblem

AugProblem holds the augmented-phylogeny data and q=4 design matrices used by the sparse engine.

make_problem

make_problem(phy, y1, y2, X1, X2, Xs1, Xs2, Xr; species = 1:phy.n_leaves) builds that problem and its root-conditioned precision from a phylogeny.

fit_q4_sparse_tmb

fit_q4_sparse_tmb(prob, Q_cond; θ0 = ..., ...) runs sparse q=4 optimisation. Starting values must be supplied as either θ0 (the full parameter vector) or β0 (the mean coefficients).

The documented q=4 marginal evaluator and general-q coevolution bindings are listed below.

DRModels.marginal_nll Function

marginal_nll(prob, Q_cond, θ; u0, n_newton) -> (nll, û, ch_H, W)

Fresh inner Newton mode + the TRUE sparse Laplace NLL L(θ). u0 warm-starts the Newton (the optimiser reuses the previous mode). Computed in whitened coordinates v = (I ⊗ L⁻¹)u (Λ = LLᵀ; Λ⁻¹ is never formed), so it stays accurate as Λ approaches singularity. û is returned in the ORIGINAL u coordinates; ch_H is the CHOLMOD factor of the WHITENED H̃ = (I⊗Lᵀ)H_uu(I⊗L) (logdet H_uu = logdet H̃ − 2N Σ log Lᵢᵢ) and W the WhitenedPrior.

source

DRModels.CoevoProblem Type

CoevoProblem

A general-q coevolution dataset bound to a tree. Y[i, :] is the q-trait vector for data row i; X is the shared n × k fixed-effect design (same β columns per trait — the common comparative-biology setup); leaf_node[i] maps row i to its column in the root-conditioned augmented node set.

source

DRModels.lc_to_cov Function

lc_to_cov(v, q) -> Matrix

Map a length-q(q+1)/2 log-Cholesky vector to a q×q SPD covariance Λ = L L', diagonal of L exponentiated. Column-major lower-triangle order (j outer, i≥j inner) — identical to lc_to_Λ at q=4. Eltype follows v (AD-friendly).

source

DRModels.joint_missing_summary Function

joint_missing_summary(fit::PreparedJointFit)

Return a copy-safe named summary of prepared-row conditional moments and their data-only masks. This summary contains no std_error: its moments condition on fitted parameters. Use imputed for native-shaped uncertainty summaries.

source