Changelog

drmSEM (development version)

Outcome Functionals & Non-Gaussian Tail Exceedance (0.5.0)

Tail Risk and Quantile Analytics (target = c("mean", "p_gt", "p_zero", "var", "quantile")): direct_effects(), indirect_effects(), total_effects(), and path_effects() can now target outcome functionals beyond the conditional mean. Supports exceedance tail probabilities $\Pr(Y > t)$ (target = "p_gt", threshold = ...), structural zero/hurdle probabilities $\Pr(Y = 0)$ (target = "p_zero"), response variance $\text{Var}(Y)$ (target = "var"), and conditional quantiles $Q_p(Y)$ (target = "quantile", prob = 0.50).
Closed-Form Analytic Evaluations (functional = "analytic"): Analytical expectation, variance, zero probability, and exceedance functions implemented for gaussian, poisson, lognormal, Gamma, nbinom2, beta, and student families.
Cross-World Natural Mediation on Functionals (effect = "natural"): 4-way decomposition ( $\Delta_{\text{total}} = \text{NDE} + \text{NIE} + \text{INT}_{\text{med}}$ ) evaluated directly on tail risk and quantiles, supported by path_effects(by = "component") to attribute risk pathways to mediator mean vs dispersion/zero-inflation channels.

Hierarchical Uncertainty & Population Integration (0.5.0)

Cluster / Block Bootstrap Uncertainty (uncertainty = "bootstrap"): Non-parametric parameter uncertainty via cluster/block resample refitting (R = 500), preserving intra-cluster correlation structures across the full piecewise SEM.
Gauss-Hermite Marginal Population Integration (population = "marginal"): Integrates out random effects to evaluate population-averaged marginal effects across nonlinear link functions using adaptive Gauss-Hermite quadrature (identity/log links evaluate via exact closed forms).

Covariance Cliques & K >= 3 Partitioning (0.5.0)

Complete Covariance Clique Detection (covary_clique(), covariance_cliques()): Bron-Kerbosch maximal clique algorithm partitions residual and higher-level covariance networks into complete cliques vs structured graphs.
Basis-Set Independence Suppression: basis_set() and dsep() automatically suppress all $\binom{K}{2}$ within-clique pairwise independence claims for complete covariance blocks, avoiding spurious test rejections.

Comprehensive Vignette Suite & Documentation Architecture (0.5.0)

Expanded Latents & Measurement Guide (vignettes/latent-variables.Rmd): Comprehensive walkthrough of formative, reflective, and MIMIC measurement constructs in piecewise SEM. Demonstrates marker variable ( $\lambda_{\text{marker}} = 1$ ) vs unit-variance ( $\text{Var}(\eta) = 1$ ) identification, construct reliability calculation (Cronbach’s $\alpha$ vs Raykov’s composite $\rho$ ), downstream distributional effect propagation to location ( $\mu$ ), scale ( $\sigma$ ), and zero-inflation ( $\text{zi}$ ) targets, and the clean separation of loadings(sem) from structural paths(sem).
Expanded Feedback Cycles & Equilibrium Guide (vignettes/feedback-cycles.Rmd): Comprehensive guide to declaring reciprocal feedback loops with drm_cycle(). Details multi-component Banach fixed-point iteration across $\mu, \sigma, \nu, \text{zi}$ , stability and contraction diagnostics (spectral radius $\rho(B) < 1$ and empirical Lipschitz constants), equilibrium total effects via total_effects(), and honest limits regarding simultaneity bias in node-wise estimation.
New Graph-Derived Piecewise Imputation Guide (vignettes/missing-data.Rmd): Full walkthrough of drm_sem(impute = "auto"). Explains how missing mediator models are derived directly from the causal DAG across continuous (gaussian, Gamma, lognormal, student, beta) and count/discrete families (poisson, nbinom2, ZIP, beta_binomial), multi-parent imputation ( $k = 2$ ), diagnostic accessors (imputation(), imputed(), check_sem()), and honest piecewise variance bounds vs joint FIML.
Refined Bivariate Nodes & Moderation Guide (vignettes/bivariate-nodes.Rmd): Complete guide to joint two-response estimation with drm_pair(), distinguishing directed paths ( $y_1 \to y_2$ ), within-observation residual correlations ( $\rho_{12}$ on $\tanh$ link), and individual-level random-effect correlations (corpairs). Demonstrates environmental moderation of residual coupling ( $x \to \rho_{12}$ ) and explains basis-set suppression in d-separation.
Modernized _pkgdown.yml Navigation: Articles reorganized into five intuitive thematic sections: “Get Started”, “Distributional Pathways”, “Latents & Measurement”, “Feedback & Covariance”, and “Missing Data & Diagnostics”, passing pkgdown::check_pkgdown() with 0 missing topics and 0 unreferenced vignettes.

Joint bivariate residual correlation (0.4)

