drmSEM overview

This vignette is the map of drmSEM: it ties together the engine/layer split, component-labelled paths, the simulation-based effect calculus, d-separation, and phylogenetic distributional SEM, and closes with a roadmap. It is deliberately prose-first; each section points to a companion vignette where the topic is worked in full.

• Distributional piecewise SEM with drmSEM (drmSEM.Rmd) — the end-to-end walkthrough.
• Effect decomposition (effect-decomposition.Rmd) — the simulation calculus in depth.
• Phylogenetic SEM (phylogenetic-sem.Rmd) — shared ancestry inside a distributional SEM.
• Comparison with other tools (comparison.Rmd) — where drmSEM sits in the landscape.
• Calibrating the d-separation test (calibration.Rmd) — cached evidence for the any-component d-separation grid.
• Effect-interval coverage and model-selection recovery (validation.Rmd) — the precomputed report scaffold for the next simulation grid.

The one idea

Most structural equation modelling treats a path X -> Y as a statement about the expected value of Y. But a predictor can do more than move a mean. It can change how variable a response is, how over-dispersed it is, how often it is a structural zero, how much among-group heterogeneity it carries, or how strongly two responses co-vary in their residuals. In ecology and evolution these are real, distinct hypotheses:

• “Warmer sites make body size more variable” — a claim about sigma, not mu.
• “Habitat raises the probability of absence, not the mean count” — a claim about zi.
• “Denser populations reduce among-site heterogeneity in survival” — a claim about sd(site).

drmSEM makes those targets first-class causal claims, keeps them labelled as such everywhere, and propagates their downstream consequences honestly.

The engine / layer split

drmSEM is a layer, not a fitting engine.

• drmTMB is the engine: it fits one distributional regression — mean, scale, shape, zero-inflation, hurdle, random-effect scale, residual correlation — per response.
• drmSEM is the graph / SEM / d-separation / path / effect-decomposition layer. It never fits its own likelihoods.

The consequences are concrete: each endogenous node is exactly one drmTMB fit, the system is piecewise (a graph of separately fitted local models, in the piecewiseSEM tradition of Shipley (2009) and Lefcheck (2016)), and the directed graph must be a DAG — cycles are an error. Every assumption about the engine’s return shapes lives in a single adapter, so drmSEM re-uses, never re-implements, drmTMB’s likelihoods, links, and family samplers.

There are two ways in, producing the identical object:

# Declarative: one drm_node() per endogenous response, a shared data frame.
sem <- drm_sem(
  size      = drm_node(drmTMB::bf(size ~ temp + habitat, sigma ~ temp),
                       family = stats::gaussian()),
  abundance = drm_node(drmTMB::bf(abundance ~ size + temp, zi ~ habitat),
                       family = drmTMB::nbinom2()),
  survival  = drm_node(drmTMB::bf(cbind(alive, dead) ~ abundance + size),
                       family = drmTMB::beta_binomial()),
  data = dat
)

# Assembler: if you have already fitted the nodes yourself, drm_psem() wires
# existing drmTMB objects into the same drm_sem object — no refitting.

This canonical chain, size -> abundance -> survival, recurs across the vignettes. Temperature is exogenous and reaches survival only indirectly, through the mean, the scale, and the zero-inflation of the intermediate nodes — exactly the situation that breaks coefficient-product mediation.

Component-labelled paths

A drmSEM edge is the tuple (from, to, component, link, term). The component is the distributional parameter of the target node that the predictor acts on:

Component	Symbol	What a path to it means
mean	mu	shifts the expected response
scale	sigma	changes residual scale / dispersion, not the mean
shape	nu	changes a shape parameter (e.g. tail/skew)
zero-infl.	zi	changes the probability of a structural zero
hurdle	hu	changes the probability of crossing the hurdle
RE scale	sd(group)	changes among-group heterogeneity
residual r	rho12	changes the residual correlation between the two responses of a bivariate node

paths() returns one row per fitted fixed-effect coefficient across all nodes, each tagged with its component and link. In the canonical model the temp path into size appears twice — once on mu, once on sigma — and they are not the same claim. The cardinal rule, enforced in paths(), plot(), the d-separation test, and the effect decomposition, is that a path to a non-mu component is never reported as a mean effect.

