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

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

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

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

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

Temporal random effects (AR1 / OU)

Wave 1 (Gaussian mean, sigma ~ 1, ML). The Julia spelling is positional because @formula cannot carry keyword arguments: drmTMB's temporal(1 | id, time = occ, structure = "ar1") is temporal(1 | id, occ, ar1) here (the R bridge accepts drmTMB's spelling).

fitted and predict are population-level (Xβ̂; drmTMB's fitted() is conditional, the conditional temporal effects are ranef(fit)[:id]). simulate and bootstrap_ci draw a fresh temporal chain per series (and a fresh (1 | id) intercept), as drmTMB's default simulate(); drmTMB refuses the temporal bootstrap, so bootstrap_ci is an extension. Wald standard errors are reported for every coordinate; drmTMB exposes only AR1 mean-coefficient Wald intervals, and no interval calibration is claimed.

DRModels.temporal Function

temporal(1 | id, time, structure)

Temporal random-intercept marker on the Gaussian mean: a stationary AR1 (ar1) or Ornstein–Uhlenbeck (ou) process within each series id, indexed by the column time, or a homogeneous Toeplitz (homtoep) within-series covariance on a complete, equally spaced panel.

drm(bf(@formula(y ~ x + temporal(1 | id, occ, ar1)), @formula(sigma ~ 1)),
    Gaussian(); data)
drm(bf(@formula(y ~ x + temporal(1 | id, elapsed, ou)), @formula(sigma ~ 1)),
    Gaussian(); data)

This is the twin of drmTMB's temporal(1 | id, time = occ, structure = "ar1"). StatsModels' @formula cannot parse keyword arguments or string literals, so the Julia spelling is positional with a bare structure name (the same convention as sd(g, phylogenetic)): temporal(1 | id, occ, ar1). Through the R bridge (drm_bridge) the exact drmTMB spelling is accepted and translated.

Semantics (as drmTMB): the first state of each series is stationary N(0, 1); AR1 needs integer time and keeps real gaps as integer powers, φ^gap (φ = tanh θ, may be negative); OU takes numeric elapsed time with correlation exp(−λ Δt) over a gap Δt (λ = exp θ > 0, in the units of time). Only gaps matter, so the time origin is irrelevant. Each (id, time) pair must be unique; rows may be in any order. The fitted quantities are the process SD (re_sd(fit)[:id]), the persistence φ or decay λ (temporal_parameters), the residual SD σ, and, when the formula also has (1 | id), the stable intercept SD (re_sd(fit)[:id_iid]).

Post-fit conventions (some differ from drmTMB):

• fitted(fit) and predict(fit, newdata) are POPULATION-level, Xβ̂ (the DRModels convention for every structured term); drmTMB's fitted() adds the conditional temporal effects. Those are in ranef(fit)[:id], data-row order (and ranef(fit)[:id_iid] per series for the ordinary intercept).

• simulate(fit) draws from the fitted marginal model: a fresh stationary chain per series (and a fresh (1 | id) intercept) plus residual noise, as drmTMB's default simulate(). bootstrap_ci uses the same draws; drmTMB refuses the temporal bootstrap, so these percentile intervals are a DRModels.jl extension with no calibration claim.

• vcov/stderror/Wald confint cover every coordinate (observed Hessian), as on DRModels' other routes; drmTMB exposes only AR1 mean-coefficient Wald intervals because the calibration of the others is not established. No calibration is claimed here either.

Wave-1 scope: Gaussian family, sigma ~ 1, ML only, intercept-only and unlabelled, at most one ordinary (1 | id) on the same id. Every other use is refused with an error.

Phylogenetic stable intercept + OU (wave 2). drmTMB's paired provider phylo(1 | species, tree = tree) + temporal(1 | species, time = elapsed, structure = "ou") is spelled

drm(bf(@formula(y ~ x + phylo(1 | species) + temporal(1 | species, elapsed, ou)),
       @formula(sigma ~ 1)), Gaussian(); data, tree = newick)

and fits y = Xβ + a_species + b_species(t) + ε with a ~ N(0, σ_a² C) (C the tree's tip correlation matrix, the scale drmTMB uses) and an independent OU path b per species. Related species share a stable baseline; their temporal departures are independent. The stable SD is re_sd(fit)[:species_phylo] (temporal_parameters(fit).sd_phylo; drmTMB sd_phylo_stable), the process SD re_sd(fit)[:species] (drmTMB sd_temporal) and the rate temporal_parameters(fit).decay (drmTMB decay_temporal). ranef(fit)[:species_phylo] holds the per-species stable modes (first-seen order) and ranef(fit)[:species] the temporal modes (data rows). As in drmTMB the pairing needs ou, an unlabelled intercept-only phylo() on the same grouping, no ordinary (1 | species), at least three species with at least two times each, at least three distinct positive lags, and tree tips that are exactly the observed species. simulate draws a fresh phylogenetic vector, fresh OU paths and fresh noise.

Homogeneous Toeplitz (wave 2). temporal(1 | id, occ, homtoep) (drmTMB structure = "homtoep") needs every id observed at the same 3–12 equally spaced integer occasions. Each series' responses then have covariance σ² R, R a free positive-definite Toeplitz correlation matrix (one correlation per lag, temporal_parameters(fit).cor, drmTMB cor_lag1…), estimated through its partial autocorrelations (coef(fit, :temporal_pac), atanh scale). sigma is the TOTAL within-series SD: as in drmTMB there is no separate process SD, residual SD or (1 | id) (not identified when every lag is free), and no latent states (ranef(fit) is empty). As in drmTMB, Wald covariance is withheld, and profile intervals and curves (confint, profile_curve) are refused for sigma and the PACs; parameter_surface, which has no drmTMB counterpart, keeps the same mean-only scope. drmTMB refuses the temporal bootstrap; as an extension, bootstrap_ci / bootstrap_summary / bootstrap_result report percentile intervals for the mean coefficients only, with no calibration claim. Rows with a missing response are dropped before the panel rules apply (as drmTMB); simulate then returns one value per data row, NaN at the dropped rows. Incomplete or unequally spaced panels, fewer than 3 or more than 12 occasions, fractional occasions and an ordinary (1 | id) are refused with drmTMB's messages.

source

DRModels.temporal_parameters Function

temporal_parameters(fit) -> NamedTuple

The fitted temporal process of a temporal(...) Gaussian fit, on natural scales: label (drmTMB's term label), structure (:ar1 / :ou), sd (process SD σ_t; drmTMB sdparsmu["temporal_sd: <label>"]),phi(AR1 persistence, one occasion apart; drmTMBcorpars$temporal) or decay (OU rate λ in the units of time; drmTMB decaypars$temporal) — the other isnothing—sd_iid(the ordinary(1 | id)SD, ornothing),sd_phylo(the stable phylogenetic SD of a pairedphylo(1 | species) + temporal(…, ou)fit, on the tip-correlation scale; drmTMBsdpars$mu["sd_phylo_stable"], or nothing) and sigma (residual SD). For the paired fit drmTMB names the process SD sd_temporal and the rate decay_temporal. For a homogeneous Toeplitz fit (homtoep) sigma is the TOTAL within-series SD, sd, phi, decay, sd_iid and sd_phylo are nothing, cor holds the lag correlations ρ_1…ρ_{K−1} (drmTMB corpars$temporal,cor_lag1…) andpacthe partial autocorrelations; both arenothing for the other structures. Errors on a fit without a temporal term.

source

Known sampling variance (meta-analysis)

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

Random intercepts on the mean combine with meta_V, e.g. bf(y ~ x + meta_V(v) + (1 | study), sigma ~ 1) or bf(y ~ x + meta_V(v) + phylo(1 | species) + (1 | study), sigma ~ 1) (with tree = …); relmat / animal take K = … / A = …. The marginal is N(Xβ, diag(v + σ²) + Σₖ sₖ² Zₖ Cₖ Zₖᵀ), drmTMB's model. Each random component needs its own grouping column; slopes, spatial, and REML are not implemented.

source

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

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

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

DRModels.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 a supported model; use the capability page to check the combinations available for your response family.

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

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

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

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

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

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

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