Feedback cycles: reciprocal causation and equilibrium effects (drm

In standard piecewise structural equation models (Shipley 2000, 2009; Lefcheck 2016), causal dependencies flow in one direction along a directed acyclic graph (DAG). However, many ecological, behavioural, and physiological systems feature reciprocal feedback loops—where two variables simultaneously cause each other:

Physiology & Behaviour: An animal’s glucocorticoid stress levels suppress its movement, while high activity elevates physiological stress (activity ⇄ stress).
Ecosystems: Predator and prey densities mutually regulate one another (predator ⇄ prey).
Plant–Soil Interactions: Plant biomass enriches microbial activity, which in turn enhances nutrient uptake and plant growth (plant ⇄ soil).

Such feedback breaks the topological ordering assumed by DAG algorithms. This vignette explains how drmSEM enables researchers to:

Declare feedback motifs explicitly using drm_cycle().
Compute equilibrium total effects via Banach fixed-point iteration across location ( $\mu$ ), scale ( $\sigma$ ), shape ( $\nu$ ), and zero-inflation ( $\text{zi}$ ) components.
Diagnose system stability using the spectral radius $\rho(B) < 1$ and empirical Lipschitz contraction constants.
Understand the honest boundary between simulation-based equilibrium calculation and simultaneity bias in node-wise estimation.

The reciprocal challenge: DAGs versus cycles

A reciprocal pair introduces two distinct statistical challenges:

Simultaneity bias during fitting: In $y_1 \rightleftarrows y_2$ , $y_1$ is a predictor in the $y_2$ equation, but $y_1$ is also a descendant of $y_2$ . Therefore, $y_1$ correlates with $y_2$ ’s error term. Ordinary node-wise maximum likelihood is inconsistent (Bollen 1989, 1996).
Equilibrium propagation: In a DAG, the total effect of an upstream driver $x$ is computed by propagating changes down a finite sequence of directed paths (topological sweep). In a cyclic system, there is no topological order: an intervention reverberates around the loop indefinitely until the system relaxes to an equilibrium steady state.

1. Declaring a feedback motif: `drm_cycle()`

By default, drmSEM enforces strict acyclicity: any undetected cycle results in a hard error. To permit feedback, you must explicitly declare the motif using drm_cycle():

# Declare a reciprocal pair between activity and stress
motif <- drm_cycle("activity", "stress")
motif
#> <feedback motif> activity <-> stress

When handed to drm_sem(..., feedback = motif) or drm_psem(..., feedback = motif): - The topological sort condenses the motif into a super-node (relaxed topological sort). - Independence claims between motif members ($y_1 \perp\mkern-10mu\perp y_2$) are suppressed from the d-separation basis set. - All undeclared cycles remain strict errors.

# Simulate data with reciprocal feedback
n <- 300
x <- rnorm(n)
# Reduced-form equilibrium simulation:
# activity = 0.4*stress + 0.5*x + e1
# stress   = -0.3*activity + 0.4*x + e2
# B = matrix(c(0, -0.3, 0.4, 0), 2, 2)
# (I - B)^(-1) is well-conditioned because spectral radius rho(B) = sqrt(0.12) = 0.346 < 1
I_minus_B_inv <- solve(diag(2) - matrix(c(0, -0.3, 0.4, 0), 2, 2))
Gamma <- matrix(c(0.5, 0.4), 2, 1)
total_mult <- I_minus_B_inv %*% Gamma

e1 <- rnorm(n, 0, 0.3)
e2 <- rnorm(n, 0, 0.3)
reduced_e <- t(I_minus_B_inv %*% rbind(e1, e2))

dat <- data.frame(
  x = x,
  activity = total_mult[1, 1] * x + reduced_e[, 1],
  stress   = total_mult[2, 1] * x + reduced_e[, 2]
)

# Fit the SEM with declared feedback:
sem_cyclic <- drm_sem(
  activity = drm_node(drmTMB::bf(activity ~ stress + x), family = stats::gaussian()),
  stress   = drm_node(drmTMB::bf(stress ~ activity + x), family = stats::gaussian()),
  data = dat,
  feedback = drm_cycle("activity", "stress")
)
#> ℹ Fitting node "activity"
#> ✔ Fitting node "stress" [116ms]
#> 
#> ℹ Fitting node "stress"
#> Warning: A declared feedback motif is fitted node-wise; ordinary ML is inconsistent
#> under simultaneity.
#> ℹ Consistent estimation (IV/2SLS or a joint likelihood) is an engine
#>   capability; equilibrium effects use the fixed-point propagator on the
#>   supplied coefficients (see docs/design/10-cyclic-feedback.md).
#> ✔ Fitting node "stress" [97ms]
#> 

# Inspect declared cycles
cycles(sem_cyclic)
#> <drmSEM feedback motifs: 1>
#> motif 1: activity <-> stress

2. Multi-component Banach fixed-point iteration

In linear Gaussian models, the steady-state equilibrium under an intervention $\Delta x$ has a closed-form matrix solution (Bollen 1989):

$y = (I - B)^{-1} \Gamma x$

where $B$ is the direct within-loop coefficient matrix (with zeros on the diagonal) and $\Gamma$ contains direct effects of exogenous variables.

Because $(I - B)^{-1} = I + B + B^2 + B^3 + \dots$ , the matrix inverse is the infinite sum over all walks around the loop.

Distributional and non-linear systems

In distributional regression, responses may follow non-Gaussian families (e.g., Gamma, Negative Binomial, Beta, Student- $t$ ) or have location-scale links ( $\log \sigma$ , $\text{logit } \text{zi}$ , $\log \nu$ ). For non-linear systems, no closed-form matrix inverse exists.

drmSEM resolves non-linear feedback using vectorized Banach fixed-point iteration (Banach 1922). Let $\mathbf{s}^{(t)}$ denote the joint state vector containing the predicted values across all active distributional components:

$\mathbf{s}^{(t)} = \begin{pmatrix} \mu_{\text{activity}}^{(t)} \\ \sigma_{\text{activity}}^{(t)} \\ \mu_{\text{stress}}^{(t)} \\ \sigma_{\text{stress}}^{(t)} \end{pmatrix}$

At each iteration $t = 1, 2, \dots, \text{max\_iter}$ : 1. Each endogenous node re-predicts all its distributional components from the current working values of its parents. 2. The state vector updates with adaptive relaxation: $\mathbf{s}^{(t+1)} = (1 - \alpha) \mathbf{s}^{(t)} + \alpha \, \mathcal{T}(\mathbf{s}^{(t)})$ 3. The iteration terminates when the Euclidean norm of change is below tolerance: $\|\mathbf{s}^{(t+1)} - \mathbf{s}^{(t)}\|_\infty < \text{tol} \quad (\text{default } 10^{-8})$

3. Stability diagnostics and Lipschitz contraction

A feedback system reaches a stable equilibrium if and only if the iteration is a strict contraction (Banach 1922; Forré and Mooij 2017). drmSEM monitors stability through two complementary diagnostics:

1. Spectral radius $\rho(B)$

For the linear Gaussian sub-system, drmSEM computes the eigenvalues of $B$ :

$\rho(B) = \max_i |\lambda_i(B)|$

$\rho(B) < 1$ : The feedback is damped. The series converges to a unique, stable equilibrium.
$\rho(B) \ge 1$ : The feedback is explosive / unstable. Trips around the loop amplify without bound, and no equilibrium effect is defined.

2. Empirical Lipschitz contraction constant

For non-linear and distributional loops, drmSEM calculates the empirical ratio of successive state deltas:

$L^{(t)} = \frac{\|\mathbf{s}^{(t+1)} - \mathbf{s}^{(t)}\|}{\|\mathbf{s}^{(t)} - \mathbf{s}^{(t-1)}\|}$

If $L < 1$ , the map is contractive on the support of the data.

4. Equilibrium total effects

To estimate the downstream consequence of an intervention in a feedback system, call total_effects():

eff_total <- total_effects(sem_cyclic, from = "x", to = "stress")
eff_total
#> <drmSEM effect>
#>  from     to    scale   mediation target estimate conf.low conf.high
#>     x stress response equilibrium   mean   0.2395   0.2104    0.2664

Key properties of equilibrium effects

mediation = "equilibrium": The output explicitly flags that the estimate was computed by fixed-point equilibrium simulation rather than a finite acyclic path product.
Honest divergence handling: If a system is unstable ( $\rho(B) \ge 1$ or divergence), drmSEM returns NA with an informative warning rather than a fabricated or infinite number.
Decomposition constraints: Because path-specific mediation decompositions require a topological sort, indirect_effects() and path_effects() are disabled on cyclic SEMs, directing users to total_effects(). direct_effects() remains valid (measuring the controlled direct effect holding mediators constant).

5. d-Separation and basis-set scoping

In acyclic graphs, conditional independence claims are evaluated using d-separation (Pearl 2009). In cyclic graphs, the correct independence criterion is $\sigma$ -separation (Forré and Mooij 2017).

Until full $\sigma$ -separation is implemented, basis_set() automatically drops conditional independence claims among members of declared feedback motifs, scoping Fisher’s C test strictly to the acyclic components of the network.

Summary of best practices

Explicitly declare cycles: Never let a feedback loop exist undeclared; use drm_cycle("y1", "y2").
Check stability: Ensure $\rho(B) < 1$ and verify that fixed-point iteration reaches convergence.
Report equilibrium effects: Use total_effects() to capture the fully reverberated steady-state response.
Be honest about simultaneity: If fitting node-wise regressions, acknowledge that standard ML coefficients in reciprocal loops are subject to simultaneity bias.

References

Banach, Stefan. 1922. “Sur Les Opérations Dans Les Ensembles Abstraits Et Leur Application Aux équations Intégrales.” Fundamenta Mathematicae 3 (1): 133–81.

Bollen, Kenneth A. 1989. Structural Equations with Latent Variables. Wiley.

Bollen, Kenneth A. 1996. “An Alternative Two Stage Least Squares (2SLS) Estimator for Latent Variable Equations.” Psychometrika 61 (1): 109–21. https://doi.org/10.1007/BF02296961.

Forré, Patrick, and Joris M. Mooij. 2017. Markov Properties for Graphical Models with Cycles and Latent Variables. arXiv:1710.08775. https://doi.org/10.48550/arXiv.1710.08775.

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. 2009. Causality: Models, Reasoning, and Inference. 2nd ed. Cambridge University Press.

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.

Feedback cycles: reciprocal causation and equilibrium effects (drm_cycle)