The rho12 component is the residual correlation between the two responses of a bivariate node (eps_y1 <-> eps_y2): a path rho12 ~ x says x changes that residual coupling, not the mean and not a directed y1 -> y2 arrow. Targeting rho12 today requires a bivariate drmTMB fit (e.g. family = c(gaussian(), gaussian()), components mu1/mu2/sigma1/sigma2/rho12) supplied via drm_psem(), after which the rho12 path is extracted and labelled like any other component path. The covariance-edge grammar now ships: covary() declares a residual or higher-level covariance edge, covariances() reports those edges separately from paths(), and d-separation knows two covarying responses should not get a y1 _||_ y2 claim. The bivariate-node declaration grammar now also ships: drm_pair() declares a joint two-response node, rho12() / corpairs() report the declared edges, and plot() draws double-headed (residual) / dashed (higher-level) arcs. Only the joint bivariate fit that reads a fitted rho12 / corpair correlation back (a non-NA estimate) remains a roadmap item (see docs/design/07-bivariate-covariance-edges.md).

paths(sem)   # component-labelled path table
plot(sem)    # the DAG, with sigma / zi edges drawn distinctly from mu edges

The effect calculus: counterfactual, not coefficient products

The classical “indirect = product of path coefficients” identity is exact only on a single linear, identity-link scale. drmSEM paths cross links (a log-link mean feeding a logit-link mean downstream) and flow through non-mean components, where a product of coefficients has no defined response-scale meaning. So drmSEM simulates the intervention instead: it sets the predictor to two contrasting values, pushes the change through every downstream node in topological order (random effects held at zero), and reads off the conditional (RE = 0) response-scale change in the target. This is counterfactual, do()-style propagation in the tradition of Pearl (2001), Pearl (2009) and Imai et al. (2010). By default the direct quantity is a controlled direct effect (Robins and Greenland 1992), and the total/indirect split is simulation-based; these coincide with natural direct and indirect effects only under linearity with no exposure–mediator interaction (VanderWeele 2014). When they diverge, indirect_effects(effect = "natural") returns the cross-world natural decomposition (natural_direct, natural_indirect, mediated_interaction) instead.

Mean vs distribution mediation. What an active mediator passes downstream is the crux. Under total_effects(method = "gcomp") it passes its expected mu; under method = "simulate" it passes a realized draw from its fitted family, using all of its response-scale parameters (sigma, zi, nu, trials, …). Only the latter carries signal that lives in a mediator’s scale, zeros, or shape.

The decomposition. indirect_effects() reports five quantities:

• total_path — the full simulated effect with mediators active;
• direct — the controlled direct effect (mediators held inactive);
• indirect — total_path - direct;
• mean_mediated — the part carried by mediator means;
• distribution_mediated — the extra effect that appears only when mediators pass realized draws, i.e. signal flowing through scale, zero-inflation, or shape (distribution total minus mean total).

So indirect is approximately mean_mediated + distribution_mediated, and the distribution_mediated row isolates exactly the contribution that no coefficient-product method can express. The default direct is a controlled direct effect, which coincides with the natural direct and indirect effects only under linearity with no exposure-mediator interaction. For the general case, indirect_effects(effect = "natural") returns the cross-world natural decomposition — natural_direct, natural_indirect, and mediated_interaction (Pearl 2001; Imai et al. 2010; VanderWeele 2014) — validated against the linear-Gaussian recovery case.

direct_effects(sem,   from = "temp", to = "survival")
total_effects(sem,    from = "temp", to = "survival", method = "simulate")
indirect_effects(sem, from = "temp", to = "survival")  # the five-row decomposition
indirect_effects(sem, from = "temp", to = "survival", effect = "natural")  # NDE / NIE / mediated interaction

Note: direct_effects() reports the mean (response) scale only. Non-mean channels (sigma, zi, …) do not appear here — they show up in paths() and in the distribution_mediated row of indirect_effects().

Outcome functionals. Because effects come from simulated populations, the target need not be a mean. total_effects(target = ...) reports the distribution-mediated effect on a response-scale functional — Pr(Y > t) ("p_gt"), the zero probability Pr(Y = 0) ("p_zero"), or the variance of Y ("var") — which is often the scientifically relevant quantity and is the headline estimand of distributional SEM. These are recovery-tested (for example, the p_zero effect recovers the Poisson zero-probability change). Coefficient uncertainty is propagated by drawing each node’s coefficients from MVN(coef, vcov) and repeating the whole propagation. See effect-decomposition.Rmd for the full treatment.

total_effects(sem, from = "temp", to = "survival", target = "p_gt", threshold = 0)
total_effects(sem, from = "temp", to = "survival", target = "p_zero")
total_effects(sem, from = "temp", to = "survival", target = "var")

d-separation: an any-component test