drm_pair() now fits one joint bivariate drmTMB model. Passing a pair to drm_sem() (or a bivariate fit to drm_psem()) no longer expands into two independent univariate nodes. Both responses share one biv_gaussian() / biv_lognormal() / biv_student() likelihood. rho12() returns the Wald table of rho12 coefficients — intercept and any rho12 ~ x predictors — with standard error and p-value on the tanh (engine: atanh_guarded) link. basis_set() drops the y1 _||_ y2 independence claim because the residual covariance edge is recorded automatically. This is still one drmTMB node, not a global multivariate SEM. Mixed-family pairs abort. Recovery: V-128 (constant ({12})), V-129 (({12} x)). See docs/design/07-bivariate-covariance-edges.md.

MAG m-separation under latent confounding

latent = on drm_sem() / drm_psem() now drives MAG m-separation. When the argument names marginalised latents, basis_set() and dsep() generate and test m-separation claims on the implied MAG instead of d-separation on the observed DAG. Conditioning is Richardson & Spirtes (2002) Corollary 5.3 anteriors (empty selection set); pairwise ⇒ global is licensed by Sadeghi & Lauritzen (2014) Theorem 3 and Lauritzen & Sadeghi (2018) Theorem 4 under a compositional-graphoid assumption. The any-component LRT is unchanged — only which claims are generated changes. Bounds, stated plainly: this is a graph projection of latent confounding onto a piecewise observed SEM, not FIML and not a reflective measurement model; selection / conditioned latents are structurally unrepresentable; parent-based Shipley & Douma separators stay unlicensed; residual OQ-16 is whether drmSEM’s any-component LRT / Fisher’s C induces a compositional graphoid. m-separation does not turn a path to sigma, zi, or nu into a mean effect. See docs/design/14-m-separation.md and D-20 (V-117 / V-118 / V-119).

Ordinal and spatial nodes

Ordinal (cumulative_logit) outcomes are now evidenced. They already worked — drmSEM imposes no family whitelist — but nothing tested them, which is worse than absent capability because the surface advertised it. paths() labels the edge mu/logit and hides the cutpoints, dsep() tests the claim, effects close additively, and a gaussian → ordinal chain recovers 0.5009 / 0.9330 against true 0.500 / 0.900. Three limitations ship as assertions, not surprises: mu is the latent linear predictor rather than E[category], so target = "mean" reports a latent-scale effect with no warning; a non-mean target such as p_gt returns exactly 0, because there is no realized-value sampler and the mean fallback collapses both scenarios; and an ordinal node cannot be distributional, so no path can target its scale.
Spatially-structured nodes are documented as relmat(1 | site, K = <matrix>). This already worked and was already tested — but only under a phylogenetic label, so it was unfindable by anyone asking about space. relmat(K = ) accepts any positive-definite matrix and is therefore strictly more flexible than drmTMB’s spatial() marker, whose only kernel is a fixed exponential with a heuristic, non-estimated range and whose mesh= is unimplemented.
Hurdle mediators no longer drop their zeros. A hurdle node reported sampler = TRUE and was then drawn as a plain truncated NB2: drmTMB folds the hurdle into model_type while family$family stays truncated_nbinom2, and drmSEM keyed on the name, so hu was never read and the entire zero component vanished from the propagated distribution. Sampling now resolves through model_type for an explicit allow-list, leaving zi_poisson/zi_nbinom2 on their existing (correct) path. Verified against drmTMB::simulate(): zero fraction 0.4444 vs 0.4446.
drm_psem() has tests. One of the package’s two documented interfaces, exported and named in the charter, previously had none: every SEM in the suite was built through drm_sem(). Its input validation and its default-data branch are now covered.

Distribution-mediated effects

The realized-value sampler list is no longer stale. drm_supported_sampler_families() listed ten families and its own comment asserted that “tweedie … has no realized-value sampler”. The engine had covered it all along: simulate.drmTMB draws rtweedie_compound(n, mu, phi = sigma^2, power = nu), which answers the exact sigma-to-dispersion question the stale comment cited as the blocker. tweedie, skew_normal, binomial and beta_binomial are now sampled, each admitted only after a moment test against drmTMB::simulate() (V-82..V-85). Consequence of the old state: a tweedie mediator silently degraded to mean propagation, losing the distribution-mediated channel that is the package’s headline claim — because of a stale vector, not an engine limitation.
binomial/beta_binomial mediators were a units bug, not just a missing feature. Both put a probability in mu while the response is a count, so the mean fallback handed downstream nodes a value on (0,1) where they had been fitted on counts. Both the sampler and drm_family_expected_mean() now use the fitted trials, so this affected mediation = "mean" as well. Where trials cannot be recovered, drmSEM warns and falls back rather than returning a probability.
The advisory list can no longer drift from the actual samplers. drm_supported_sampler_families() only feeds check_sem()’s sampler column; the load-bearing list is the switch() in drm_sample_family(). V-86 asserts every family named in the vector actually draws and that an unnamed one warns, so widening one without the other now fails a test instead of making check_sem() lie.
Parameterizations are not restated in drmSEM: the new branches call drmTMB’s own generators (rtweedie_compound, rskew_normal_public, drm_beta_shapes), so a mapping cannot drift between the two packages.

Missing data

