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

DRM.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(y ~ x1, 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.

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(y ~ x + phylo(1 | sp), 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) on the Woodbury spine (#439), (c) Location–Scale–Scale (LSS) models (sd(g) ~ z, sd(species, phylogenetic) ~ z, and multi-component LSS models; #558), and (d) the bivariate q=4 PLSM Laplace engine (reml_q4).

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

Missing response handling

Incomplete responses (missing or NaN in y) are supported under the observed-rows pattern (matching response = "include" in the R bridge). For Location-Scale-Scale models (#559), the group index and scale design Z_g are parameterised over all G levels, while the likelihood is evaluated on observed rows.

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. Re-verified live 2026-08-25 against the installed drmTMB 0.7.0:

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, ...),
       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's first slice, method = :REML is not implemented for this family: the residual-only cell has no random effects to integrate out, and extending REML to the structured cells is a later slice; use method = :ML (the default).

d = (; y1 = exp.(0.4 .+ randn(200)), y2 = exp.(0.1 .+ randn(200)), x = randn(200))
fit = drm(bf(@formula(y1 ~ x), @formula(y2 ~ x),
             @formula(sigma1 ~ 1), @formula(sigma2 ~ 1), @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

DRM.DrmFit Type

DrmFit

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

source

Coefficients and (co)variance

DRM.fixef Function

fixef(fit) -> Vector{Pair}

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

source

DRM.re_sd Function

re_sd(fit) -> Dict{Symbol,Float64}

Estimated random-effect (random-intercept) standard deviations, keyed by grouping factor. A mean-axis random intercept (y ~ x + (1|g)) is keyed by the bare group name and is on the response scale. 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.

source

DRM.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

Wald 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 carry an interpretable Wald test only for blocks whose working-scale-zero null is a meaningful hypothesis: the location blocks :mu/:mu1/:mu2 (coefficient = 0) and :rho12 (atanh ρ12 = 0 ⇔ no correlation). For :sigma/:sigma1/:sigma2/:resd/:resid/:recov/ :phylocov the zero-on-working-scale null is not the scientific one — e.g. log σ = 0 means σ = 1, not the σ = 0 variance boundary — so those rows show z and Pr(>|z|) as NaN rather than a misleading test of an arbitrary scale reference (issue #320). To test a variance component against 0 use a boundary-corrected likelihood-ratio test (lrt_boundary); the estimate and SE for those blocks are still reported. A boundary / singular direction (Inf SE) also reports NaN z / p rather than a spurious z = 0, p = 1 (issue #323.2).

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

DRM.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 population-level scale coefficients (on the working/log scale) use coef(fit, :sigma) instead.

source

DRM.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

DRM.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

DRM.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

DRM.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

DRM.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

DRM.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.

Currently populated for the Gaussian closed-form RE paths (exact GLS conditional means). Returns an empty Dict for models without random effects. Non-Gaussian GLMM posterior modes (GHQ/Laplace) are not yet wired — see issue #73.

source

Prediction and simulation

StatsAPI.predict Function

predict(fit, newdata; type = :response, se = false) -> Vector / Dict / NamedTuple

Population-level prediction on newdata (a NamedTuple / column table), random / structured effects integrated out. type = :response (default) returns the response-scale mean — the family inverse link applied to Xβ̂ (exp for Poisson/Gamma, logistic for Beta/Binomial, identity for Gaussian); type = :link returns Xβ̂. In-sample, predict(fit, data) ≈ fitted(fit). Univariate returns a vector; bivariate returns Dict(:mu1 => …, :mu2 => …).

se = false (default) is the point prediction above. se = true adds delta-method standard errors for the mean (glmmTMB/drmTMB se.fit parity):

• univariate → a NamedTuple (; prediction, se);

• bivariate → a NamedTuple (; prediction::Dict, se::Dict) keyed :mu1, :mu2.

The SE uses the μ-block of vcov(fit): link scale se_i = sqrt(xᵢ' Vμ xᵢ) (with Vμ = vcov(fit)[r, r], r the :mu coef range from fit.blocks); response scale multiplies by the inverse-link derivative |dμ/dη| at η̂ (identity → 1, exp → exp(η), logistic → μ(1−μ)).

source

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

DRM.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), random / structured effects integrated out (exactly like predict). 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).

type = :response (default) applies each parameter's inverse link, so in-sample it reproduces marginal_parameters (i.e. fit.means[:mu], fit.scales[...]). 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

DRM.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).

In-sample these equal predict_parameters(fit, data) (response scale).

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

DRM.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

DRM.simulate Function

simulate(fit; nsim = 1, rng = default_rng())

Draw parametric (residual-level) replicate response(s) from the fitted model — the building block of a parametric bootstrap and posterior-predictive checks. Each draw uses the fitted per-observation mean μ̂ and the fitted dispersion / scale parameters for the family; for random-effect models the draw is conditional on the random effects being zero (population level).

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, 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
Y   = simulate(fit; nsim = 100)   # nobs × 100 Matrix

source

StatsAPI.fitted Function

fitted(fit)

Fitted mean(s). Univariate / random-effect models return the μ vector (the population/marginal mean Xβ̂); bivariate models return Dict(:mu1=>…, :mu2=>…).

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 — randomized quantile residuals (Dunn & Smyth; DHARMa / glmmTMB style). For observation i with fitted distribution F_i, r_i = Φ⁻¹(u_i) where u_i is the (randomized, for discrete families) probability-integral transform of y_i. Under a correct model the r_i are i.i.d. standard normal. Univariate only.

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

• 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, NegBinomial2, TruncatedNegBinomial2, Binomial, BetaBinomial, CumulativeLogit (ordinal);

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

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

source

Predicting distributional parameters

Prediction interpretation

The embedded predict_parameters docstring above uses legacy wording about integrating effects out. The current implementation sets random and structured effects to zero. With a nonlinear link, these are different predictions: for example, exp(η) differs from averaging exp(η + b) over a non-degenerate random effect b. Use the fixed-effect interpretation here.

predict returns the response (mean) prediction, but a distributional regression also models the scale and—bivariately—the 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, with random / structured effects set to zero. This is an inverse-link transformation of the fixed-effect linear predictor, not response-scale integration 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 DRM, 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; fenced for v1.0 (D-181). API and numerics may change; not covered by the R-parity scoreboard.

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. See the engine-internals reference for the prepared-state design and route limits.

DRM.CategoricalLogit Type

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

source

DRM.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

DRM.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

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 (DRM.jl#631). 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. Canonical location–scale profiling remains serial in this release, including when threads = true.

• 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.

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

DRM.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 (DRM.jl#631);

• 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.

source

DRM.bootstrap_ci Function

bootstrap_ci(formula, family; data, B = 300, level = 0.95, rng = default_rng(), threads = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)
bootstrap_ci(fit; data, B = 300, level = 0.95, rng = default_rng(), threads = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)

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 / 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. If you already have fit = drm(...), pass the fit directly to avoid refitting the base model before the bootstrap replicates.

source

DRM.bootstrap_summary Function

bootstrap_summary(formula, family; data, B = 300, level = 0.95, rng = default_rng(), threads = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)
bootstrap_summary(fit; data, B = 300, level = 0.95, rng = default_rng(), threads = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)

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, any failed replicate errors after all failures are recorded. Set failures = :skip to compute summaries from successful replicates; call bootstrap_result to inspect the skipped failures.

source

DRM.bootstrap_result Function

bootstrap_result(formula, family; data, B = 300, level = 0.95, rng = default_rng(), threads = false, failures = :error, check_converged = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)
bootstrap_result(fit; data, B = 300, level = 0.95, rng = default_rng(), threads = false, failures = :error, check_converged = false, K =, A =, tree =, coords =, algorithm = :auto, g_tol = 1e-8)

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;

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

• 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 failures and then errors if any replicate failed. failures = :skip computes summaries from successful replicates and keeps the failure records in the return value. Set check_converged = true to treat non-converged refits as failed replicates. Passing an existing DrmFit reuses that point estimate as the bootstrap seed fit and starts directly with the B simulated refits. Gaussian bootstrap refits pass algorithm and g_tol through to drm(...); this is useful for large structured models where :auto selects a sparse route and the tolerance is part of the benchmarked workflow.

source

Information criteria

DRM.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

DRM.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

DRM.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.

A convention gap on the bivariate q=2/q=4 routes (#477)

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) now report the same normalised scale (#477, 2026-08-25).

They previously omitted the (n_β/2)·log(2π) constant while these univariate routes included it, so reml_loglik(fit) meant different things depending on which route produced the fit. For the q=4 phylo layout with n_β = 6 the gap was 3·log(2π) ≈ 5.51 — large enough to read as a real disagreement between engines rather than a labelling difference, which is exactly how it misled this project once (see the corrected note in test/parity/q4-reml/biv-q4-phylo-reml/expected.toml).

Every REML route in DRM.jl now reports the normalised form, matching lme4, glmmTMB, TMB and drmTMB. See fit_q4_reml's docstring in src/reml_q4.jl for the derivation and for the evidence: the q=4 parity gate's atol_loglik fell from 5.5436 to 0.03 once the constant was no longer being absorbed.

source

DRM.estimation_method Function

estimation_method(fit) -> Symbol

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

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

DRM.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 (#136).

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

DRM.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 (DRM.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. 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 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), the tested null (variance = 2. 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 this case and emits a one-time warning; use lrt_boundary (or a parametric bootstrap) for a boundary-correct 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

DRM.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) -> Vector{Float64}

Prior (per-observation) weights used in the fit — drmTMB / glmmTMB's weights(). DRM.jl fits do not currently store prior weights, so this returns ones(nobs(fit)) (every observation weighted equally). Extends StatsAPI.weights.

Example

fit = drm(bf(@formula(y ~ 1 + x), @formula(sigma ~ 1)), Gaussian(); data)
weights(fit) == ones(nobs(fit))   # all-ones prior weights

source

DRM.update Function

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

Refit fit's model with a new formula (a bf bundle), reusing the fitted family — the convenience refit verb, mirroring R's update. Equivalent to drm(formula, family(fit); data = data, kwargs...).

data must be supplied: a DrmFit does not retain its data, so update cannot reuse the original observations. Any extra keyword arguments (K, A, tree, coords, g_tol, …) are forwarded to drm.

Example

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

source

Diagnostics and accessors

DRM.check_drm Function

check_drm(fit) -> NamedTuple

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

• converged — the optimiser's convergence flag.

• 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 (penalty = drm_phylo_penalty(...)). Such a fit reports standard errors from the penalized curvature, which are credible-interval-shaped rather than frequentist, and loglik is the unpenalized data log-likelihood.

• ok — true when converged, the gradient is small, and the covariance is PD. On a penalized fit the gradient criterion is dropped: the stored objective is unpenalized, so its gradient is non-zero at the MAP optimum by construction and scoring it would report a correct fit as broken. It 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

DRM.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

DRM.is_converged Function

is_converged(fit::DrmFit) -> Bool

Whether this fit may be trusted: the optimiser reported convergence and the optimum is not degenerate. A false here means the reported estimates / standard errors should not be trusted.

This is deliberately STRICTER than the raw fit.converged flag. 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 to +Inf. Optim.converged only asks whether the gradient test was met, and at such a point it returns true.

Measured 2026-08-24 on a one-row-per-species phylo fit (#461): sd_phylo = 22980, sigma = 7.5e-15, loglik = 6.8e13, converged = true — and 25% of parametric-bootstrap replicates landed on such a point. Downstream that is WORSE than an outright failure, because every consumer treats the fit as usable and a percentile interval silently inherits the nonsense.

Checked here, at the single public accessor, rather than at the ~30 DrmFit construction sites across 20 family files — one place that every consumer already goes through. fit.converged still exposes the raw optimiser flag for anyone who wants it.

source

StatsAPI.deviance Function

deviance(fit::DrmFit) -> Float64

Deviance of the fitted model, -2 · loglik(fit) — drmTMB's deviance(). 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

Fit summaries and route inventories

Base.summary Method

summary(fit::DrmFit)

Coefficient table for a fitted model — the DRM.jl analogue of drmTMB's summary(). Returns the same CoefTable as coeftable (estimates, SEs, z, p, CIs).

source

DRM.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 DRM.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 any approximated or borrowed number (#466).

Coverage by family

Wired (reports Optim.iterations(res) from the LBFGS run that produced θ̂): 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) are wired 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.

Still -1 (no single outer LBFGS call to attribute the count to, or not yet wired): 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.

source

DRM.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, :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

DRM.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), 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

DRM.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

DRM.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

DRM.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).

• 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

DRM.repeatability Function

repeatability(fit; component = nothing, level = 0.95, method = :delta) -> 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.

source

DRM.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.

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. The R ↔ Julia bridge defines its admitted cells and refusals; these docstrings do not expand that contract.

DRM.drm_bridge Function

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

Fit a DRM.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 (#492). These materialised columns are NOT (yet) reconstructed for newdata — a formula using them combined with newdata fails loudly (missing column) rather than silently mismodelling; see #467.

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.

source

DRM.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 — the same primitive DRM.profile_result / DRM.bootstrap_result calls the R bridge previously had to reach by calling DRM.jl's underscore-prefixed marshalling internals directly (see #475); this kwarg is the supported route that replaces that qualified-internal call. 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

DRM.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 (#575) reached through the SAME marshalling-friendly boundary drm_bridge uses — one SUPPORTED entry point for the drmTMB R shim (drm_julia_reml_objective_at(), R/julia-bridge.R), replacing its previous dependency on five private DRM.jl names (_bridge_data, _bridge_formula, _bivariate_q4_marker, _design, _phylo_species_index) reached by qualified name. 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 DRM.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 DRM.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 is a DIAGNOSTIC: it selects nothing, fits nothing, and promotes no capability-ledger row. See #575 for the cross-engine mode-finder-vs- objective-translation question it exists to answer.

source

DRM.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(y ~ x, sigma ~ x), data)        # @warn names the dropped rows
fit   = drm(bf(y ~ x, sigma ~ x), Gaussian(); data = clean)

Why this exists: DRM.jl has no native missing-data path, so a raw missing/NaN response or predictor makes drm error with an opaque message. drm_listwise gives a guided, documented complete-case route and emits a clear @warn reporting how many rows (and which columns) were removed.

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 tracked as follow-up under issue #49 (report/fiml-missing-data-design.md); this helper is the 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 input columns are returned unchanged.

source

Staged pair association diagnostics

DRM.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

DRM.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.

DRM.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. DRM.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.

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(mu = @formula(y ~ x + phylo(1 | species)), sigma = @formula(~1)),
          Gaussian(); data = dat, tree = tree, penalty = pen)
fit.phylo_penalty        # the penalty at the optimum
loglik(fit)              # UNPENALIZED data log-likelihood

source

DRM.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(mu = @formula(y ~ x + phylo(1 | species)),
           sigma = @formula(~ 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

DRM.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

DRM.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.

DRM.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 class: gaussian × binomial (gaussian_bernoulli). The other four reviewed classes need a rectangle probability and are not implemented yet; they error rather than silently approximating.

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

DRM.latent_normal Function

latent_normal()

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

source

DRM.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

DRM.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.

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];

• level — the confidence level;

• method — the method used.

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

source

DRM.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].

source

DRM.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

DRM.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 fit_mixed_family result. The cross-family bivariate route is experimental and not release-ready: single-fixture 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.

DRM.marginal_nll Function

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

Fresh inner Newton mode + the TRUE sparse Laplace NLL L(θ). u0 warm-starts the Newton (the optimiser reuses the previous mode). Returns the CHOLMOD factor ch_H and the sparse prior P so callers can reuse them for the gradient.

source

DRM.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

DRM.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