A piecewise SEM is tested by the independence claims its DAG implies (Pearl 2009): every pair of non-adjacent nodes should be conditionally independent given their parents. drmSEM sharpens the usual mean-only version (Shipley 2000, 2009; Lefcheck 2016). Because each node models a whole conditional distribution, a missing arrow X -> Y asserts that X predicts none of Y’s modelled components — mu, sigma, zi, and the rest. Adjacency is correspondingly any-component: an arrow into zi(Y) alone makes X adjacent to Y.

basis_set() enumerates the claims. dsep() tests each by refitting Y’s node with X added to every component sub-formula and comparing to the base node by a likelihood-ratio test, with degrees of freedom equal to the number of added coefficients. fisher_c() combines the claim p-values into Fisher’s C, C = -2 * sum(log p), on 2k degrees of freedom for k claims (Shipley 2000, 2009); a small p-value flags a missing path — including a missing path into a non-mean component, which a mean-only tool could never detect.

basis_set(sem)   # the implied independence claims
dsep(sem)        # any-component LRT per claim
fisher_c(sem)    # combined goodness-of-fit

The any-component definition is a deliberate, novel research choice: it treats “X is irrelevant to Y” as a statement about Y’s whole conditional distribution. The calibration grid (mean-only, distributional, and cross-link DGPs; n from 100 to 1000; 200 replicates per cell) met the pre-specified Type-I, augmented-component-count, Fisher-uniformity, and power criteria. Broader Fisher’s C and d-separation settings remain claim-scoped until separately tested. Claims whose augmented refit cannot be formed are flagged and dropped from Fisher’s C rather than silently scored.

Phylogenetic distributional SEM

Shared ancestry enters a drmSEM model as a structured random effect on a node. drmTMB exports the markers phylo(), animal(), relmat(), and spatial(), and drmSEM recognises them as markers: they are stripped from the causal edge set, so a node can carry a phylogenetic random effect today while paths(), dsep(), fisher_c(), and the effect engine operate on the fixed-effect DAG.

sem_phylo <- drm_sem(
  y1 = drm_node(drmTMB::bf(y1 ~ x + phylo(1 | species, tree = phy)),
                family = stats::gaussian()),
  y2 = drm_node(drmTMB::bf(y2 ~ y1 + x + phylo(1 | species, tree = phy)),
                family = stats::gaussian()),
  data = dat
)
paths(sem_phylo)
indirect_effects(sem_phylo, from = "x", to = "y2")

This is the genuinely new contribution. Existing phylogenetic SEM tools (phylopath (Bijl 2018), phylosem (Thorson and Bijl 2023)) focus on causal relationships among trait means under phylogenetic covariance; drmSEM extends this to distributional SEM, where shared ancestry can be incorporated while direct and indirect effects are estimated for changes in means, variances, zero probabilities, and other response-distribution components. Practical notes: phy is an ape “phylo” object and must be ultrametric (rescale a raw tree with ape::compute.brlen(phy, "Grafen")), and every species level must be a tip label. Phylopath-style model-set comparison (drm_dag() / drm_model_set() / compare() / best() / average(), Phase 2) and evolutionary covariance models (drm_phylo_cov() for BM / lambda / OU / kappa on a fixed grid, Phase 3) ship today; phylogenetic paths into distributional components and jointly estimated evolutionary covariance remain roadmap items. See phylogenetic-sem.Rmd.

Roadmap