drm_sem() now has a row-alignment policy. Previously drmSEM had none at all: one data frame went to every node and each drmTMB fit dropped incomplete rows by its own rules, so a piecewise SEM could be fitted on several different samples with nothing said. The new na_action argument takes "warn" (default — fit, but report which node used how many rows), "common" (fit every node on the shared complete-case set), or "fail". Per-node counts appear in attr(x, "alignment_issues"), in print(), and as an nobs column in check_sem(). "common" intersects the per-node row sets, so columns no node models cannot cost observations.
A silent wrong-number bug is fixed. drm_fixed_design() built its design matrix with model.matrix(), which honours na.action = na.omit, so with missing predictors it returned fewer rows than requested. That surfaced either as "number of items to replace is not a multiple of replacement length" or — when the row counts divided evenly — as silent recycling into a scrambled design matrix. The model frame now passes NAs through, so incomplete rows become NA predictions instead of wrong numbers. This also fixes standardize(method = "latent").
d-separation no longer compares two different samples. A claim whose augmented refit used a different number of observations (because the added variable has missing values) now reports status = "n_mismatch" and is excluded from Fisher’s C, instead of contributing an invalid likelihood ratio.
drm_sem(impute = "auto") derives each missing predictor’s imputation model from the causal graph. In x -> m -> y with m incomplete, node y needs a missing-predictor model for m — and that model is node m’s own formula and family. drmSEM emits mi(m) plus impute_model(m ~ x, family = <node m's family>) automatically; the user never writes an impute_model(). The conditioning set is therefore derived from the DAG rather than guessed by folklore. (A congeniality argument follows from imputer and analysis model sharing one graph, but nothing in the package measures congeniality and it carries no validation number — it is a design rationale, not a validated claim.) Inspect what was derived with imputation(); inspect fitted missing values with imputed(). Opt-in, because imputation asserts a missing-at-random assumption. Limit, stated plainly: the SEM is piecewise, so the downstream node re-estimates the parent’s model inside its own likelihood rather than sharing its estimates — imputation uncertainty is propagated within a node but not across nodes, and this is never FIML across the SEM. A node with two incomplete Gaussian parents emits two independent mi() terms (drmTMB 0.7.0); a Gamma response with one Bernoulli parent also emits (drmTMB #1088, mp-gamma-bernoulli, V-122); a lognormal response with one Bernoulli parent also emits (drmTMB #1092, mp-lognormal-bernoulli, V-123); a beta_binomial response with one Bernoulli parent also emits (drmTMB #1094, mp-beta-binomial-bernoulli, V-124); a student response with one Bernoulli parent also emits (drmTMB #1096 5fdf834c1, mp-student-bernoulli, V-125); an nbinom2 response with one Gaussian parent also emits (drmTMB #1095, mp-nbinom2-gaussian, V-126); a zero-inflated Poisson node (family = poisson() plus zi ~ 1) with one Bernoulli parent also emits (drmTMB #1097 9d4b63400, mp-zi-poisson-bernoulli, V-127; D-23: mi() in mu only, ZIP mixture, not a poisson alias). Leftover zi_nbinom2 and mi() on zi stay refused. k > 2, a non-Gaussian k = 2, and those leftovers still fail loud. This is not a general missing-data SEM and not FIML. imputation() / imputed() branch on uncertainty_status, never on is.na(std_error). See docs/design/13-missing-data.md. |

Maturity hardening

Effect objects now expose partial-uncertainty diagnostics. Parametric effect calls warn and attach attr(x, "uncertainty_issues") when a node is non-converged or lacks a usable fixed-effect covariance, so the canonical beta-binomial survival node no longer contributes zero parameter uncertainty silently. Non-finite effect draws are counted in attr(x, "value_issues"), summaries return NA rather than NaN when every draw is unusable, and log-link predictions are clamped below floating-point overflow instead of creating Inf values that are later dropped.
path_effects() now has a live fitted-SEM integration test covering both per-mediator and per-component attribution through a fitted Gaussian mediator with sigma ~ x into a log-link Poisson outcome.
CRAN-prep hygiene now includes a copyright-holder role in Authors@R and redirect-free pkgdown URLs. drmSEM is still not CRAN-submittable while drmTMB and symbolizer are GitHub Remotes; that remains an explicit external release blocker.

symbolizer integration (`Suggests`)

symbolize.drm_sem() walks a fitted drm_sem in topological order, calls symbolizer::symbolize() on each node’s underlying drmTMB fit, and returns a symbolized_drm_sem. The renderer S3 methods as_latex.symbolized_drm_sem(), equations.symbolized_drm_sem(), and assumption_table.symbolized_drm_sem() collate the per-node output with topological-order section headers, so a whole SEM gets publication-quality equations directly from the fitted object — no hand-maintained LaTeX in vignettes or papers.
Covers the model-specification layer only: per-node distributional equations, families, links, assumptions. The SEM layer (DAG, d-separation, Fisher’s C, the simulation-based effect calculus) stays with drmSEM itself, because the distribution-mediated effect has no closed-form symbolic representation by construction.
symbolizer is a Suggests dependency; the integration is gated on requireNamespace("symbolizer"). See vignette("equations-via-symbolizer").

Documentation: rigorous referencing audit

Added structured @references blocks (Rdpack \insertRef{}) to every user-facing help page (?dsep, ?fisher_c, ?direct_effects, ?indirect_effects, ?total_effects, ?path_effects, ?compare, ?drm_sem, ?drm_psem, ?drm_node, ?covary, ?drm_pair, ?drm_cycle, ?drm_composite, ?drm_phylo_cov, ?as_lavaan, ?standardize, …) so every borrowed concept is attributed at the point of use.
All 12 vignettes now carry bibliography: ../inst/REFERENCES.bib and link-citations: true; inline mentions (“Pearl”, “Imai”, “VanderWeele”, “Shipley”, “Lefcheck”, “Felsenstein”, “Pagel”, “Hansen”) have been converted to [@Key] citations. Each vignette ends with a ## References section, rendered by Pandoc.
New inst/REFERENCES.bib and an expanded paper.bib carry the canonical sources for piecewise SEM, d-separation, counterfactual mediation, distributional regression, phylogenetic comparative methods, information criteria, and simulation methodology. DESCRIPTION now imports Rdpack with RdMacros: Rdpack so \insertRef{} resolves at install time.
The three drmSEM constructions built on the cited foundations — the any-component d-separation test, the distribution-mediated effect row, and the BIC-style CBIC information criterion — are explicitly flagged as package-specific contributions (with the underlying components cited) in ?compare, ?dsep, ?indirect_effects, vignette("model-selection"), and vignette("effect-decomposition").
README.md gains a “Methodological background” section listing the canonical entry points and pointing to inst/REFERENCES.bib.

drmSEM 0.5.0

This release closes the cyclic / feedback-graph milestone (roadmap §0.5) and ships the dev-line surface accumulated since 0.2.0. The DAG-only restriction is lifted for declared feedback motifs, with an equilibrium estimand and a pure-R fixed-point propagator; undeclared cycles remain a hard error. The feature highlights below are grouped by area. Consistent feedback estimation (IV/2SLS or a joint likelihood), full sigma-separation, distributional feedback equilibria, and the joint bivariate fit remain engine-dependent and are carried forward to the live-drmTMB lane (see docs/memory/CODEX_HANDOFF.md).

Sampler and propagation fixes

drm_sample_family() and effect propagation now match current drmTMB parameterization for the common sampler families in live recovery tests. Default fitted dpars such as sigma are carried into prediction engines even when no explicit sigma ~ ... formula is declared, and lognormal nodes now use mu = meanlog, sigma = sdlog, with mean mediation propagating exp(mu + sigma^2 / 2) (OQ-1, V-57..V-60).

DAG plot: faithful legend + readable edges

plot.drm_sem() now builds its legend from the components actually drawn (sourced from the same style function as the edges, so the two can never drift) instead of always listing all seven distributional components — the hero/landing-page DAG previously showed nu/hu/sd(.)/rho12 swatches for paths that did not exist. Covariance rows are added only for the classes truly present.
Parallel paths between one pair (e.g. a mu and a sigma arrow on the same edge) are now fanned onto separate arcs instead of overlapping into a single line, and a layout = matrix (optionally row-named) can be supplied for a fixed, crossing-free diagram. A node-fill legend (endogenous response vs exogenous predictor) is drawn, and hu gets a distinct linetype so it no longer relies on colour alone to separate from zi (colour-blind safety).
The legend construction is now a tested pure helper (drm_path_legend()), so a regression that re-introduces phantom legend entries fails CI.
Composite measurement edges: plot.drm_sem(show = "all") now draws each [drm_composite()] construct’s indicators pointing into the construct as steel-blue arrows (indicators shown as distinctly-filled nodes), so a formative measurement model reads apart from the structural paths (OQ-15).

Outcome functionals across the effect API (OQ-11)

All three effect functions now report the effect on a chosen functional of the outcome distribution, not just the mean: target = "mean" / "p_gt" / "p_zero" / "var" / "quantile" (new, with a prob argument). target already rode direct_effects() / total_effects(); it now also rides indirect_effects() (effect = "controlled"), where every leg reports the contrast on the functional and the mean-/distribution-mediated split still closes (indirect = mean_mediated + distribution_mediated). This is where distribution-mediated paths earn their keep — a path into sigma/zi/nu can move a tail probability or quantile while leaving E[Y] nearly unchanged.
The "quantile" target reports the prob-quantile of the simulated outcome (kernel-validated: a path into sigma shifts the upper quantile but not the median).
direct_effects() / total_effects() gain functional = c("simulate", "analytic"). "analytic" evaluates the functional in closed form from the predicted parameters (no Monte-Carlo noise) for the gaussian and poisson families — exact var/p_gt/p_zero/quantile. Other families abort with a pointer back to "simulate" (their sigma↔︎dispersion scale is the OQ-1 open item); analytic needs mean mediation (method = "gcomp").
Fix: the functional engine now honours the mediator-propagation mode ("mean" vs "distribution") instead of always simulating the mediator, so the controlled decomposition is non-degenerate for a non-mean target. effect = "natural" remains mean-only (the cross-world functional contrast is open, OQ-8/OQ-11); a feedback SEM stays mean-only (the equilibrium response).

Validation wave 2 harness + newcomer docs

inst/validation/generate.R + a validation article now provide the wave-2 coverage/calibration harness (effect-CI coverage against a known-effect linear-Gaussian DGP, and model-selection recovery rate), mirroring the OQ-6 calibration pattern (cached .rds + a vignette that renders with or without the cache). Full replicate runs happen in the live lane; see docs/design/12-coverage-calibration.md.
The README quick-start and the intro vignette now open with a concrete biological question and show illustrative (clearly-marked, engine-free) paths() and indirect_effects() output so a newcomer can see the component-labelled and distribution_mediated rows and learn how to read them.

Simulation-based recovery grid (validation wave 1)

A campaign of numerical-recovery tests on real drmTMB fits now exercises the machinery end-to-end, not just for finiteness (V-45..V-73; see the new docs/design/11-validation-matrix.md). The effect decomposition is recovered across the family×link grid (gaussian, poisson, nbinom2, binomial, beta_binomial, beta, Gamma, lognormal) — mean-mediated equals the fitted-coef product / a predict_parameters() do-contrast, total = direct + indirect closes, and distribution_mediated matches the Jensen-gap magnitude from fitted params (the V-7 live-fit follow-up). Each family’s drm_sample_family() mean and variance are checked against drmTMB::simulate() (closing OQ-1 gaps), outcome functionals (p_zero/var/p_gt) are recovered, and the standardization sigma_E pipeline, composite-as-response, feedback equilibrium vs the fitted reduced form, and natural NDE/NIE are validated on live fits. A nonlinear feedback fixed point is added at the kernel tier (V-73). Tweedie / zero_one_beta inflation / student nu samplers stay flagged for the live lane. Wave 2 (effect-CI coverage, d-sep Type-I/power, model-selection recovery rate) is the calibration layer, tracked in CODEX_HANDOFF.md.

Standardization: GLM mean-path `sigma_E` (OQ-4)

The latent standardization of a mu path on a constant-variance link now divides by sqrt(Var(eta) + sigma_E^2), adding the link’s theoretical latent-scale error variance — logit pi^2/3, probit 1, cloglog pi^2/6 (Grace et al. 2018; piecewiseSEM’s latent.linear). This corrects the earlier mild over-standardization of GLM mean paths. Identity-link mu and non-mu components (sigma/zi/sd(*)) are unchanged; the log-link families’ mean-dependent variance term remains deferred. Validated in closed form (V-44).

Interop (graph interchange)

Graph interchange, not a fitting bridge. A new pure-R interop layer (R/interop.R) translates a drmSEM component-labelled graph to and from the neighbouring ecosystems’ text formats. drmSEM still never fits its own likelihoods, and lavaan/brms fitting interop stays out of the 0.x scope.
as_lavaan(sem) (and as_lavaan(dag)) emits a lavaan model-syntax string: one y ~ x1 + x2 regression per endogenous node (the mean structure) and one y1 ~~ y2 line per declared covariance edge (covariances()).
Honesty: lavaan syntax cannot express a distributional-component path (an arrow into sigma, zi, nu, hu, sd(group), rho12). as_lavaan() therefore collapses to the mean structure and reports every dropped non-mu path — both as a dropped attribute and via a one-time cli message. A non-mean path is never silently misrepresented as a lavaan mean regression.
from_lavaan(syntax) parses lavaan syntax back into a drmSEM graph skeleton: ~ regressions become per-response node formulas in a drm_dag(), and ~~ lines become covary() declarations. Reflective measurement (=~) lines are ignored with a warning (reflective measurement needs a joint likelihood, out of 0.x scope). Pure string parsing — nothing is evaluated or fitted, so from_lavaan(as_lavaan(sem)) round-trips the directed mean structure and the covariance edges.
as_dot(sem) (and as_dot(dag)) exports the component-labelled DAG as a Graphviz DOT string: one labelled edge per typed edge, with non-mean paths dashed/greyed. Unlike lavaan, DOT keeps every component path.

Feedback / cyclic motifs (0.5.0, grammar + equilibrium engine)

drm_cycle("y1", "y2") declares a feedback motif; drm_sem() / drm_psem() gain a feedback = argument that accepts it. Cycles remain a hard error unless declared — a declared motif is condensed into one topological layer, so the DAG check still rejects every undeclared cycle. cycles(sem) lists the declared motifs.
Honest fitting. Node-wise ML of a declared cycle is inconsistent under simultaneity, so drm_sem() warns: consistent estimation (IV/2SLS or a joint likelihood) is an engine capability, not something drmSEM fakes.
d-separation drops independence claims among a motif’s nodes (DAG d-separation does not hold across a cycle; full sigma-separation is deferred).
Equilibrium total effects (0.5.x). total_effects() now reports the equilibrium response of a feedback SEM, iterating the mean-propagation map to its fixed point (the mediation column reads "equilibrium"); if the feedback diverges (no stable equilibrium, spectral radius >= 1) the estimate is NA with a warning — never a fabricated number. direct_effects() (the controlled direct effect, which does not traverse the cycle) also works. The mean/distribution decomposition through a cycle is out of scope, so indirect_effects() / path_effects() refuse a feedback SEM and point to total_effects(). The internal propagate_fixedpoint() carries a spectral-radius / max-iter stability guard; closed-form tests confirm it recovers the linear reduced form (I − B)⁻¹ Γ (V-42) and that the equilibrium total effect equals the reduced-form total effect of the exposure (V-43). Design: docs/design/10-cyclic-feedback.md.
New vignette “Feedback cycles: reciprocal causation and the equilibrium effect (drm_cycle)” works a reciprocal activity ⇄ stress pair: why a cycle breaks the DAG machinery (simultaneity bias in fitting; an equilibrium, not a path product, for effects), the opt-in drm_cycle() / feedback = grammar and cycles(), the equilibrium estimand ((I − B)⁻¹ Γ = the walk sum, generalized to the fixed point; stability needs ρ(B) < 1), and the total_effects() equilibrium output (mediation = "equilibrium", target = "mean" only, NA on divergence) with indirect_effects() / path_effects() refusing a feedback SEM. It is prominently honest that node-wise ML of a declared cycle is inconsistent and drmSEM does not fake consistency. Engine-dependent chunks are illustrative-only (eval = has_engine, FALSE).

Effect decomposition: paired Monte-Carlo and honest framing

Bug fix (intervals). indirect_effects() now computes the controlled-direct, mean-mediated, and distribution-mediated legs from one shared coefficient draw per replicate (drm_decomp_legs()), mirroring the natural-effect branch. Previously the three legs were drawn independently, so with the default seed = NULL the mean_mediated / distribution_mediated intervals subtracted unrelated parameter draws and were inflated (not a valid paired contrast). Point estimates are unchanged; the reported intervals are now common-random-numbers (paired) contrasts that isolate the propagation mode.
New asserted tests lock the shipped decomposition path (not only the kernels): the additive identity indirect = mean_mediated + distribution_mediated (V-36), the lognormal Jensen-gap closed form and its sign flip through drm_decomp_legs() (V-37), the linear-outcome zero (V-38), a two-mediator chain (V-39), seed reproducibility (V-40), and an end-to-end live-fit check that the distribution-mediated channel is real and the decomposition closes (V-41).
Honest framing. Documentation now describes distribution_mediated as a Jensen-gap / interventional-mediation term — non-zero only when a higher moment of the mediator responds to the exposure and the outcome is curved in that mediator — identified under the usual mediation assumptions plus a correctly specified mediator distribution. The novelty is positioned as implementational (a mediator’s sigma/zi/shape as a first-class causal target), with the estimand credited to the interventional/distributional mediation literature (Pearl; Imai et al.; VanderWeele 2015; Vansteelandt & Daniel 2017).

Bivariate nodes (0.4, grammar layer)

drm_pair() declares a bivariate (joint two-response) node — two response formulas, two families, an optional rho12 ~ x residual-correlation model (a directed path into the rho12 component), and an auto-detected higher-level corpair edge where the two responses share a grouping level. It is the bivariate counterpart of drm_node().
drm_expand_pair() bridges a pair onto the shipped covariance-edge grammar (two drm_node() sub-nodes + covary() edges), so the residual (rho12) and higher-level (corpair) arcs flow through covariances() / basis_set() / dsep() unchanged — the documented hook point for the 0.4 joint fit.
rho12() and corpairs() accessors report the declared residual and higher-level correlation edges (of a drm_pair or a drm_sem), kept separate from paths(). The estimate is NA: drmSEM never fits its own likelihoods, so reading a fitted correlation back needs a live bivariate drmTMB fit (the 0.4 engine deliverable). The declaration is never a fabricated estimate. See docs/design/07-bivariate-covariance-edges.md.
plot(sem) now draws covariance edges as double-headed arcs — solid grey for a residual rho12, dashed grey for a higher-level corpair — so the three edge classes (directed path, residual covariance, higher-level covariance) are visually distinct. A new show = c("all", "paths") argument toggles the arcs off ("paths" draws the directed structural edges only).
New vignette “Bivariate nodes: two responses, one correlation (drm_pair)” walks an animal-personality example (activity & boldness) through drm_pair(), print(), rho12(), corpairs(), and drm_expand_pair(), motivating the difference between a directed activity -> boldness path and the residual rho12 / higher-level corpair covariance arcs. It is emphatic that the correlation estimates are NA by construction — declaring a bivariate node is not fitting one, and the joint bivariate drmTMB fit is the 0.4 engine step. See docs/design/07-bivariate-covariance-edges.md.

Path-specific effects (OQ-5)

path_effects(effect = "natural") reports each mediator’s cross-world natural indirect effect, with an identified column flagged FALSE under a recanting witness (another mediator that is a descendant of from and an ancestor of this one), after Avin, Shpitser & Pearl (2005). The detection is pure graph logic; the cross-world computation reuses the validated drm_natural_target kernel.

Latent constructs (0.3)

drm_composite() now reports reliability (Cronbach’s alpha, an internal-consistency measure for reflective indicator sets — unclamped, NA for a single indicator) and accepts standardize = TRUE (mean-0 / sd-1 score). A new summary.drm_composite() shows the loadings table, first-PC variance, and reliability. A composite may now be used as a response (a node modelling the construct), not only a predictor.
New vignette “Latent constructs: composites and reliability” walks through formative vs reflective-flavoured (PCA) constructs, construct quality, and construct-as-response; reflective measurement-model SEM stays out of scope (needs a joint likelihood → 0.4), documented in docs/design/09-latent-variables.md.

Calibration (OQ-6)

The live-drmTMB lane generated inst/calibration/calibration-results.rds (drmTMB 0.1.3.9000, Git SHA 17b1321, 14,400 calibration replicates). All five pre-specified acceptance checks pass, so V-17 is promoted to validated for the OQ-6 mean-only / distributional / cross-link calibration grid only; broader Fisher’s C / d-sep settings remain claim-scoped until separately tested.

drmSEM 0.2.0

Second release. Post-0.1 work: a unified effect-API surface, first-class covariance-edge and composite-construct grammars, per-mediator and per-component path attribution, finalized standardization conventions, and a scaffolded Fisher’s C calibration study. The any-component d-separation calibration remains experimental until its (compute-heavy) study is run; everything else is CI-validated against a live drmTMB or kernel-validated by closed-form tests.

Documentation

New vignette “Covariance edges, composites, and path attribution” walks through covary()/covariances(), drm_composite()/loadings(), and path_effects(). Every exported function now carries an example.

Path-specific effects (OQ-5, per-mediator)

path_effects(object, from, to, through=) decomposes the indirect effect into a per-mediator contribution: an inclusion effect (T({Mj}) - direct, only Mj responds) and an exclusion effect (T(all) - T(all \\ Mj), Mj’s marginal given the rest), plus total_indirect and an explicit interaction_remainder. The pieces sum to the total only in the additive case; the remainder is reported, never forced to zero. Model-based attribution, not a nonparametric path-specific identification claim. The cross-world natural variant is the OQ-5 follow-up.
path_effects(by = "component") splits each mediator’s effect into a mean_channel and one channel per non-mean component (sigma_channel, zi_channel, … – the drop when that component is frozen at its reference value), plus a component_remainder for the part that does not separate cleanly under a nonlinear outcome.

Composite latent constructs (0.3, first increment)

drm_composite(name, indicators, weights, method = c("fixed", "pca")) declares a composite (formative) construct — a weighted sum or first-PC index of observed indicator columns. It is materialized as an ordinary column before fitting, so a node formula can use it as a predictor or response with no engine change.
drm_sem() / drm_psem() gain a composites = argument; loadings(sem) reports the indicator-to-construct loadings, kept separate from paths().
Reflective (measurement-model) latent variables remain out of scope (they need a joint likelihood); see docs/design/09-latent-variables.md.

Inference hardening (towards 0.2)

Analytic effect cross-checks are now asserted tests (test-analytic-effects.R, pure-R, no engine): the Gaussian identity-link mean-mediated effect equals the coefficient product a*b*w; a non-mean (sigma) path contributes exactly nothing to the mean channel and the distribution-mediated effect goes to zero when the outcome is linear in the mediator; the distribution-mediated effect matches the lognormal closed form (and flips sign) across a downstream nonlinearity; natural and controlled effects diverge under an exposure-mediator interaction by the predicted amount; and the Poisson Pr(Y>0) and Gaussian Var(Y) outcome-functional effects match their closed forms.
Standardization conventions finalized and documented (OQ-4; see docs/design/08-standardization.md and ?standardize): standardized coefficients are reported on the link scale only; factor predictors keep SD = 1 (raw per-contrast effect, lavaan std.nox convention); the latent divisor is per-component, so sigma/zi paths standardize on their own link scale. A Gelman 2-SD opt-in for continuous-vs-factor comparability and a theoretical- variance term for GLM mean paths are noted as tracked refinements. ## Covariance edges: rho12 and corpair (OQ-14, grammar layer)
covary(y1, y2, level = ) declares a covariance edge between two responses — a residual correlation (rho12, within-observation, level = NULL) or a higher-level random-effect correlation (corpair, between-unit, a grouping level). Covariance edges are double-headed arcs: they carry no direction and no mediated effect.
drm_sem() / drm_psem() gain a covariances = argument that takes covary() declarations and validates them against the node records.
covariances(sem) reports residual and higher-level edges separately, kept out of paths() (which stays directed-only, including any x -> rho12 path).
basis_set() / dsep() are now covariance-aware: a declared rho12 / corpair edge between y1 and y2 drops the y1 _||_ y2 | predictors independence claim (Shipley’s bidirected-edge rule).
This is the pure-R grammar/d-separation layer. drm_pair() (joint bivariate fitting), rho12() / corpairs() accessors that read a live fit, and double-headed-arc plotting need a bivariate drmTMB fit and remain on the roadmap (OQ-14).

Unified effect-API surface (OQ-12)

direct_effects(), total_effects(), and indirect_effects() now share one argument vocabulary: uncertainty = c("parametric", "none", "bootstrap"), nsim (inner distributional realizations), and population = c("conditional", "marginal"). total_effects() selects mediation with method = c("gcomp", "simulate").
The previous mediation, draw, and n_sim arguments are deprecated aliases — they still work but emit a deprecation warning, and the new argument wins when both are supplied. B (the number of uncertainty replicates) is unchanged. No simulation kernel changed.
Not-yet-implemented choices fail fast with a pointer to the tracking question: uncertainty = "bootstrap" (refit bootstrap, OQ-10) and population = "marginal" (marginalizing over the random-effect distribution, OQ-9).
direct_effects() gains target / threshold, so a controlled direct effect can be read on an outcome functional (p_gt, p_zero, var) as well as the mean, matching total_effects().

drmSEM 0.1.0

First public release of drmSEM — a distributional piecewise structural equation modelling layer built on the drmTMB fitting engine. drmSEM does not fit its own likelihoods; each endogenous node is one drmTMB fit and the system is piecewise over a DAG. Causal paths can target any distributional component (mean, scale, zero-inflation, shape, random-effect scale, residual correlation), and effects are estimated by Monte-Carlo counterfactual propagation rather than coefficient products.

Building a SEM

drm_node() specifies one endogenous node (a drmTMB formula + family).
drm_sem() assembles a SEM declaratively, fitting each node, and drm_psem() assembles one from already-fitted drmTMB objects. Both return the same object.

Component-labelled paths

paths() returns one row per fitted coefficient, labelled by the distributional component it targets (mu, sigma, nu, zi, hu, sd(*), rho12). A path to sigma or zi is never reported as a mean effect.
plot() draws the DAG with edges styled by component.

d-separation (any modelled component)

basis_set(), dsep(), and fisher_c() test missing arrows. A missing arrow asserts X has no effect on any modelled component of Y, via a likelihood-ratio test of Y’s node augmented with X on every component sub-model. Fisher’s C combines the claim p-values (C = -2*sum(log p), 2k df).

Simulation-based effects

direct_effects(), indirect_effects(), and total_effects() propagate a do-style intervention through the fitted DAG by Monte-Carlo, on the response scale. Coefficient products are never used across non-Gaussian or cross-link paths.
The indirect decomposition reports mean_mediated and distribution_mediated parts, so effects flowing through a mediator’s scale, zero-inflation, or shape are visible rather than collapsed into a mean effect.
indirect_effects(effect = "natural") adds the cross-world natural decomposition — natural_direct, natural_indirect, and mediated_interaction (Pearl; Imai, Keele & Yamamoto) — alongside the default controlled-direct / simulation-indirect split, validated on the linear-Gaussian recovery case.
total_effects(target = c("p_gt", "p_zero", "var")) reports distribution-mediated effects on outcome functionals beyond the mean — Pr(Y > threshold), Pr(Y = 0), and Var(Y) — with the p_zero effect recovery-tested against the Poisson zero-probability change.
standardize() reports standardized path coefficients on the component’s link scale.

Model comparison (confirmatory model sets)

drm_dag() captures one unfitted candidate causal model (a set of node formulas), and drm_model_set() collects named candidates into a comparison set — drmSEM’s analogues of phylopath::define_model_set().
A new model-selection article explains the CBIC-default workflow, shows how to request the CICc support ranking, and summarizes the cached C-3 recovery results.
compare() now defaults to criterion = "CBIC", the BIC-style ranking that passed the model-selection recovery check. The comparison table still reports CICc and CBIC deltas and weights side by side; use criterion = "CICc" for the phylopath-style support ranking.
best() returns the top-ranked fitted SEM under the selected criterion; average() returns criterion-weighted (conditional) model-averaged standardized path coefficients.

Phylogenetic covariance

drm_phylo_cov() builds a phylogenetic relatedness matrix from an ape tree under a fixed evolutionary model ("BM", "lambda", "OU", "kappa"), ready to feed a node via relmat(1 | species, K = K). The evolutionary parameter is fixed by the caller (a grid), not jointly estimated.
dsep() augment-refits phylogenetic nodes correctly, evaluating each refit in the SEM’s captured fitting environment so structured-effect objects (the tree / relatedness matrix) resolve.

Plots

plot.drm_effect() draws an effect (forest) plot for an effect decomposition.

Diagnostics

check_sem() reports unsupported surfaces, convergence, and sampler coverage.

Status

Early and experimental; APIs may change before a stable release. Kernel logic is validated by tests that run without the engine; the full drmTMB-integration path is validated in CI where drmTMB is compiled. See the validation ledger in docs/memory/VALIDATION_LEDGER.md.