Structured-effect markers

Status — Reference

Mirrors drmTMB's Structured-effect markers (6 in drmTMB). These markers wrap a random-effect term inside a bf formula to give it a known correlation structure (phylogeny, space, pedigree, an arbitrary relatedness matrix) or a known sampling-variance (meta-analysis).

Correlation-structured random effects

DRM.phylo Function

phylo(1 | species)

Phylogenetic structured random intercept on the Gaussian mean: pass the tree via drm(...; tree = tree) (an AugmentedPhy from random_balanced_tree / augmented_phy, or a Newick string). The phylogenetic correlation is built from the tree (sigma_phy_dense) and the marginal is fit in closed form. (The q=4 phylogenetic location-scale model — a structured effect on log σ too — uses the verified sparse-Laplace engine instead; see HANDOVER.md.)

source

DRM.spatial Function

spatial(1 | site)

Coordinate-spatial structured random intercept on the Gaussian mean. Pass site coordinates via drm(...; coords = coords) (a G×2 matrix, one row per site level in first-seen order). The spatial correlation K(ρ) = exp(-d / ρ) is built from pairwise distances and the range ρ is estimated jointly. Closed-form Gaussian marginal (K is rebuilt each evaluation since it depends on ρ).

source

DRM.animal Function

animal(1 | id)

Animal-model structured random intercept. Supply the additive-relatedness matrix over the levels of id via drm(...; A = A). Reuses the closed-form structured-Gaussian engine (same as relmat).

source

DRM.relmat Function

relmat(1 | id)

Structured random-intercept marker with a user-supplied relatedness matrix:

drm(bf(y ~ x + relmat(1 | id), sigma ~ 1), Gaussian(); data, K = K)

K is the correlation/relatedness matrix over the levels of id, ordered as they first appear in data. The marginal stays Gaussian (closed-form).

source

Known sampling variance (meta-analysis)

DRM.meta_V Function

meta_V(v)

Formula marker for Gaussian meta-analysis: v is the data column of known sampling variances. Use inside a μ formula, e.g. bf(y ~ x + meta_V(v), sigma ~ 1). The residual (between-study) heterogeneity SD is the σ parameter (the between-study τ of classical meta-analysis). (Diagonal known variances; dense/bivariate sampling covariance is planned.)

source

DRM.meta_vcov_bivariate Function

meta_vcov_bivariate(v1, v2; cov12 = nothing, cor12 = nothing) -> MetaVcovBivariate

Build the known row-paired sampling covariance for a bivariate meta-analysis — drmTMB's meta_vcov_bivariate(). One entry per study: sampling variances v1, v2 and either sampling covariances cov12 or sampling correlations cor12 (scalar or vector; at most one of the two, default independence).

Consumed by the bivariate Gaussian route as a fit-call keyword:

V = meta_vcov_bivariate(v1, v2; cor12 = 0.6)
fit = drm(bf(mu1 = @formula(y1 ~ x), mu2 = @formula(y2 ~ x),
             sigma1 = @formula(sigma1 ~ 1), sigma2 = @formula(sigma2 ~ 1),
             rho12 = @formula(rho12 ~ 1)),
          Gaussian(); data = dat, V = V)

The fitted sigma1 / sigma2 are the between-study heterogeneity SDs and rho12 the residual (heterogeneity) correlation — the known sampling covariance is added on top per row, never absorbed into them. (This mirrors the univariate meta_V contract: sigma is heterogeneity alone, V separate.)

drmTMB spelling

drmTMB writes bf(mu1 = y1 ~ x + meta_V(V = V), …). A keyword argument inside a formula is not representable in StatsModels' @formula, so DRM.jl takes the object via the V = keyword instead — the same place tree / K / A live. Writing meta_V(...) inside a bivariate formula errors with a pointer to this keyword.

source

DRM.MetaVcovBivariate Type

MetaVcovBivariate

Known row-paired sampling covariance for bivariate meta-analysis — the object meta_vcov_bivariate builds and drm(...; V = …) consumes. Fields v1, v2, cov12, one entry per study.

Matrix(V) materialises drmTMB's dense 2n × 2n block-diagonal form (stacking order y1[1], y2[1], y1[2], …); the constructor also accepts such a matrix back, refusing any cross-study (off-block) entry.

source

Location–scale–scale submodel markers

DRM.sd Function

sd(group)

Formula marker for a location–scale–scale model. Used only on the left-hand side of a bf component formula,

bf(y ~ x + (1 | g),              sigma ~ x, sd(g) ~ z)                  # iid RE SD
bf(y ~ x + phylo(1 | species),   sigma ~ x, sd(species, phylogenetic) ~ z)  # phylo SD