Capability	Status
Component-labelled paths (mu, sigma, nu, zi, hu, sd(group), rho12)	now (a x -> rho12 path requires a bivariate drmTMB fit via drm_psem())
Bivariate covariance-edge grammar: covary() declaration + covariances() accessor + covariance-aware d-separation (drops the y1 _||_ y2 claim)	now
Bivariate-node declaration: drm_pair() grammar, rho12() / corpairs() accessors (declared edges), double-headed / dashed-arc plotting	now (0.2.x)
Bivariate joint fit: rho12() / corpairs() read-back of a fitted (non-NA) correlation from a live bivariate drmTMB fit	roadmap (0.4 engine)
Any-component d-separation + Fisher’s C	now (calibration grid validated; broader settings claim-scoped)
Simulation-based direct / mean-mediated / distribution-mediated / total effects	now
Natural direct / indirect / mediated-interaction effects (indirect_effects(effect = "natural"))	now (cross-world split, validated on the linear-Gaussian case)
Outcome functionals beyond the mean (Pr(Y>t), Pr(Y=0), var via direct_effects() / total_effects(target = ...))	now (first set; more functionals + bootstrap CIs on roadmap)
Unified effect-API surface (method, uncertainty, nsim, population shared across the effect functions)	now (mediation / draw / n_sim kept as deprecated aliases; bootstrap + marginal are roadmap)
Path-specific attribution: per-mediator (inclusion/exclusion), per-component (mean/sigma/zi), and natural variant with recanting-witness flag (path_effects())	now
Composite (formative / PCA) latent constructs with reliability + loadings (drm_composite() / loadings())	now
Component-labelled DAG plot + effect forest plot	now
Phylogenetic correction inside piecewise SEM (Phase 1, intercept-only)	now
Phylopath-style model-set comparison (best / averaged model)	now (Phase 2; drm_dag / drm_model_set / compare / best / average)
Evolutionary covariance models (BM / OU / Pagel’s lambda / kappa)	now (Phase 3; drm_phylo_cov, fixed grid — not jointly estimated)
Jointly estimated evolutionary covariance (λ / OU / κ)	roadmap
Phylogenetic paths into distributional components	roadmap (Phase 4)
Reflective latent measurement models (latent estimated jointly with the structural model)	out of scope for 0.x (needs a joint likelihood; composite constructs ship — see above)
Full joint (single-likelihood) multivariate SEM	out of scope for 0.x
Declared feedback motifs: drm_cycle() / feedback= grammar, relaxed toposort, d-sep suppression, equilibrium propagator	now (0.5.0)
Equilibrium total effect of a feedback SEM (total_effects(), fixed-point, NA on divergence)	now (0.5.x; the mean/distribution decomposition through a cycle stays out of scope)
Undeclared cyclic / feedback graphs	an error (declare the motif with drm_cycle() to allow it)
New likelihoods	out of scope (drmSEM never fits its own)

Recap

• drmSEM is a layer on the drmTMB engine: one fit per endogenous node, DAG-only, observed variables only.
• Paths are component-labelled; a non-mu path is never a mean effect.
• d-separation is the any-component LRT, combined by Fisher’s C.
• Effects are counterfactual and simulation-based, never coefficient products, and the indirect decomposition exposes the distribution-mediated channel explicitly.
• Phylogeny enters as a structured random effect, giving a first slice of phylogenetic distributional SEM.

For the full walkthrough see drmSEM.Rmd; for effects, effect-decomposition.Rmd; for phylogeny, phylogenetic-sem.Rmd; for positioning, comparison.Rmd.

References

Bijl, Wouter van der. 2018. “phylopath: Easy Phylogenetic Path Analysis in R.” PeerJ 6: e4718. https://doi.org/10.7717/peerj.4718.

Imai, Kosuke, Luke Keele, and Teppei Yamamoto. 2010. “Identification, Inference and Sensitivity Analysis for Causal Mediation Effects.” Statistical Science 25 (1): 51–71. https://doi.org/10.1214/10-STS321.

Lefcheck, Jonathan S. 2016. “piecewiseSEM: Piecewise Structural Equation Modelling in R for Ecology, Evolution, and Systematics.” Methods in Ecology and Evolution 7 (5): 573–79. https://doi.org/10.1111/2041-210X.12512.

Pearl, Judea. 2001. “Direct and Indirect Effects.” Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence (UAI), 411–20.

Pearl, Judea. 2009. Causality: Models, Reasoning, and Inference. 2nd ed. Cambridge University Press.

Robins, James M., and Sander Greenland. 1992. “Identifiability and Exchangeability for Direct and Indirect Effects.” Epidemiology 3 (2): 143–55. https://doi.org/10.1097/00001648-199203000-00013.

Shipley, Bill. 2000. “A New Inferential Test for Path Models Based on Directed Acyclic Graphs.” Structural Equation Modeling 7 (2): 206–18. https://doi.org/10.1207/S15328007SEM0702_4.

Shipley, Bill. 2009. “Confirmatory Path Analysis in a Generalized Multilevel Context.” Ecology 90 (2): 363–68. https://doi.org/10.1890/08-1034.1.

Thorson, James T., and Wouter van der Bijl. 2023. “phylosem: A Fast and Simple R Package for Phylogenetic Inference and Trait Imputation Using Phylogenetic Structural Equation Models.” Journal of Evolutionary Biology 36 (10): 1357–64. https://doi.org/10.1111/jeb.14234.

VanderWeele, Tyler J. 2014. “A Unification of Mediation and Interaction: A 4-Way Decomposition.” Epidemiology 25 (5): 749–61. https://doi.org/10.1097/EDE.0000000000000121.