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

DRModels.DrmFit Type

DrmFit

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

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 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. 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 and all non-Gaussian Laplace/GLMM phylo routes), 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

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

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 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 the structured (phylo/relmat/animal) routes. Returns an empty Dict for a fit with no random-effect block at all. A non-Gaussian GLMM fitted by the Gauss–Hermite/Laplace marginal routes (e.g. Binomial/Poisson with (1|g)) never computes conditional modes — those routes integrate the random effect out of the marginal likelihood without ever forming a posterior mode — so calling ranef on such a fit throws an ArgumentError instead of silently returning an empty Dict.

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

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

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

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

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())

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, 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

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 — 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 is exact when the relatedness matrix has a unit diagonal (for example an ultrametric tree of height 1). 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 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 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. This route is not available through the R bridge.

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

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 (DRModels.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

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 (DRModels.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

DRModels.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. A homogeneous Toeplitz temporal() fit returns the mean-coefficient rows only (drmTMB refuses the temporal bootstrap and defers σ and lag-correlation intervals).

source

DRModels.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. A homogeneous Toeplitz temporal() fit returns the mean-coefficient rows only, as bootstrap_ci.

source

DRModels.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. For a homogeneous Toeplitz temporal() fit, summary holds the mean-coefficient rows only, as bootstrap_ci.

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

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.

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

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 (#575): 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 (#477); 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 (#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

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

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

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

Prior (per-observation) weights used in the fit — drmTMB / glmmTMB's weights(). DRModels.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

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

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:

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

• 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" (still converged). Does not affect ok.

• 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

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

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 (estimates, SEs, z, p, CIs).

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

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

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, coords = nothing, newdata = nothing,
           options = Dict())

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 (#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

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 — the same primitive DRModels.profile_result / DRModels.bootstrap_result calls the R bridge previously had to reach by calling DRModels.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

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 (#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 DRModels.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 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 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

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

Why this exists: DRModels.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 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 input columns are returned unchanged.

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.

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

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(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

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

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.

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

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