Engine internals and advanced constructors

Developer reference

This appendix records docstrings for the sparse, coevolution, bridge, and marginal-likelihood implementation surface. Names beginning with _ are nonpublic implementation helpers. Exported advanced constructors are included for source coverage, not as a claim that every combination is a supported drm(...) route. Start with the model guides and capability matrix for ordinary fitting.

Marginal-method and association internals

DRM.AGHQ Type

AGHQ <: MarginalMethod

1-D Liu–Pierce adaptive Gauss–Hermite quadrature (#448). Opt-in via marginal = :AGHQ on Poisson (1 | g) only (nAGQ=5 default). This is plumbing — k=1 ≡ 1-point Laplace; k≈5 nll agrees with GHQ-32. Not a recovery headline, not a capability-chip flip, not tensor AGHQ on phylo Laplace. :REML is not wired to :AGHQ this slice.

source

DRM.Laplace Type

Laplace <: MarginalMethod

Laplace marginal: Gaussian approximation at the posterior mode. The default. On Poisson (1 | g) the public :LA path is non-adaptive GHQ-32, not 1-point Laplace and not AGHQ.

source

DRM.Variational Type

Variational <: MarginalMethod

Gaussian-variational (VA/ELBO) marginal — opt-in alternative to Laplace for bias-sensitive random-effect models (#136). Public drm path (Experimental): Poisson / Binomial / NegBinomial2 / Gamma / Beta random intercept (1 | g) via marginal = :VA (scale families need sigma ~ 1). Phylo, crossed, correlated slopes, ZI/hu, and 136e stay unwired; #136 stays open.

source

DRM.MarginalMethod Type

MarginalMethod

How a model's random-effect integral is approximated. Subtypes: Laplace (mode + curvature; the default, what drmTMB/TMB use), Variational (maximize an ELBO over a Gaussian q; opt-in, steadier on dispersion/shape — #136), and AGHQ (1-D Liu–Pierce adaptive Gauss–Hermite; opt-in Poisson (1 | g) only — #448).

source

DRM.LatentNormal Type

LatentNormal

Association kernel for associate_pairs — drmTMB's latent_normal(). The two margins are coupled through a bivariate standard normal latent with correlation eta.

source

DRM._aghq_1d_logint Function

_aghq_1d_logint(logf, mode, hess; k=5)
_aghq_1d_logint(logf, mode, hess, z, w)

Log of ∫ exp(logf(x)) dx by Liu–Pierce (1994) adaptive Gauss–Hermite quadrature, wrapping _gauss_hermite.

mode maximises logf; hess = logf''(mode) must be negative. Nodes and weights (z, w) are the physicists' rule ∫ h(x) e^{-x²} dx ≈ Σ wₖ h(zₖ). The adaptive map is x = mode + √2 σ z with σ = 1/√(−hess), so

∫ exp(logf) ≈ √2 σ Σₖ wₖ exp(logf(xₖ) + zₖ²).

k=1 (node 0, weight √π) is exactly the 1-point Laplace approximation logf(mode) + ½ log(2π) − ½ log(−hess). That is a plumbing identity, not a quadrature or recovery claim.

A non-scalar mode fails loud (dim ≠ 1). Multi-d / tensor AGHQ is out of scope for #448.

source

DRM._aghq_require_1d Function

_aghq_require_1d(dim)

Throw ArgumentError unless dim == 1. 1-D Liu–Pierce AGHQ only (#448).

source

DRM._assoc_logdiffexp Function

_assoc_logdiffexp(a, b)

log(exp(a) − exp(b)) for a ≥ b, computed without forming either exponential. Returns -Inf when the two are indistinguishable rather than a negative argument to log.

source

DRM._poisson_group_aghq_logint Function

_poisson_group_aghq_logint(y, η0, lf, idx, σb, k)

Log ∫ L(y_g | b) φ(b; 0, σb²) db for one Poisson (1 | g) group via 1-D Liu–Pierce AGHQ.

The group posterior is strictly concave. The inner mode is a fixed 12-step Newton unroll (same discipline as the VA inner solve) so the outer nll stays ForwardDiff-smooth. k=1 recovers 1-point Laplace of this integrand — plumbing, not a recovery claim.

source

Profile nuisance-status internals (no stability guarantee)

DRM._ProfileNuisanceResult Type

Internal result of one fixed-coordinate nuisance solve.

accepted is deliberately stronger than a finite optimizer-reported minimum: the optimizer must have terminated successfully, its minimizer must be finite, and the objective is evaluated again at that minimizer. We do not impose a separate score threshold here; Optim termination is not a proof of stationarity.

source

Augmented phylogeny and sparse linear algebra

DRM.AugmentedPhy Type

AugmentedPhy{T}

Augmented-state sparse phylogenetic precision for a rooted tree, including multifurcations.

Fields

• n_leaves::Int – number of tip species (p).

• n_total::Int – leaves + internal ancestor nodes.

• Q_topology::SparseMatrixCSC – (n_total × n_total) topology contribution to the sparse precision. The actual phylogenetic precision is Q_topology / σ²_phy. A positive-length tree stores 3n_total - 2 nonzeros.

• leaf_indices::Vector{Int} – maps a leaf k ∈ 1:p to its row/col in the augmented state. Ordering matches the order leaves were encountered in the Newick string (left-to-right).

• leaf_names::Vector{String} – species names parsed from the Newick.

• branch_lengths::Vector{T} – the n_total - 1 branch lengths in the order the parser walked the tree.

• root_index::Int – which augmented row is the root.

Q_topology is positive semi-definite (rank n_total - 1). The all-ones vector is its sole zero eigenvector — fixing the root removes the degeneracy. The sparse log-likelihood path adds a positive contribution to the leaf diagonals (proportional to λ_phy² / d_total) which renders the active solve matrix positive definite without any explicit ridge.

source

DRM.Vinv_mul Function

V^{-1} z via Woodbury: (z - S (M^{-1} (1/σ²) S' z)) / σ²

source

DRM.build_M Function

Build M = P + (1/σ²) S'S. Return (P, M, chM, chP). chP = chol(P + ridge) for logdet P; chM = chol(M) for the Woodbury solves.

source

DRM.exact_traces Function

EXACT Tr(Q_cond M^{-1}) and Tr(S M^{-1} S') via Takahashi selected inverse. Returns (tr_QM, tr_SMS).

source

DRM.logdetV_val Function

logdet V = n log σ² + logdet M − logdet P

source

DRM.make_phy Function

make_phy(edges::AbstractVector{<:Tuple}, n_leaves::Integer;
         root_index::Integer = -1) :: AugmentedPhy{Float64}

Convenience constructor: build an AugmentedPhy from a list of edges (parent_id, child_id, branch_length) with contiguous integer node ids 1:n_total. Leaves must be exactly 1:n_leaves; internal nodes follow them. The root may have any internal id. Every internal node must have at least two children, and every non-root node must have one parent.

If root_index < 0 it is auto-detected as the unique node that is not a child in any edge.

This bypasses the Newick parser — useful for tests and for trees that arrive from another tool already as edge lists.

source

DRM.make_problem_from_Q Function

make_problem_from_Q(Q_cond, y1, y2, X1, X2, Xs1, Xs2, Xr; group) -> (prob, Q)

Build an AugProblem from a known G×G group-level precision (relmat / animal / fixed-range spatial) rather than from an augmented phylogeny. group[i] maps data row i to a 1-based level in 1:G. Returns (prob, sparse(Q_cond)) for fit_q4_sparse_tmb.

This is the #189 level-indexed entry point; the phylo path keeps make_problem.

source

DRM.q4_marginal_diagnostic Function

q4_marginal_diagnostic(prob, Q_cond, theta; u0=nothing, n_newton=40,
                       gradient=false)

Developer diagnostic for the q=4 ML Laplace objective. It evaluates the same pieces used by marginal_nll in order and returns a NamedTuple with ok, first_nonfinite, stages, and (when finite) nll.

This is intentionally diagnostic-only: the optimizer keeps its usual Inf barrier, while hard large-data failures can call this helper to see whether the first non-finite component is the parameter vector, among-axis covariance, sparse prior precision, inner mode, Laplace components, or exact-gradient assembly.

source

DRM.takahashi_diag Function

takahashi_diag(ch::SparseArrays.CHOLMOD.Factor) -> Vector{Float64}

Convenience: return ONLY diag(Q⁻¹) (a length-n vector, in the ORIGINAL ordering) via the Takahashi recursion. Same cost as takahashi_selinv but without materialising the full sparse output (a small allocation win when only the diagonal is needed, as in the EM E-step's per-trait variance).

source

DRM.takahashi_selinv Function

takahashi_selinv(ch::SparseArrays.CHOLMOD.Factor) -> SparseMatrixCSC

Compute the Takahashi selected inverse of the matrix Q whose sparse Cholesky factor is ch (P · Q · Pᵀ = L · Lᵀ with P = I[ch.p, :]). Returns a SparseMatrixCSC holding Q⁻¹ (in the ORIGINAL un-permuted ordering) at the union sparsity of Pᵀ (L + Lᵀ) P. Entries outside that pattern are NOT computed (and are NOT zero in general).

Cost: O(nnz(L)) arithmetic + O(nnz(L)·log(max_col_nnz)) for the symmetric lookups (with constant max_col_nnz on a tree this is O(nnz(L)) overall).

source

DRM.theta_len Function

Total θ length = sum(beta widths) + 10 (log-Cholesky Λ).

source

DRM.unpack_theta Function

Slice θ into a β NamedTuple and the 10-vector lc (eltype follows θ).

source

DRM.pack_theta Function

Build a Float64 θ from a β NamedTuple and a Λ matrix.

source

Newick and topology parser internals (no stability guarantee)

DRM._phy_validate_topology Function

Validate a rooted acyclic tree and return (root, children) without Q assembly.

source

DRM._parse_label! Function

Parse one Newick label, retaining quoted Unicode/control characters exactly.

source

q=4 and D·R·D separation evaluators

The fz_ prefix and fit_q4_sparse_fisherz are historical names: the correlation map is spherical/LKJ, not Fisher-z. The warning on DRM.fit_q4_sparse_fisherz below gives the correct back-transform.

DRM.beta_widths Function

Column widths (k1,k2,ks1,ks2,kr) of the five design matrices.

source

DRM.fit_q4_reml Function

fit_q4_reml(prob, Q_cond; beta0, Lambda0, [phi0], g_tol, ...) -> NamedTuple

Fit REML objective over phi = (beta_rho, lc). beta_mu AND beta_sigma (the location and scale fixed effects) are profiled out internally; only beta_rho stays outer.

Returns NamedTuple: (phi, beta, Lambda, reml_loglik, ml_loglik, converged, iterations, g_residual, f_calls, u_hat)

Automatic warm restart (#484)

On some cells the REML LBFGS's first line-search step from the ML warm start fails outright (zero accepted steps — a starting-value problem, not slow convergence, so a bigger iterations or looser g_tol cannot fix it). That exact stall (!converged with the minimizer still sitting at phi0) is detected automatically and retried by re-deriving a coarser ML warm start and continuing from there — judged at the SAME g_tol the caller passed, never a loosened one. A REML fit whose first attempt already moves at all — converged or not — is completely unaffected; see the block above phi_hat for the mechanism. Needs beta0/Lambda0 to fire (skipped if the caller supplied phi0 directly, since there is then nothing to re-derive a coarser start from).

Normalisation convention (#477)

reml_loglik reports the normalised Patterson–Thompson restricted log-likelihood: the raw objective ℓ_ML(θ, β̂) − ½ logdet(S) plus (n_β/2)·log(2π), the constant from integrating the flat prior over the n_β marginalised fixed effects (n_β = the combined width of the beta_mu1, beta_mu2, beta_s1, beta_s2 designs — exactly the Schur complement's dimension; beta_rho is never marginalised, so it does not count). That matches lme4, glmmTMB and TMB, so reml_loglik is directly comparable across engines.

Changed 2026-08-25 (#477). It previously reported the unnormalised form, while DRM.jl's own univariate REML routes — _fit_fixed_gaussian_reml (gaussian_core.jl), the Gaussian mean (1 | g) route (gaussian_ranef.jl) and location_only.jl — already added the constant. So one package reported two different scales under one name, and reml_loglik(fit) meant different things depending on which route produced the fit. That was an inconsistency rather than a convention choice: the convention had already been made on the univariate side and the bivariate routes simply had not followed it.

The evidence is the q=4 parity gate. Its atol_loglik was 5.5436, of which 5.513631 was this constant — a tolerance that existed almost entirely to absorb the offset, and therefore tested almost nothing. It is now 0.03, the cross-optimum spread alone: a 185× tightening, verified 33/33. A constant cannot move the argmax, so the optimisation is untouched; only the reported value moved.

reml_q2.jl carries the same change and the same derivation, but has no parity fixture of its own — it is verified only by sharing this one's arithmetic.

source

DRM.fit_q4_sparse_fisherz Function

fit_q4_sparse_fisherz(prob, Q_cond; β0, Σa0 / Λ0, ...) -> NamedTuple

Fit the q=4 PLSM in the separation (D·R·D) OUTER parameterization, with Σ_a = D R D, D = diag(exp.(d)) and R = C Cᵀ a correlation matrix. The inner engine is the UNTOUCHED marginal_and_exact_grad. Returns Σ_a and the 6 among-axis correlations directly, alongside the usual fit diagnostics.

The name says Fisher-z; the map is not Fisher-z

fisherz here, and the fz_ prefix throughout, is historical. R is built from a spherical/LKJ correlation-Cholesky whose angles are α = π·(tanh θ + 1)/2 ∈ (0, π), so a fitted angle-real θ back-transforms as ρ = cos(π·(tanh θ + 1)/2) — not the Fisher-z bijection ρ = tanh θ. Inverting with tanh returns the wrong correlation. Only the PD/interior property carries over from the verified q=2 Fisher-z path, not the link itself. The 6 among-axis correlations are returned already on the ρ scale, so prefer them to inverting θ_R by hand.

source

DRM.fz_DRD Function

Σ_a = D R D (4×4 PD) from φ_a = [d (4 log-SDs); θ_R (6 angle-reals)].

source

DRM.fz_R Function

4×4 correlation matrix R = C Cᵀ from the 6 angle-reals θ_R (always PD, unit diag).

source

DRM.fz_R_chol Function

Lower-tri correlation-Cholesky C (4×4) from the 6 angle-reals θ_R (R = C Cᵀ).

source

DRM.fz_init_from_Sigma Function

Initial φ_a = [d; θ_R] reproducing a starting Σ_a (used to seed the optimiser).

source

DRM.fz_marginal_and_grad Function

fz_marginal_and_grad(prob, Q_cond, ψ; u0, n_newton) -> (nll, g_ψ, û, chH)

TRUE sparse Laplace NLL and its EXACT gradient in the D·R·D OUTER parameters ψ = [β; φ_a]. Calls the UNTOUCHED marginal_and_exact_grad on θ = [β; lc(φ_a)], then chain-rules the Σ_a block: g_φ = Jᵀ g_lc with J = ∂lc/∂φ_a. β passes through.

source

DRM.fz_marginal_nll Function

D·R·D outer marginal NLL only (for FD verification / line search).

source

DRM.fz_psi_len Function

Length of ψ = nβ + 10 (the D·R·D outer parameter vector).

source

DRM.fz_psi_to_theta Function

ψ (D·R·D outer) → engine θ (β + lc). Pure reparameterization of the Σ_a block.

source

DRM.fz_unpack_psi Function

Split ψ into the β-block (Float-passthrough) and φ_a (10).

source

DRM.marginal_and_exact_grad Function

marginal_and_exact_grad(prob, Q_cond, θ; u0, n_newton) -> (nll, grad, û, ch_H)

TRUE sparse Laplace NLL and its EXACT 17-dim gradient (cheap + implicit), with the inner mode û FROZEN. All CHOLMOD-blocked logdet θ-derivatives use the Takahashi selected inverse of H (O(p)); all other pieces are single-level AD. Returns û and the factor so the caller can warm-start the next evaluation.

source

Coevolution and phylogenetic interaction kernels

DRM.coevo_marginal_cov Function

coevo_marginal(prob, Q_cond, β, Λ, σ_res) -> (ℓ, û, ch_H, P)

EXACT Laplace (= Gaussian) marginal log-likelihood at the given parameters, plus the inner mode û, the CHOLMOD factor of H_uu, and the sparse prior P. β is k × q (trait-major columns), Λ is q×q SPD, σ_res a length-q vector of residual SDs.

source

DRM.coevo_pack Function

Pack (β::k×q, Λ::q×q, σ_res::q) → θ (Float64).

source

DRM.coevo_unpack Function

Unpack θ → (β::k×q, Λ::q×q, σ_res::q). Eltype follows θ.

source

DRM.cov_to_lc Function

cov_to_lc(Λ) -> Vector{Float64}

Inverse of lc_to_cov: log-Cholesky vector of a SPD Λ (q inferred from size). Column-major lower-triangle order.

source

DRM.fit_coevolution Function

fit_coevolution(prob, Q_cond; β0, Λ0, σ0, g_tol, iterations) -> NamedTuple

Fit the general-q coevolution model by maximising the EXACT (conjugate-Gaussian) Laplace marginal over (β, log-Cholesky Λ, log σ_res). Uses LBFGS with a central finite-difference gradient: the marginal is a single sparse solve, so FD is cheap, and a quasi-Newton method converges far faster than a simplex on the k·q + q(q+1)/2 + q-dimensional θ. ForwardDiff cannot flow through the CHOLMOD factor (the same constraint the q=4 engine hits), hence FD rather than AD here. Returns (; β, Λ, σ_res, loglik, converged, iterations, θ).

source

DRM.fit_coevolution_q2_reml Function

fit_coevolution_q2_reml(prob, Q_cond; β0, Λ0, σ0, rho0, g_tol, iterations) -> NamedTuple

REML (Patterson–Thompson) fit for the bivariate q=2 structured residual- correlation coevolution model — the REML counterpart of fit_coevolution_q2_residual. Marginalises beta_mu1 and beta_mu2 (the two axes carrying the structured random effect Λ); beta_sigma1, beta_sigma2, and beta_rho12 stay outer parameters, since sigma1/sigma2/ rho12 are intercept-only fixed effects with no random-effect axis to integrate against. See this file's header for why exactly those axes (and not, e.g., all of them as in the q=4 engine).

Because mu1/mu2 enter the leaf residual linearly and (Λ, D) are held fixed while profiling beta, the profile step is an EXACT single Newton step (not an iterative/alternating scheme) — see _q2_profile_and_schur.

Returns (; β, Λ, residual_cov, σ_res, rho12, reml_loglik, ml_loglik, loglik, converged, iterations, u_hat), mirroring fit_coevolution_q2_residual's field names (loglik is set to reml_loglik, matching the q=4 REML fitter's convention in _fit_bivariate_q4_structured).

Normalisation convention (#477)

Like reml_q4.jl's fit_q4_reml, reml_loglik here reports the normalised Patterson–Thompson restricted log-likelihood: the raw objective ml_ll - 0.5*logdet(S) plus the (n_β/2)·log(2π) that lme4/glmmTMB/TMB add when integrating the flat prior over the marginalised fixed effects. Here n_β = length(β̂) (beta_mu1 + beta_mu2; beta_sigma1/beta_sigma2/ beta_rho12 stay outer and are never marginalised, and rho12 comes from D̂ rather than from β̂, so it is correctly not counted).

Changed 2026-08-25 (#477), together with the q=4 route — see fit_q4_reml's docstring for the derivation and the evidence. Note this route has no parity fixture of its own: it is verified only by sharing the q=4 route's arithmetic, which the q=4 gate confirmed by tightening from 5.5436 to 0.03. A q=2 REML parity fixture would be the way to check it directly.

source

DRM.fit_coevolution_q2_residual Function

fit_coevolution_q2_residual(prob, Q_cond; ...) -> NamedTuple

Fit the q=2 phylogenetic coevolution model with a bivariate residual covariance. This is the same target as drmTMB's bivariate Gaussian q2 location route: two phylogenetic random-effect axes (mu1, mu2) plus residual rho12. It is ML-only and complete-response only at the caller level.

source

DRM.fit_phylo_interaction Function

fit_phylo_interaction(y, X, C_A, C_B; pμ = size(X, 2),
                      munames = ["(Intercept)"], group = :interaction,
                      g_tol = 1e-8) -> DrmFit

Fit the Gaussian bipartite two-tree interaction model by maximum likelihood. C_A (n_A×n_A) and C_B (n_B×n_B) are the host/parasite phylo correlation matrices; observations are in Kronecker order (cell (a, b) at row (a−1)·n_B + b), so length(y) == n_A·n_B. Returns a DrmFit whose loglik is the marginal log-likelihood; the interaction SD σ and the residual SD σ_e are recovered via re_sd(fit)[group]andexp(coef(fit, :sigma)[1]).

A dense n×n covariance is assembled and factored each evaluation — O(n³), intended for modest n.

source

DRM.lc_len Function

Number of free log-Cholesky parameters for a q×q covariance: q(q+1)/2.

source

DRM.lc_metric Function

lc_metric(prob, Q_cond, θ, u0; h=1e-5) -> Matrix{Float64}

Observed-information (Fisher-scoring) metric on the 10 log-Cholesky parameters of the q=4 among-axis covariance Λ.

Builds a 10×10 Hessian of the marginal NLL w.r.t. θ[8:17] by central finite-differences of the exact lc-gradient from marginal_and_exact_grad, symmetrises, and ridge-projects to SPD so the metric is usable as a preconditioner / natural-gradient operator (AI-REML / Fisher scoring).

Notes

• Infrastructure for REML / non-Gaussian RE gradient work (#11 / #165) — not a public ML solver. The natural-gradient EM that used this metric failed the #13 MLE-parity gate against fit_q4_sparse_tmb (see docs/dev-log/plans/2026-08-01-natgrad-decision-gate.md).

• Cost: 20 exact-gradient evaluations per call (warm-started from u0).

source

DRM.make_coevo_problem Function

make_coevo_problem(phy, Y, X; species) -> (prob::CoevoProblem, Q_cond)

Build a CoevoProblem and the root-conditioned sparse tree precision Q_cond (q-agnostic; the SAME object the q=4 engine consumes). species[i] maps data row i to a tip 1:phy.n_leaves; default one row per tip.

source

DRM.make_coevo_problem_from_covariance Function

make_coevo_problem_from_covariance(K, Y, X; group) -> (prob, Q)

Build a CoevoProblem from a known structured covariance/relatedness matrix over observed levels. The returned precision is K^-1, kept sparse for the shared exact-Gaussian coevolution fitter.

source

DRM.make_coevo_problem_from_precision Function

make_coevo_problem_from_precision(Q, Y, X; group) -> (prob, Q)

Build a CoevoProblem from a known structured precision matrix over observed levels rather than from an augmented phylogeny. group[i] maps row i to a 1-based level in Q. This is the direct exact-Gaussian q2/q route for known relatedness/covariance fixtures such as relmat and animal matrices.

source

DRM.phylo_interaction_nll Function

phylo_interaction_nll(θ, y, X, CAB; pμ) -> Real

Marginal negative log-likelihood of the bipartite two-tree interaction model at parameters θ = [βμ(pμ); log σ_e; log σ], with CAB = C_A ⊗ C_B (the dense n×n interaction correlation) precomputed. The marginal is exactly Gaussian: V = σ²·CAB + σ_e²·I. A non-PD step returns a large finite penalty (so a line-search probe never throws / returns Inf). AD-safe (eltype(θ) generic).

source

DRM.simulate_coevolution Function

simulate_coevolution(phy, β, Λ, σ_res; nrep, rng) -> (; Y, X, species)

Simulate nrep observations per tip from the general-q coevolution model on phy: tip random effects vec(B) ~ N(0, Q_cond⁻¹ ⊗ Λ) drawn via the sparse Cholesky, plus Xβ and N(0, diag(σ_res²)) residuals. β is k × q. X = [1 x] with a single standard-normal covariate (k = 2).

source

Bridge and mixed-family payload helpers

DRM._bridge_dpars Function

_bridge_dpars(fit)

Per-observation distributional parameters on the response scale, keyed by drmTMB dpar name ("mu", "sigma", …).

This is what drmTMB's post-fit surface actually consumes. fitted_distribution() — the hub for qq_plot(), worm_plot(), centile_chart() and exceedance() — builds its d/p/q closures from fitted_distribution_params(), which calls predict_parameters(object, dpar = <all dpars>, type = "response") and needs one column per dpar. The R side supplies the density machinery from its own family tables; the only thing it cannot derive is the fitted parameter values.

Covers the in-sample case (R's newdata = NULL). Fresh-data prediction goes through predict_parameters(fit, newdata), which is a separate payload.

A dpar is not fitted(). For a mixture family the two differ, and feeding the wrong one produces a wrong density silently because both are in range. drmTMB's mu dpar for zero_one_beta is the interior beta component mean plogis(eta_mu), which it feeds to drm_beta_shapes(mu, sigma); DRM.jl stores that as beta_mu and puts the unconditional mean (1 - zoi) * mu + zoi * coi — the right answer for fitted() — in means[:mu]. The override below repairs that one family.

Checked against drmTMB's full dpar table (R/family-dpq.R): every other family DRM.jl implements already agrees, including truncated NB2, whose means[:mu] is the untruncated mean and so is already the correct dpar. DRM.jl has no zi/hurdle families, the other place this trap lives.

source

DRM._bridge_dpars_newdata Function

_bridge_dpars_newdata(fit, newdata)

Distributional parameters on the response scale for fresh rows.

R's predict_parameters(object, newdata = ..., type = "response") is what fitted_distribution(object, newdata = ...) calls, and the julia-engine vignette records the gap this closes: fresh-data Julia prediction is currently limited to location parameters.

Unlike the in-sample block this reads the FORMULA rather than the stored means/scales, so each parameter comes back as its own linear predictor pushed through its link — which is already the dpar drmTMB wants (for zero_one_beta, mu here is plogis(eta_mu), the interior beta mean, with no override needed).

source

DRM._bridge_meta_parts Function

_bridge_meta_parts(fit)

For a meta-analysis fit (meta_V(v) on mu), the between-study heterogeneity tau and the known sampling variances V_known; nothing otherwise.

Why this is a recovery rather than a stored field. gaussian_meta.jl stores scales[:sigma] = sqrt(V + tau^2) — the TOTAL per-study SD, which is what simulate() needs. But drmTMB's meta sigma dpar is the heterogeneity ALONE, with V_known supplied separately and the density forming sqrt(V_known + sigma^2) itself. Emitting the total as sigma and a V_known would double-count the sampling variance.

Both are recoverable exactly from what the fit already holds — tau from the sigma coefficient (already the right dpar) and V = total^2 - tau^2, verified to 1.1e-16 — so no field, no struct change, and no change to sigma()'s public contract is required.

Returns nothing when the sigma block carries predictors, because per-row tau then needs the design matrix and DrmFit does not retain the data. That is a declared boundary, not a silent approximation.

source

DRM._bridge_trials Function

_bridge_trials(fit)

Per-row binomial denominator, or nothing when the family has none.

fitted_distribution_params() attaches params$trials for the binomial and beta_binomial model types; without it the R side cannot evaluate those densities on a Julia fit.

source

DRM._bridge_mixture_family_hint Function

Actionable refusal text for a zero-inflated / hurdle family SPELLING.

Zero-inflation (zi) and hurdle (hu) are FORMULA parts here, not families: the count family stays poisson / nbinom2 and a keyed zi ~ … or hu ~ … entry adds the mixture (see _bridge_formula above, test_zi.jl, test_hurdle.jl). drmTMB spells them the same way – its drm_family_type() returns "poisson" for a zero-inflated Poisson, and "zi_poisson" is only its post-fit model_type – so the R bridge never sends these tags. A caller writing drm_bridge by hand reasonably might, and a bare "unsupported family" does not tell them the spelling that works.

Deliberately NOT an alias. Mapping "zi_poisson" to Poisson() would fit a PLAIN Poisson, with no error, whenever the caller omitted the zi ~ part – a silent wrong answer in place of a loud refusal.

Returns nothing (so the caller falls through to the plain message) unless the tag is a recognised mixture prefix on a count family this bridge accepts.

source

DRM.drm_bridge_q2_known_precision Function

drm_bridge_q2_known_precision(; Y, X, group, Q,
                              structured_type = "relmat",
                              precision_source = nothing,
                              options = Dict())

Private diagnostic boundary for a restricted q2 known-precision provider payload. This consumes Q as a precision matrix directly through make_coevo_problem_from_precision; it does not invert or relabel Q as a covariance matrix, and it does not imply formula, slope, REML, or interval support. Provider identity is deliberately narrow: structured_type = "relmat" records precision_source = "Q", while structured_type = "animal" records precision_source = "Ainv".

source

DRM.drm_bridge_q2_phylo Function

drm_bridge_q2_phylo(; Y, X, species, tree, options = Dict())

Private diagnostic boundary for the restricted q2 phylogenetic coevolution point export. This intentionally bypasses the public formula bridge because the general-q coevolution model is diagonal-residual Gaussian evidence, not the full bivariate rho12 q2 route. The return payload is the same primitive dictionary shape consumed by the row-contract tests.

source

DRM.fit_mixed_family Function

fit_mixed_family(; y1, X1, fam1, y2, X2, fam2,
                   trials1=ones(n), trials2=ones(n),
                   Xsigma1=ones(n,1), Xsigma2=ones(n,1), K=32, g_tol=1e-6)

Fit the cross-family bivariate model (shared per-observation latent) and return a NamedTuple with fixed effects β1/β2, loadings λ1/λ2, per-axis scale σ1/σ2 (NaN for dispersionless axes), dispersion sub-model coefficients βσ1/βσ2 (log σ scale, empty for dispersionless axes), link-scale variances v1/v2, the latent-scale correlation rho_latent, loglik, converged, iterations, and the two family instances fam1/fam2 (carried so the post-fit accessors in mixed_family_postfit.jl can recover each axis's inverse link).

fam1/fam2 are DRM family instances. Supported: Gaussian, Poisson, Binomial, NegBinomial2, Beta, Gamma. Dispersion-carrying families (Gaussian/Beta/Gamma/NB2) all carry the scale σ in a per-observation log σ SUB-MODEL log σ_{k,i} = (Xσ_k · β_σk)_i with size(Xσ_k, 2) coefficients; Poisson/Binomial are dispersionless. Each family maps σ to its own dispersion (Beta φ = 1/σ², Gamma α = 1/σ², NB2 size θ = 1/σ²), matching the univariate fitters and drmTMB (#315/#316). Xsigma1/Xsigma2 default to a single intercept column (ones(n,1)), reproducing the intercept-only scalar dispersion. For dispersion-carrying families σ1/σ2 returns a representative scale exp(mean_i log σ_{k,i}) (recover the NB2 size as θ = 1/σ1²); the full sub-model is in βσ1/βσ2. trials* are Binomial denominators (ignored otherwise).

source

DRM._mf_nparams Function

_mf_nparams(fit) -> Int

Number of free parameters in a fit_mixed_family fit: the two fixed-effect blocks β1/β2, the two latent loadings λ1/λ2, and the two dispersion sub-model blocks βσ1/βσ2 (empty for dispersionless axes). This equals the length of the internal θ vector the optimiser minimised.

source

Bridge label metadata internals (no stability guarantee)

DRM._BridgeFormulaLabels Type

Internal, typed provenance for bridge-only public coefficient labels.

atoms maps collision-safe materialised Julia symbols to their exact R term spellings. function_labels retains source spelling for every admitted scalar function term, including redundant parentheses. Neither changes a formula, matrix, or fit.

source

Location-scale and structured-fit kernels

DRM._fit_bivariate_q4_structured Function

_fit_bivariate_q4_structured(...)

q=4 PLSM front end for level-indexed structured providers (relmat / animal / spatial) — issue #189. Reuses fit_q4_sparse_tmb via make_problem_from_Q; does not rewrite the verified Laplace engine.

Spatial uses a fixed range (spatial_range, default = mean pairwise distance). Non-tree bootstrap_sigma_a is deliberately unsupported.

source

DRM._fit_cumulative_phylo_laplace Function

_fit_cumulative_phylo_laplace(fam, y, Xμ, K, labels, tree, nmμ, grp, g_tol; se)

Sparse-Laplace CumulativeLogit() fit with a phylogenetic random intercept phylo(1 | grp) on the mean linear predictor η (#563 S8 follow-on). Cutpoints stay ordinary fixed effects (shared across the tree, exactly as in _fit_cumulative_ranef's GHQ route); only the mu intercept varies by tip. See the file-level comment above for the exact-gradient design.

source

DRM._fit_corr_locscale Function

_fit_corr_locscale(fam, kind, rk, y, Xμ, Xψ, xs, gidx, G, nmμ, nmσ, grp;
                   link, se, g_tol, respobs, trials, Q)

Fit a correlated (rk == :corr) or independent (rk == :slope) random intercept/slope on the mean of a non-Gaussian family through the unified q2 location–scale Laplace core, and wrap it as a DrmFit with the GHQ :recov block. Xψ is the scale fixed-effect design (zeros(n,0) for the mean-only Poisson). link is the mean inverse link (:log/:logit/:identity). Q is the G×G group-level precision; default sparse(I,G,G) (i.i.d. groups, cluster 1 byte-identical behaviour). A structured Q (from _general_cov_setup or _locscale_phylo_setup) routes the kron(Q,Λ⁻¹) prior through the same spine (cluster 3: structured non-Gaussian random slopes).

source

DRM._fit_fixed_gaussian_reml Function

_fit_fixed_gaussian_reml(fam, y, Xμ, Xσ, nmμ, nmσ, g_tol) -> DrmFit

REML fit of the fixed-effect Gaussian location–scale model. Profiles out β_μ by weighted least squares and optimises the restricted log-likelihood over β_σ; β̂_μ is recovered as the final WLS estimate. The returned DrmFit carries estim_method = :REML,reml_loglik, andml_loglik(the plain ML log-lik at the REML parameters, for reference). Internal — reached viadrm(...; method = :REML).

source

DRM._fit_gaussian_locscale_phylo Function

_fit_gaussian_locscale_phylo(fam, y, Xμ, Xψ, gidx, G, Q, nmμ, nmσ, grp;
                              coupled, asymmetric, se, g_tol) -> DrmFit

Fit a UNIVARIATE Gaussian location-scale model with a phylogenetic random effect on BOTH axes (by default SEPARATE / uncorrelated) or on the σ axis only (asymmetric).

Modes controlled by kwargs:

• coupled = false (default): SEPARATE block Λ = diag(L11², L22²), L21 ≡ 0. Returns DrmFit with :mu, :sigma, :resd_mu (logL11), :resd_sigma (logL22).

• coupled = true: FREE L21, full 2×2 Λ with mean↔σ correlation. Returns DrmFit with :mu, :sigma, :recov (logL11, logL22, L21).

• asymmetric = true: σ-phylo only (mean fixed effects, no mean phylo RE). Returns DrmFit with :mu, :sigma, :resd_sigma (logL22).

The kernel is Val(:gaussian_mean): η = mean, ψ = log σ, integrating the Gaussian location-scale likelihood through the q=2 augmented-state Laplace spine.

source

DRM._fit_locscale Function

_fit_locscale(kind, y, Xμ, Xψ, gidx, G, Q; ...)

Fit the q=2 non-Gaussian location–scale model. Returns a named tuple with the packed estimate θ, the fixed effects beta_mu / beta_psi, the 2×2 group-level covariance Lambda, the marginal nll, and a converged flag.

source

DRM._fit_ranef_gaussian Function

_fit_ranef_gaussian(..., reml=false) -> DrmFit

Gaussian location–scale with one mean random intercept (1 | g) on the Woodbury spine. reml=false (default) is ML, byte-for-byte with the historical nll. reml=true (#439) adds ½ logdet(Xμ′ V⁻¹ Xμ) − ½ pμ log(2π) to that nll (Patterson–Thompson). Reached via drm(...; method = :REML) for a single intercept only — σ-RE, slopes, and multi-ranef stay rejected. Worked example: test/test_reml_ordinary_ranef.jl (not in the default suite yet).

source

DRM._fit_sigma_axis_re Function

_fit_sigma_axis_re(fam, kind, y, Xμ, Xψ, gidx, G, nmμ, nmσ, grp; link, se, g_tol,
                   obs_prop, trials)

Fit a standalone σ-axis random intercept sigma ~ 1 + (1|g) via the q=2 location–scale Laplace core. Mean is fixed-effects only (no mean RE). The latent is placed on the scale axis: Zη = 0, Zψ[:,1] = 1. Only logL11 (= log τ) is optimized; the unused latent axis-2 is pinned to a fixed tiny variance ε = 1e-6. Returns a DrmFit with a :resd block (log τ), so re_sd(fit)[:g_logsigma] = τ (the _logsigma suffix marks this SD as living on the log-σ scale axis).

source

DRM._general_cov_setup Function

_general_cov_setup(C, labels) -> (Q, leaf_node)

Resolve a user-supplied per-group covariance/relatedness matrix C (G×G, SPD) to the (precision Q, leaf_node) pair the sparse-Laplace spine consumes. C is first rescaled to a unit-diagonal correlation R (so the recovered :resd block is the random-effect SD σ_b, with prior b ~ N(0, σ_b² R), matching the phylo convention), then Q = R⁻¹. leaf_node maps each observation to its group level via _group_index — the relatedness/animal-model/spatial analogue of the tree's leaf-to-node map. Mathematically identical to the phylo path with the tree precision swapped for an arbitrary PD precision.

source

DRM._ls_components Function

_ls_components(Λ) -> NamedTuple

Named group-level summaries of the 2×2 covariance Λ: the mean-axis SD, the scale-axis SD, and the mean↔scale correlation ρ_a.

source

DRM._ls_marginal_grad Function

_ls_marginal_grad(kind, y, Xμ, Xψ, gidx, G, Q, θ) -> Vector

Exact gradient of _ls_fit_nll at the packed θ = [βμ; βψ; λ(3)]. Returns an all-NaN vector if the inner Laplace mode fails to converge (rare; this is the infeasibility signal that pairs with the _ls_fit_nll value sentinel — a gradient-based optimiser rejects a NaN step rather than mistaking a zero gradient for stationarity, see #314). O(p) in the number of groups.

source

DRM._ls_marginal_nll Function

_ls_marginal_nll(kind, y, η0, ψ0, gidx, G, P; a0=nothing)

Laplace-approximate marginal negative log-likelihood for the q=2 location–scale model. Returns (nll, â, ok): the marginal NLL, the inner mode, and a success flag (the inner Newton solve can fail at extreme parameters).

source

DRM._ls_obs_information Function

_ls_obs_information(kind, y, Xμ, Xψ, gidx, G, Q, θ; h=1e-5) -> Symmetric

Observed information ∂²M/∂θ² of the Laplace marginal at θ, as the symmetrised central finite-difference Jacobian of the exact gradient _ls_marginal_grad.

source

DRM._ls_profile_ci Function

_ls_profile_ci(kind, y, Xμ, Xψ, gidx, G, Q, θ̂; kwargs...)

Compatibility wrapper returning only (lower, upper) from _ls_profile_ci_result. Use the structured result to inspect endpoint failure, no-crossing, and nuisance diagnostics.

source

DRM._ls_profile_ci_result Function

_ls_profile_ci_result(kind, y, Xμ, Xψ, gidx, G, Q, θ̂; idx, level=0.95,
                      nll_min=nothing, se=nothing)

Profile-likelihood CI for packed parameter idx, inverting 2(ℓ̂ − ℓ_profile) = χ²₁(level) by a Venzon–Moolgavkar-style guarded-Newton root-find (the profile deviance and its envelope-theorem slope, bracket-safeguarded). A directional ±Inf may mean either that no finite crossing was found in the bounded search range or that an endpoint failed; _ls_profile_ci_result distinguishes them. se (a Wald SE for idx) seeds the bracket width; if omitted it is derived from the observed information.

Zη/Zψ are the mean/scale latent loadings and MUST match the model that was fit: the profiler reconstructs the marginal from (kind, y, Xμ, Xψ, gidx, G, Q) plus these loadings, so passing the wrong loadings profiles a different model (#325.4). They default to the canonical loadings (Zη=[1 0], Zψ=[0 1]) — the only case wired through LocScaleObjective today; a non-canonical (slope-axis) fit must pass its own Zη/Zψ here for the CI to be correct.

source

DRM._phylo_mean_laplace_hetero_fg Function

_phylo_mean_laplace_hetero_fg(kind, aux_from, n, Xμ, Xσ, leaf_node, Q, logdetQ, θ; …)

Heteroscedastic generalisation of _phylo_mean_laplace_nuisance_fg (#164): the single scalar dispersion nuisance θσ is replaced by a per-observation log-dispersion linear predictor ησ = Xσ·βσ, so θ = [βμ(pμ); βσ(pσ); logσ]. aux_from(ησ::Vector) returns a per-observation aux whose dispersion is exp.(ησ) (a vector kind, e.g. Val(:nb2_hetero)).

The σ-axis kernel derivatives (nval, nr, nw) are taken w.r.t. each observation's log r_i, so the βσ gradient is the Xσ-chained version of the scalar nuisance gradient: where the scalar code accumulates one gnuis, this accumulates a pσ-vector gν with each per-observation contribution weighted by Xσ[i, k], and the implicit cross-term crossν becomes a q × pσ matrix — structurally identical to how the mean axis already chains r/w through Xμ. A one-column constant Xσ reproduces _phylo_mean_laplace_nuisance_fg exactly.

source

DRM._phylo_mean_leaf_index Function

_phylo_mean_leaf_index(phy::AugmentedPhy, labels) -> Vector{Int}

Map per-row group labels to the tree's raw leaf index (1:phy.n_leaves) — the convention make_loc_problem's species= keyword expects. Mirrors the three-tier match _poisson_phylo_setup already uses for the Laplace phylo routes, most to least specific: 2. leaf NAME (phy.leaf_names) — subset-tolerant: labels may name only SOME of the tree's leaves (e.g. after a caller drops missing-response rows upstream — #482). A leaf that never appears in labels simply gets no observation and stays in the phylo PRIOR only, exactly like the σ-phylo route's documented "species whose every row is missing stays in the prior with no likelihood term" (gaussian_core.jl, the σ-phylo missing-response block).

3. integer tip index 1:phy.n_leaves.

4. positional (_group_index, first-seen order) — the only tier with no name information to anchor to, so it requires labels to visit every tip exactly once (G == phy.n_leaves); this is the pre-#482 behaviour, kept as a last-resort fallback.

Throws an ArgumentError naming the mismatch if none of the three tiers apply (e.g. group labels that are neither tree tip names nor valid integer indices, and do not even form a complete one-level-per-tip bijection).

source

DRM._reml_normalise Function

_reml_normalise(reml_ll, n_beta)

Add the (n_beta/2)·log(2π) normalising constant to an unnormalised Patterson–Thompson restricted log-likelihood, so DRM.jl's bivariate REML routes report on the same scale as lme4, glmmTMB, TMB — and as DRM.jl's own univariate REML routes, which have always added it (#477).

n_beta counts only the marginalised fixed effects. Non-finite input passes through unchanged so -Inf barriers and NaN sentinels keep their meaning.

source

DRM._vcov_from_hessian Function

_vcov_from_hessian(H; context = "")

Covariance matrix from an observed-information (Hessian) matrix, guarded against boundary degeneracy.

Symmetrises H, then inverts it — unless it is numerically singular, in which case it falls back to the Moore–Penrose pseudo-inverse and warns, naming the parameter coordinates that are flat. The decision is made from the eigenvalues rather than from whether inv happens to throw, so the result does not depend on the LAPACK build or CPU.

The pseudo-inverse keeps a fit usable, but standard errors for the flagged coordinates are not trustworthy: at a variance boundary the sampling distribution is not the usual normal approximation. Use method = :profile or the bootstrap entry points for those targets (see src/inference.jl), or the χ̄² boundary machinery in src/chibar.jl.

source

DRM.q2_reml_phi_len Function

Length of the q=2 REML outer parameter vector phi = (lc(Λ), logσ1, logσ2, ηρ).

source

Paired whitening internals (no stability guarantee)

These private helpers implement the transformed latent-state route used by coupled non-Gaussian location-scale fitting, observed information, and marginal bootstrap simulation. They are documented here so their contracts remain visible to Documenter's completeness check; they are not public API.

DRM._LSWhitenedSeed Type

Owned transformed warm start; acceptance is intentionally not stored.

source

DRM._ls_whitened_eval Function

_ls_whitened_eval(kind, y, Xμ, Xψ, gidx, G, Q, theta, Zη, Zψ;
                  seed=nothing, gradient=true)

Private paired value/gradient route in whitened latent coordinates. A returned seed is an owned initial guess only: every invocation re-solves and certifies at its current theta, Q, data, and loadings.

source

DRM._ls_whitened_information Function

_ls_whitened_information(kind, y, Xμ, Xψ, gidx, G, Q, theta,
                          Zη=canonical, Zψ=canonical; h=1e-5, seed=nothing)

Observed information for the private paired whitening route. This has the same central-difference and symmetrisation contract as _ls_obs_information, but every theta ± h e_k call receives a separately owned initial-guess seed. Failure of either certified paired evaluation makes the whole information matrix unusable (all NaN), rather than mixing paired and legacy gradients.

source

DRM._ls_whitened_vcov Function

_ls_whitened_vcov(args...; kwargs...)

Wald covariance from private paired observed information. Singular information returns nothing, as in _ls_vcov; failed/non-finite paired information is also refused explicitly instead of passing it to matrix inversion.

source

DRM._ls_bootstrap_effect Function

Draw group-by-axis effects using a frozen sparse triangular precision factor.

source

Location-scale inner-mode acceptance

For the non-Gaussian location-scale Laplace engine, a positive-definite latent Hessian alone does not certify a mode. The inner solver must also return finite coordinates and a fresh finite gradient satisfying norm(gradient) <= tol * (1 + norm(a)). The default remains tol = 1e-9 with at most 200 iterations. Failure to meet this condition is a failed inner solve, not an accepted likelihood evaluation.

Near a mode, objective rounding can hide a useful Newton step. A separate local polishing check permits a full, undamped step with an objective increase of at most four units in the last place (ULPs), provided the actual displacement is small, the gradient strictly decreases and meets the same stationarity criterion, and the trial Hessian is finite and positive definite.

For NB2 and Gamma, a larger rounding discrepancy can trigger a second check under the same safeguards. It estimates the objective change from directional derivatives and a quadrature integral, avoiding subtraction of two nearly equal objective values. The prior contribution retains multiplication and summation residuals; its discarded rounding terms are tracked alongside the existing data and quadrature error estimates. The estimated change plus its numerical error margin must be negative. This fallback is unavailable at or across the kernels' predictor-clamp boundaries, or when the tracked arithmetic risks overflow or underflow. An unavailable or inconclusive estimate leaves ordinary backtracking in place. The margin is an engineering estimate, not a proven error bound.

Neither polishing check guarantees exact objective descent or a global optimum. Other steps retain the ordinary descent check; coordinate-identical trials do not count as progress.

This helper is internal and has no stability guarantee.

DRM._ls_inner_estimated_change Function

_ls_inner_estimated_change(kind, y, η0, ψ0, gidx, G, P, Zη, Zψ, a, trial)

Estimate the local joint-objective change for an otherwise certified Newton trial when direct Float64 subtraction is unreliable. This is a guarded numerical estimate, not a proof of descent: unsupported families, clamp-boundary segments, and non-finite arithmetic return nothing. margin = Q8 + E must be negative before the caller may use the estimate.

source

Phylogenetic group-index internals (no stability guarantee)

DRM._lss_phylo_group_index Function

_lss_phylo_group_index(tree, labels, grp) -> phy, gidx, G

Resolve a Gaussian LSS phylogenetic grouping column against its tree. Integer labels are positional leaf indices 1:p; text labels match phy.leaf_names exactly. The full input must contain every leaf before any missing response is removed, so an all-missing response tip remains a prior-only tree state instead of changing the covariance dimension. This helper deliberately accepts only the shipped AugmentedPhy and Newick-string tree contracts; arbitrary sigma_phy_dense providers are not qualified for labelled LSS identity.

source

Non-Gaussian sparse-Laplace kernels

DRM._fit_beta_relmat_laplace Function

_fit_beta_relmat_laplace(fam, y, Xμ, Xσ, C, labels, nmμ, nmσ, grp, g_tol; se)

Beta sparse-Laplace fit with a general user-supplied PD covariance C on the mean (logit) random intercept (relmat/animal/spatial(1 | grp)), with the precision φ (via the sigma slot, σ = 1/√φ) a fixed nuisance parameter. Reuses the verified phylo nuisance spine via _general_cov_setup; only the prior precision differs (C⁻¹ vs the tree topology). Exact O(p) gradient carries over (#167).

source

DRM._fit_betabinomial_crossed_laplace Function

_fit_betabinomial_crossed_laplace(fam, s, ntr, Xμ, comps, nmμ, nmσ, g_tol; se)

Beta-binomial sparse-Laplace fit with two crossed random intercepts on the logit mean, e.g. (1 | g) + (1 | h), constant-σ (overdispersion) only. Reuses the verified Beta-family crossed nuisance Laplace spine (_fit_crossed_mean_laplace_nuisance) with the :betabinomial_fixed kernel (#166; see docs/dev-log/plans/2026-08-02-166-betabinomial-kernel-design.md).

source

DRM._fit_betabinomial_phylo_laplace Function

_fit_betabinomial_phylo_laplace(fam, s, ntr, Xμ, labels, tree, nmμ, nmσ, grp, g_tol; se)

Beta-binomial sparse-Laplace fit with a phylogenetic random intercept phylo(1 | grp) on the logit mean, constant-σ (overdispersion) only. Reuses the verified Beta-family nuisance Laplace spine (_fit_phylo_mean_laplace_nuisance) with the :betabinomial_fixed kernel — shifted digamma/trigamma/polygamma arguments (s+a, n-s+b, n+a+b) replace Beta's (a, b) (#166; see docs/dev-log/plans/2026-08-02-166-betabinomial-kernel-design.md).

source

DRM._fit_gamma_relmat_laplace Function

_fit_gamma_relmat_laplace(fam, y, Xμ, Xσ, C, labels, nmμ, nmσ, grp, g_tol; se)

Gamma sparse-Laplace fit with a general user-supplied PD covariance C on the mean random intercept (relmat/animal/spatial(1 | grp)), with the shape α (via the sigma slot, σ = 1/√α) a fixed nuisance parameter. Reuses the verified phylo nuisance spine via _general_cov_setup; only the prior precision differs (C⁻¹ vs the tree topology). Exact O(p) gradient carries over (#167).

source

DRM._fit_nb2_relmat_laplace Function

_fit_nb2_relmat_laplace(fam, y, Xμ, Xσ, C, labels, nmμ, nmσ, grp, g_tol; se)

NB2 sparse-Laplace fit with a general user-supplied PD covariance C on the mean random intercept (relmat/animal/spatial(1 | grp)), with the dispersion θ (the sigma slot) a fixed nuisance parameter. Reuses the verified phylo nuisance spine via _general_cov_setup: the only difference from the phylo route is that the prior precision comes from C⁻¹ instead of the tree topology. Exact O(p) gradient and Takahashi log-det derivatives carry over unchanged (#167).

source

DRM._fit_poisson_crossed_laplace Function

_fit_poisson_crossed_laplace(fam, y, Xμ, comps, nmμ, g_tol)

Internal engine-lane fitter for Poisson random-intercept GLMMs with one or more independent scalar random-effect components. comps is a vector of (w, gidx, G, label) tuples, matching the Gaussian multi-RE fitter.

source

DRM._fit_poisson_phylo_laplace Function

_fit_poisson_phylo_laplace(fam, y, Xμ, labels, tree, nmμ, grp, g_tol)

Internal sparse-Laplace fitter for Poisson() with a phylogenetic random intercept phylo(1 | grp) on the mean. The tree is represented by the root-conditioned augmented precision, so the latent state has one effect per non-root tree node and the mode/logdet computations stay sparse.

source

DRM._fit_poisson_relmat_laplace Function

_fit_poisson_relmat_laplace(fam, y, Xμ, C, labels, nmμ, grp, g_tol; se)

Poisson sparse-Laplace fit with a general user-supplied PD covariance C on the mean random intercept (relmat/animal/spatial(1 | grp)). Reuses the verified phylo Laplace spine via _general_cov_setup: the only difference from the phylo route is that the prior precision comes from C⁻¹ instead of the tree topology. Exact O(p) gradient and Takahashi log-det derivatives carry over unchanged.

source

DRM._fit_poisson_spatial_coord Function

_fit_poisson_spatial_coord(fam, y, Xμ, labels, coords, nmμ, grp, g_tol; se)

Poisson spatial(1 | grp) fit with a coordinate-based exponential-kernel spatial covariance C(ρ) = exp(-d/ρ) and the range ρ estimated jointly (#270). The site coordinates (coords, a G×2 matrix, one row per group level in first-seen order) give the pairwise distances d; the random intercept is b ~ N(0, σ_b² C(ρ)) on log λ. The outer parameter vector is θ = [β_μ; log σ_b; log ρ], so ρ is a genuine hyperparameter whose gradient flows through C(ρ) into the Laplace marginal — not a fixed-range or profile-over-grid approximation.

The marginal is the Poisson Laplace approximation (same convention as the verified sparse general-covariance path), written in AD-traceable operations so ForwardDiff supplies the exact outer gradient and the covariance Hessian. Recovered blocks: :mu (fixed effects), :resd (log σ_b), :range (log ρ).

Note

The range ρ is only weakly identified from a single spatial realization; :resd (the spatial SD) and the fixed effects are the robustly recoverable pieces. For a precomputed spatial covariance, use spatial(1 | grp) with K = … instead (the verified O(p) sparse path).

source

DRM._poisson_crossed_laplace_fg Function

_poisson_crossed_laplace_fg(y, Xμ, gidx, G, hidx, Hh, lf, θ; grad, b0, newton_tol, newton_maxiter)

Module-level f/g evaluation of the Poisson crossed-random-intercepts Laplace marginal NLL at a single θ, used both by _fit_poisson_crossed_intercepts_laplace and by the standing FD-vs-analytic gradient gate (#165). Hoisting it out of the fit closure lets the gate drive a controlled, tightly-converged, warm-started inner mode (the same recipe the q4 Q-gate / the Poisson-phylo gate use to reach ≤ 1e-6). Returns (val[, grad], b, ok); on a non-PD inner solve ok = false.

source

DRM._poisson_phylo_laplace_fg Function

_poisson_phylo_laplace_fg(y, Xμ, leaf_node, Q, logdetQ, lf, θ; grad, b0, newton_tol, newton_maxiter)

Module-level f/g evaluation of the Poisson phylogenetic-Laplace marginal NLL at a single θ, used both by _fit_poisson_phylo_laplace and by the standing FD-vs-analytic gradient gate (#165). Hoisting it out of the fit closure lets the gate drive a controlled, tightly-converged, warm-started inner mode (the same recipe marginal_and_exact_grad / the q4 Q-gate use to reach ≤ 1e-6).

The marginal is

L(θ) = data(b̂) + ½ σ⁻² b̂'Q b̂ + q logσ − ½ logdet Q + ½ logdet H,
H = σ⁻² Q + diag(Σ_{i∈leaf} μ_i).

Because b̂ solves ∂(data+prior)/∂b = 0, the total θ-gradient is the explicit part (b̂ frozen) plus the implicit logdet correction through db̂/dθ = −H⁻¹ ∂²(data+prior)/∂b∂θ. For Poisson the data Hessian weight and its η-derivative coincide (d²ℓ = d³ℓ = μ), so a single μ drives both the logdet trace and the cross term — this is the family-specific third-derivative term the IFT introduces. Returns (val[, grad], b, ok); on a non-PD inner solve ok = false.

source

Helpers referenced by source docstrings

These nonpublic helpers have no source docstrings; the descriptions below keep references from the documented kernels navigable without creating new APIs.

_group_index

Maps observation labels to integer group indices in first-seen order and returns the indices together with the number of distinct groups.

_fit_phylo_mean_laplace_nuisance

Prepares a tree's root-conditioned precision and observation-to-leaf map, then delegates to the general-precision mean-effect Laplace fitter.

_fit_crossed_mean_laplace_nuisance

Fits the crossed mean-effect Laplace objective for a supplied family and nuisance-parameter specification, using the prepared group indices.

Prepared missing-predictor development route

Experimental

Exported for evaluation; fenced for v1.0 (D-181). API and numerics may change; not covered by the R-parity scoreboard.

Limited developer interface

This prepared-array interface covers a Gaussian response with one Gaussian or Bernoulli predictor, or two independent Gaussian predictors. The likelihood has two development frontends: the joint formula frontend and drmTMB(..., engine = "julia") through drm_bridge_joint. They admit a Gaussian identity-link response, one or two bare additive mi() terms, and complete fixed-effect exogenous designs. Grouped predictors, further predictor families, random or structured effects, REML, and all other missing-predictor models remain outside this route.

The R preparation route supports response = "drop" or "include" with predictor = "model"; its response-drop preprocessing is deliberately not presented as native-TMB response-policy parity. Profile and bootstrap intervals are unavailable. Gaussian predictor-SD Wald intervals use a delta transformation on the natural SD scale and may cross zero; they are not a claim of native interval parity or interval coverage. Conditional variance at fixed parameters is not the native R imputed() standard error.

Design matrices must be complete. Only the modelled predictors and y may contain missing. For the one-predictor interface, the parameter order is mean coefficients, the coefficient of x, residual log-SD coefficients, predictor coefficients, and (Gaussian only) predictor log-SD. This last coordinate is a log-SD, not R's natural-scale sigma_mi_x.

using DRM
x = Union{Missing,Float64}[0.8, missing, 1.0, missing]
y = Union{Missing,Float64}[1.7, -0.2, missing, missing]
z = [-0.6, 0.3, 0.8, -0.1]
X = hcat(ones(4), z)
model = prepared_joint_model(y, x, X, ones(4, 1), X;
    predictor = :gaussian, mu_names = ["(Intercept)", "z"],
    predictor_names = ["(Intercept)", "z"], original_row = [17, 4, 81, 29])
theta = [0.2, -0.35, 0.7, 0.1, 0.0, 0.25, log(0.9)]
row_loglik = prepared_joint_rowloglik(model, theta)
moments = prepared_joint_conditional_moments(model, theta)
@assert row_loglik[4] == 0.0  # Both variables missing: integral equals one.
@assert moments.mean[1] == x[1]
(model.original_row, row_loglik, moments)

([17, 4, 81, 29], [-2.6077661253412465, -1.1969370291436168, -1.208639745941908, 0.0], (mean = [0.8, -0.046752557737715306, 1.0, -0.025], variance = [0.0, 0.6113418698207638, 0.0, 0.81], status = [:observed, :gaussian_posterior, :observed, :predictor_only]))

This example evaluates a supplied parameter vector; it does not fit its four rows. On an identifiable dataset, fit_prepared_joint(model) estimates the parameters. Inspect joint_missing_summary(result).optimizer_status and .covariance_status separately. Its .uncertainty_status is :not_computed: this low-level summary does not calculate standard errors. Use imputed for native-shaped imputation summaries and their own uncertainty status.

Two independent Gaussian predictor models

The two-predictor array interface uses an n × 2 predictor matrix and a tuple of two complete predictor designs. It shares the prepared likelihood operations above. Direct formula and R bridge admission use this same kernel; the bounded route does not establish full native-R parity.

Independent predictor models can produce correlated conditional imputations: when both predictors are missing, the observed response informs their joint values. The returned covariance retains this dependence.

using DRM
x = Union{Missing,Float64}[0.8 missing; missing missing; 1.0 0.3; missing missing]
y = Union{Missing,Float64}[1.7, -0.2, missing, missing]
z = [-0.6, 0.3, 0.8, -0.1]
X = hcat(ones(4), z)
model = prepared_joint_model(y, x, X, ones(4, 1), (X, X);
    predictor_variables = (:x1, :x2),
    mu_names = ["(Intercept)", "z"],
    predictor_names = (["(Intercept)", "z"], ["(Intercept)", "z"]))
# beta, b1, b2, residual log-SD, alpha1, logtau1, alpha2, logtau2
theta = [0.2, -0.35, 0.7, -0.4, 0.1, 0.0, 0.25, log(0.9), -0.2, 0.15, log(0.8)]
moments = prepared_joint_conditional_moments(model, theta)
@assert prepared_joint_rowloglik(model, theta)[4] == 0.0
@assert moments.covariance[2, 1, 2] > 0  # Opposite response slopes.
(moments.mean, moments.covariance[2, :, :])

([0.8 -0.4087367219407032; -0.0599370185518045 -0.09407605511594014; 1.0 0.3; -0.025 -0.21500000000000002], [0.62316412815904 0.08435623137792908; 0.08435623137792908 0.6019132359210763])

This is a parameter-point calculation, not a fit of four observations. On an identifiable dataset, fit_prepared_joint(model) fits the same prepared model. Select a predictor explicitly with imputed(result; variable = :x1) or :x2. Its Gaussian imputation SE combines conditional variance with first-order parameter uncertainty using the full fitted covariance; it does not establish interval coverage or provide multiple-imputation draws.

DRM.PreparedTwoJointGaussianModel Type

Validated, formula-free Gaussian-response model with two independent Gaussian predictors.

source

DRM.PreparedTwoJointGaussianFit Type

Fitted exact prepared model for two independent Gaussian missing predictors.

source

DRM.JointTwoMissingMetadata Type

Row-level conditional moments for two Gaussian missing predictors.

source

Ordinal and categorical predictor models

The finite-state prepared interface accepts one ordinal or categorical predictor and a Gaussian response. It integrates over every possible state when the predictor is missing. Supply an n × K × p array containing the complete mean design for each row and state; ordinal contrasts and categorical dummy variables must already be encoded in this array. The direct mi() formula frontend and R bridge construct this state design for their bounded finite-state routes; they do not establish full native-R prediction or accessor parity.

using DRM
levels = ["low", "middle", "high"]
x = Union{Missing,String}["middle", missing, "high", missing]
y = Union{Missing,Float64}[0.4, -0.3, missing, missing]
z = [-0.6, 0.3, 0.8, -0.1]
Xstate = zeros(4, 3, 2)
for i in 1:4, k in 1:3
    Xstate[i, k, :] = [1.0, k - 2.0]
end
model = prepared_joint_model(y, x, Xstate, ones(4, 1), reshape(z, 4, 1);
    predictor = :ordinal, levels = levels, variable = :severity)
# Mean coefficients, log-SD, predictor slope, first cut, log cut spacing.
theta = [0.1, 0.5, log(0.7), 0.3, -0.6, log(1.2)]
moments = prepared_joint_conditional_moments(model, theta)
@assert prepared_joint_rowloglik(model, theta)[4] == 0.0
@assert moments.mean[1] == 2.0
moments.probabilities

4×3 Matrix{Float64}:
 0.0       1.0       0.0
 0.445737  0.332944  0.221319
 0.0       0.0       1.0
 0.361237  0.291253  0.347511

This example evaluates parameters, without fitting four rows. For ordinal predictors, cumulative probabilities use logistic(cutpoint - linear_predictor); the predictor design must not span an intercept. Raw parameters are mean coefficients, residual log-SD coefficients, predictor coefficients, then the first cutpoint and log positive spacings. The raw covariance uses these same coordinates.

For predictor = :categorical, the first declared level is the baseline. Predictor coefficients are ordered by nonbaseline level, then design term; there is no cutpoint block. The prepared constructor accepts two or more states; that low-level admission does not replace native frontend restrictions.

On identifiable datasets, fit_prepared_joint(model) returns state-weighted fitted means and retains posterior probabilities. Ordinal imputed() reports expected scores and conditional score SDs, without adding parameter uncertainty. Categorical imputed() reports the first modal category code; a metric SE is unavailable, with an explicit status. Fit covariance failures take precedence. Neither summary is a multiple-imputation draw or an interval-coverage claim.

The direct formula frontend uses this same kernel. Supply declared levels for textual ordered data; declaring nominal levels fixes the baseline explicitly:

using DRM, Random
rng = MersenneTwister(563)
labels = ["low", "medium", "high"]
codes = repeat(1:3, 30)
z = randn(rng, 90)
x = Union{Missing,String}[labels[k] for k in codes]
y = Union{Missing,Float64}[0.2 + 0.3*z[i] + 0.4*codes[i] + 0.5*randn(rng) for i in 1:90]
x[7:7:84] .= missing
y[14] = missing
data = (; y, x, z)
form = bf(@formula(y ~ z + mi(x)), @formula(sigma ~ 1))
predictor = impute_model(@formula(x ~ z);
    family = CumulativeLogit(), levels = ["low", "medium", "high"])
fit = drm(form, Gaussian(); data,
    impute = (x = predictor,),
    missing = miss_control(response = "include", predictor = "model"))
@assert isfinite(loglik(fit))
(cutpoints(fit), imputed(fit; rows = :missing))

([-0.6448749468671564, 0.7663893239060184], (variable = ["x", "x", "x", "x", "x", "x", "x", "x", "x", "x", "x", "x"], original_row = [7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84], model_row = [7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84], observed = Bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], estimate = [1.5878066037624508, 1.8476362104006432, 2.099280036744192, 1.602459056043538, 1.27093062431851, 2.8087306740147815, 1.9781840546158644, 2.192020405389711, 2.2789818524991605, 1.5960820897921224, 1.6817631106546826, 2.6331273987237713], std_error = Union{Missing, Float64}[0.6956960016045987, 0.804221353284999, 0.7691884654524983, 0.7043458301730751, 0.5108139243443377, 0.4349392010420196, 0.7717706083404391, 0.7569741326725772, 0.7392093490159926, 0.6985809058286164, 0.72669024160341, 0.5843048184212348], source = ["conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score", "conditional_expected_score"], uncertainty_status = ["ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok", "ok"]))

Use CategoricalLogit() for a nominal predictor. In direct Julia, coef(fit) and vcov(fit) retain the prepared raw parameter order, including ordinal cutpoint coordinates; cutpoints(fit) returns constrained cutpoints separately. This differs from R's public coefficient table, which omits predictor cutpoints. The fitted likelihood is shared, but full accessor parity and the remaining native default-fit discrepancies are still programme requirements. No-intercept means follow R's first-factor coding: the first categorical main effect receives full indicators, and later factors use reduced contrasts. An ordinal missing predictor uses polynomial contrasts unless it is that first factor. Complete numeric, plain-string, symbol and Boolean covariates are admitted. Generated R state-expanded designs verify plain-string and Boolean factors and their interactions, including coefficient names; symbol-valued factors have not yet been separately checked against R. Package-specific categorical/ordered value types require a typed contrast contract and are refused when used in the mean formula. They remain required parity work. Interactions containing mi() remain unsupported. R-prepared bridge designs retain R's coding.

For new rows, predict(fit, newdata) uses the fitted design and its factor contrasts. Supply a known predictor state; missing or unknown states are refused. The prediction does not condition on a supplied new response. Training fitted values remain available with fitted(fit).

newdata = (; z = [0.0, 0.5], x = ["low", "high"])
predict(fit, newdata)
# Residual SD needs only the columns used by the sigma formula.
predict(fit, (; z = [0.0, 0.5]); dpar = :sigma)

2-element Vector{Float64}:
 0.5279078289955729
 0.5279078289955729

Use type = :link for the linear predictor or type = :response for the response scale (the default). se = true returns predictions and delta-method standard errors when the retained observed-information covariance is finite, symmetric and positive definite. It refuses unavailable or invalid covariance; point predictions with se = false remain available. This is not a prediction interval, profile interval, or bootstrap. General missing-state new-data integration and complete R accessor parity remain separate requirements.

DRM.PreparedFiniteJointModel Type

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

source

DRM.PreparedFiniteJointFit Type

Fitted exact prepared ordinal or categorical missing-predictor model.

source

DRM.JointFiniteMissingMetadata Type

Row-level posterior state summaries for a finite-state prepared fit.

source

DRM._finite_joint_ordinal_logprobabilities Function

Log ordinal probabilities from c1, logspacing..., without sigmoid subtraction.

source

Finite-state retained-prediction internals (no stability guarantee)

DRM._joint_finite_state_prediction_plan Function

_joint_finite_native_state_design(formula, data, variable, levels, predictor)

Construct the row-then-state mean design by evaluating the complete mean formula after replacing its bare mi(variable) marker with each declared state. The single expanded StatsModels schema chooses no-intercept factor coding for all terms together: the first categorical term is full rank and later terms use treatment contrasts. Ordinal marker columns are then replaced by the fixed polynomial contrasts unless that marker itself is the full-rank first factor.

source

DRM.PreparedJointModel Type

PreparedJointModel

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

source

DRM.PreparedJointFit Type

PreparedJointFit

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

source

DRM.prepared_joint_model Function

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

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

source

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

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

source

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

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

source

DRM.prepared_joint_rowloglik Function

prepared_joint_rowloglik(model, theta)

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

source

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

source

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

source

DRM.prepared_joint_conditional_moments Function

prepared_joint_conditional_moments(model, theta)

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

source

prepared_joint_conditional_moments(model::PreparedTwoJointGaussianModel, theta)

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

source

prepared_joint_conditional_moments(model, theta)

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

source

DRM.fit_prepared_joint Function

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

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

source

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

source

DRM.joint_missing_summary Function

joint_missing_summary(fit::PreparedJointFit)

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

source

DRM.JointMissingMetadata Type

Original-row and conditional-moment metadata for a prepared joint fit.

source

DRM.PreparedJointGaussian Type

Tag for the prepared Gaussian-predictor prototype.

source

DRM.PreparedJointBernoulli Type

Tag for the prepared Bernoulli-predictor prototype.

source

DRM.prepared_joint_initial Function

A finite, neutral initial parameter vector in the documented fixed order.

source

DRM._has_joint_mi Function

Whether any formula axis contains the mi() marker.

source

DRM._fit_joint_formula Function

_fit_joint_formula(f, data; impute, missing, ...)

Build a prepared exact joint model from one Gaussian response formula containing one or two bare additive mi(x) mean terms. This deliberately partial frontend admits fixed Gaussian ML responses with complete remaining covariates. The one-predictor route permits Gaussian or Bernoulli predictor models; the two-predictor route requires two independent Gaussian predictor models. It neither drops rows nor fills modelled predictors before fitting.

source

DRM.drm_bridge_joint Function

drm_bridge_joint(payload)

Fit the shared prepared Gaussian-response joint likelihood from a versioned primitive payload. R owns formula parsing and exogenous-design validation. Observation masks, not working placeholder values, identify missing data. The result retains native design-column order, raw covariance coordinates, original row IDs and conditional imputation uncertainty. This internal bridge does not add new family, weighting, REML, profile or bootstrap admissions.

source

prepared_joint_nll

prepared_joint_nll(model, theta) returns the negative sum of row log-likelihoods, retaining predictor-only observations and the exact zero contribution of rows where both response and predictor are missing.