to put a linear predictor on the log standard deviation of a random effect. Plain sd(g) targets the iid (1 | g) intercept (log σ_b,k = Z_k' α, one row per group level); sd(group, phylogenetic) targets the per-species phylogenetic SD (drmTMB: sd(group, level = "phylogenetic"); @formula cannot parse keyword arguments, so the level is a bare symbol). Predictors must be constant within each group level. Univariate Gaussian only; coefficients appear as the :sd (iid) or :sd_phylo (phylogenetic) block.

Capabilities

• Estimators: supports both Maximum Likelihood (method = :ML, default) and Restricted Maximum Likelihood (method = :REML).

• Missing response: incomplete responses (missing or NaN in y) are fit via the observed-rows pattern (response = "include" in R bridge), preserving group and phylogenetic structures across all G levels.

• Sparse scaling: phylogenetic LSS models scale to large trees in O(p) time via sparse = true (or algorithm = :sparse_lbfgs), automatically selected when G > 500 species.

• Multi-component models: combines multiple grouping factors or iid + phylogenetic random-effect SD models.

source

DRM.sd_phylo Function

sd_phylo(group)

DEPRECATED legacy spelling for the phylogenetic SD submodel — use sd(group, phylogenetic) ~ … (drmTMB: sd(group, level = "phylogenetic")). Both put a linear predictor on the log per-species SD of the phylogenetic random effect: a ~ MVN(0, D_a K D_a) with D_a = Diagonal(exp.(Z * α)) and K the Brownian-motion phylogenetic correlation from tree. Kept working, exactly as the twin keeps its legacy spelling, and canonicalised to the same :sd_phylo block; new code should write the sd(group, …) form.

source

Advanced tree preparation helpers

These exported helpers prepare or inspect phylogenetic covariance inputs for advanced workflows. They do not make every tree or structured-effect combination an admitted fitting route; use the capability matrix for that boundary.

DRM.augmented_phy Function

augmented_phy(newick::AbstractString) :: AugmentedPhy{Float64}

Parse a minimal Newick string and return the augmented-state sparse precision representation.

Restrictions

• Rooted multifurcating trees are admitted: every internal node must have at least two children. Unary nodes are currently unsupported.

• Leaf names may be old-compatible unquoted labels or lossless single-quoted labels. Doubled apostrophes inside a quoted label represent one apostrophe. Internal labels are tolerated but discarded.

• Literal NUL bytes are invalid Newick input and are rejected before parsing.

• Non-root branch lengths must be finite, > 0, and have finite reciprocal.

• The root has no parent branch; the optional root length in (…):0.0; is read but does not enter Q.

Example

phy = augmented_phy("((A:0.1,B:0.2):0.3,C:0.5);")
phy.n_leaves      # 3
phy.n_total       # 5  (3 leaves + 2 internal)
length(phy.branch_lengths)   # 4
nnz(phy.Q_topology)          # 13  (5 diagonal + 8 off-diagonal entries)

source

DRM.random_balanced_tree Function

random_balanced_tree(p::Integer; branch_length::Real = 0.1) :: AugmentedPhy

Build a near-balanced binary tree with p leaves. All branch lengths equal branch_length. Used in benchmarks and scaling tests.

When p is a power of 2 this is perfectly balanced. Otherwise the left-over leaf at each level is carried up one extra step (so the tree remains binary, just with slightly uneven depths). Branch lengths stay uniform — the goal is a representative sparse-tree topology, not an ultrametric one.

source

DRM.random_caterpillar_tree Function

random_caterpillar_tree(p::Integer; branch_length::Real = 0.1) :: AugmentedPhy

Build a caterpillar (maximally unbalanced "ladder") binary tree with p leaves: leaves 1 and 2 form the first cherry, that node joins leaf 3 at the next internal node, and so on — a chain of p-1 internal nodes of depth p-1.

The shape is the worst case for sparse-Cholesky fill-in: O(p) on a balanced tree (log-depth) does not by itself prove O(p) on a deep caterpillar, so this generator feeds the multi-shape scaling sweep (#16). All branch lengths equal branch_length.

source

DRM.phylo_tree_height Function

phylo_tree_height(phy::AugmentedPhy) -> Float64

Maximum root-to-tip path length of phy. For an ultrametric tree this is the common tip variance its covariance implies when σ²_phy = 1; for a non-ultrametric tree it is the largest tip variance.

O(p): a breadth-first walk over the sparse topology, where an edge's length is recovered as -1/Q_topology[i, j]. sigma_phy_dense would give the same number but inverts a dense matrix, so it is unusable as a routine check.

Why this matters. DRM.jl builds its phylogenetic covariance from the branch lengths as supplied. On an ultrametric tree of height h, the fitted sd_phylo carries a factor sqrt(h) relative to unit-tip-variance correlation scale. R's drmTMB instead standardises via ape::vcv(tree, corr = TRUE), whose tips always have variance 1. A non-ultrametric conversion is tip-wise, not one scalar. See drm_phylo_penalty, where these choices change what sd_u means.

For an ultrametric tree, divide all branch lengths by h before fitting to match drmTMB's unit-tip-variance correlation scale. A non-ultrametric tree needs the tip-wise standardization D^{-1/2} Σ D^{-1/2}, not one scalar branch rescaling; this raw-branch constructor deliberately performs neither transform.

source

DRM.augmented_tree_precision Function

augmented_tree_precision(phy::AugmentedPhy) -> (Q, leaf_pos, q)

Return the root-conditioned augmented topology precision Q = Q_topology[keep, keep] — a sparse, O(p)-nnz, positive-definite matrix over the q = n_total - 1 non-root augmented nodes — together with the map leaf_pos[t] from leaf t ∈ 1:p to its row/column in Q, and q. This is the sparse precision the end-to-end O(p) Gaussian path feeds DIRECTLY as Qₖ, bypassing the dense leaf-correlation inversion.

source

DRM.sigma_phy_dense Function

sigma_phy_dense(phy::AugmentedPhy; σ²_phy::Real = 1.0) :: Matrix

Build the dense (p × p) leaf covariance Σ_phy = σ²_phy · (S Q_cond⁻¹ S') where Q_cond is phy.Q_topology with the root row/col removed and S selects leaves. This is what the existing dense path expects; used by verification tests to compare sparse vs. dense.

This is O(p³) in storage and time — only intended for small trees in tests. Do NOT call it on the real workload; the entire point of AugmentedPhy is to avoid materialising Σ_phy.

source

DRM.phylo_correlation Function

phylo_correlation(tree) -> Matrix

Leaf correlation matrix from a phylogeny. tree is an AugmentedPhy (e.g. from random_balanced_tree / augmented_phy) or a Newick string. The phylo covariance is built with sigma_phy_dense and rescaled to unit diagonal.

source