Effect Types: dma vs crumble

# Optional: enable parallel cross-fitting
library(future)
plan(multisession)

Overview

dma supports four mediation effect decompositions from the general framework of Liu, Williams, Rudolph & Díaz (2024):

Effect Estimands Requires moc
Natural (N) NDE, NIE, ATE No
Organic (O) ODE, OIE No
Randomized Interventional (RI) RDE, RIE Yes
Recanting Twins (RT) 4-pathway decomposition Yes

This vignette demonstrates each effect type on the weight_behavior dataset from the mma package and compares results with crumble, which uses the same estimation framework but with mlr3superlearner (ensemble of mean, GLM, MARS, and random forest) instead of engression for outcome regressions. Crumble results are taken from its README using crossfit_folds = 1L, epochs = 20L and learners = c("mean", "glm", "earth", "ranger").

Setup

library(dma)
library(mma)

data(weight_behavior)
weight_behavior <- na.omit(weight_behavior)

set.seed(2345)

The shift functions define the modified treatment policies — here, static interventions setting everyone to sports = 1 (control) or sports = 2 (treated):

d0 <- \(data, trt) factor(rep(1, nrow(data)), levels = c("1", "2"))
d1 <- \(data, trt) factor(rep(2, nrow(data)), levels = c("1", "2"))

Natural Effects

Natural direct and indirect effects (Pearl, 2001). The NDE captures the effect of shifting treatment while holding the mediator at its natural value under the control regime; the NIE captures the remaining indirect pathway.

No mediator-outcome confounders (moc) are needed — the cross-world counterfactual is handled by the nested pseudo-outcome cascade.

dma(
  data = weight_behavior,
  trt = "sports",
  outcome = "bmi",
  covar = c("age", "sex", "tvhours"),
  mediators = c("exercises", "overweigh"),
  d0 = d0, d1 = d1,
  effect = "N",
  control = dma_control(crossfit_folds = 1L, num_epochs = 100L, riesz_epochs = 20L)
)
#> ✔ Fitting outcome regressions... 1/1 folds [40.8s]
#> ✔ Computing alpha n density ratios... 1/1 folds [40.5s]
#>
#> ══ Results `dma()` ═════════════════════════════════════════
#>
#> ── Direct Effect
#>       Estimate: -0.043
#>     Std. error: 0.18
#> 95% Conf. int.: -0.396, 0.31
#>
#> ── Indirect Effect
#>       Estimate: 0.835
#>     Std. error: 0.232
#> 95% Conf. int.: 0.379, 1.29
#>
#> ── Average Treatment Effect
#>       Estimate: 0.792
#>     Std. error: 0.305
#> 95% Conf. int.: 0.195, 1.389

Comparison with crumble:

Estimand dma crumble Difference
ATE 0.792 (0.305) 1.029 (0.280) -0.237
NDE -0.043 (0.180) 0.017 (0.172) -0.060
NIE 0.835 (0.232) 1.012 (0.218) -0.177

Both packages agree on the qualitative story: a small, insignificant direct effect and a large, significant indirect effect. Point estimates differ modestly — expected since engression (distributional regression) and mlr3superlearner (ensemble of classical learners) are fundamentally different learning strategies, and both runs use only 1 cross-fitting fold.

Organic Effects

Organic direct and indirect effects (Lok, 2015). Unlike natural effects, organic effects use the observed mediator distribution rather than a counterfactual one, avoiding cross-world assumptions.

dma(
  data = weight_behavior,
  trt = "sports",
  outcome = "bmi",
  covar = c("age", "sex", "tvhours"),
  mediators = c("exercises", "overweigh"),
  d0 = d0, d1 = d1,
  effect = "O",
  control = dma_control(crossfit_folds = 1L, num_epochs = 100L, riesz_epochs = 20L)
)
#> ✔ Fitting outcome regressions... 1/1 folds [41s]
#> ✔ Computing alpha n density ratios... 1/1 folds [40.1s]
#>
#> ══ Results `dma()` ═════════════════════════════════════════
#>
#> ── Organic Direct Effect
#>       Estimate: -0.085
#>     Std. error: 0.177
#> 95% Conf. int.: -0.431, 0.262
#>
#> ── Organic Indirect Effect
#>       Estimate: 0.898
#>     Std. error: 0.233
#> 95% Conf. int.: 0.442, 1.354

Comparison with crumble:

Estimand dma crumble Difference
ODE -0.085 (0.177) 0.011 (0.177) -0.096
OIE 0.898 (0.233) 1.028 (0.223) -0.130

Again, qualitative agreement — insignificant direct effect, significant indirect effect — with comparable standard errors.

Randomized Interventional Effects

Randomized interventional direct and indirect effects (Vansteelandt & Daniel, 2017). These effects are identified even in the presence of mediator-outcome confounders (moc) affected by treatment, at the cost of a weaker decomposition that uses a randomized mediator draw from the population.

The moc argument specifies intermediate confounders — here, snack. Estimation requires Z’ permutation via Rsymphony.

dma(
  data = weight_behavior,
  trt = "sports",
  outcome = "bmi",
  covar = c("age", "sex", "tvhours"),
  mediators = c("exercises", "overweigh"),
  moc = "snack",
  d0 = d0, d1 = d1,
  effect = "RI",
  control = dma_control(crossfit_folds = 1L, num_epochs = 100L, riesz_epochs = 20L)
)
#> ✔ Permuting Z-prime variables... 1/1 tasks [2.1s]
#> ✔ Fitting outcome regressions... 1/1 folds [1m 0.7s]
#> ✔ Computing alpha r density ratios... 1/1 folds [1m 16.3s]
#>
#> ══ Results `dma()` ═════════════════════════════════════════
#>
#> ── Randomized Direct Effect
#>       Estimate: 0.014
#>     Std. error: 0.182
#> 95% Conf. int.: -0.343, 0.371
#>
#> ── Randomized Indirect Effect
#>       Estimate: 0.951
#>     Std. error: 0.245
#> 95% Conf. int.: 0.47, 1.431

Comparison with crumble:

Estimand dma crumble Difference
RDE 0.014 (0.182) 0.016 (0.177) -0.002
RIE 0.951 (0.245) 1.030 (0.230) -0.079

The closest agreement of any effect type — direct effects are nearly identical, and indirect effects are within one standard error.

Recanting Twins

Recanting twins pathway decomposition (Vo et al., 2024). This is the finest decomposition, separating the total effect into four causal pathways plus an intermediate confounding term:

dma(
  data = weight_behavior,
  trt = "sports",
  outcome = "bmi",
  covar = c("age", "sex", "tvhours"),
  mediators = c("exercises", "overweigh"),
  moc = "snack",
  d0 = d0, d1 = d1,
  effect = "RT",
  control = dma_control(crossfit_folds = 1L, num_epochs = 100L, riesz_epochs = 20L)
)
#> ✔ Permuting Z-prime variables... 1/1 tasks [2.1s]
#> ✔ Fitting outcome regressions... 1/1 folds [1m 46.8s]
#> ✔ Computing alpha n density ratios... 1/1 folds [1m 17.1s]
#> ✔ Computing alpha r density ratios... 1/1 folds [1m 17.5s]
#>
#> ══ Results `dma()` ═════════════════════════════════════════
#>
#> ── Path: A -> Y
#>       Estimate: 0.241
#>     Std. error: 0.191
#> 95% Conf. int.: -0.133, 0.615
#>
#> ── Path: A -> Z -> Y
#>       Estimate: 0.074
#>     Std. error: 0.02
#> 95% Conf. int.: 0.035, 0.113
#>
#> ── Path: A -> Z -> M -> Y
#>       Estimate: -0.019
#>     Std. error: 0.019
#> 95% Conf. int.: -0.055, 0.018
#>
#> ── Path: A -> M -> Y
#>       Estimate: 0.866
#>     Std. error: 0.232
#> 95% Conf. int.: 0.412, 1.32
#>
#> ── Intermediate Confounding
#>       Estimate: -0.208
#>     Std. error: 0.052
#> 95% Conf. int.: -0.311, -0.106

Comparison with crumble:

Pathway dma crumble Difference
A → Y 0.241 (0.191) 0.037 (0.184) 0.204
A → Z → Y 0.074 (0.020) -0.020 (0.024) 0.094
A → Z → M → Y -0.019 (0.019) -0.001 (0.010) -0.018
A → M → Y 0.866 (0.232) 1.051 (0.216) -0.185
Intermediate conf. -0.208 (0.052) -0.013 (0.026) -0.195

The dominant pathway (A → M → Y) is consistent across both packages. The finer decomposition shows more variation between learners, which is expected — the RT decomposition involves the most nuisance parameters (both natural and randomized cascades), so differences in the outcome regression learner compound through more levels. Both packages agree that the A → M → Y path is by far the largest and most significant component.

Summary

Point estimates on weight_behavior data (crossfit_folds = 1)
Effect Estimand dma crumble
Natural ATE 0.792 1.029
NDE -0.043 0.017
NIE 0.835 1.012
Organic ODE -0.085 0.011
OIE 0.898 1.028
Rand. Interventional RDE 0.014 0.016
RIE 0.951 1.030
Recanting Twins A->Y 0.241 0.037
A->Z->Y 0.074 -0.020
A->Z->M->Y -0.019 -0.001
A->M->Y 0.866 1.051
Int. Conf. -0.208 -0.013
ATE 0.954 1.054

Both packages target the same estimands using the same semiparametric framework, but differ in how nuisance parameters are estimated:

Qualitative conclusions are identical across all four effect types: a large, significant indirect effect operating through the mediators (exercises, overweigh), with a small and generally insignificant direct effect. Point estimate differences reflect the different learners and the inherent variability from using a single cross-fitting fold.

Distributional View: Oracle Potential Outcomes

The weight_behavior results above are point estimates from real data. To illustrate what each effect type captures distributionally, we simulate from a DGP with a mediator-outcome confounder (Z) where all potential outcomes are known.

DGP: \(W \sim N(0,1)\), \(A \in \{0,1\}\) equiprobable, \(Z = 0.3A + 0.4W + \varepsilon_Z\) (mediator-outcome confounder), \(M = 0.5A + 0.3W + 0.2Z + \varepsilon_M\), \(Y = 0.8A + M^2 + 0.5Z + 0.4W + \varepsilon_Y\).

The non-linear \(M^2\) creates distributional effects that go beyond location shifts: the indirect pathway changes the shape of the outcome distribution.

set.seed(42)
n_oracle <- 10000
W_o <- rnorm(n_oracle)
eps_Z <- rnorm(n_oracle, sd = 0.3)
eps_M <- rnorm(n_oracle, sd = 0.4)
eps_Y <- rnorm(n_oracle, sd = 0.5)

# Helper: potential outcomes given (a_Y, a_M, a_Z) for the treatment entering each equation
po <- function(a_Y, a_M, a_Z) {
  Z <- 0.3 * a_Z + 0.4 * W_o + eps_Z
  M <- 0.5 * a_M + 0.3 * W_o + 0.2 * Z + eps_M
  Y <- 0.8 * a_Y + M^2 + 0.5 * Z + 0.4 * W_o + eps_Y
  list(Y = Y, M = M, Z = Z)
}

# Standard four counterfactuals (treatment enters all equations the same way)
Y00 <- po(0, 0, 0)$Y  # Y(0, M(0), Z(0))
Y11 <- po(1, 1, 1)$Y  # Y(1, M(1), Z(1))

# Cross-world: Y(1, M(0)) — treat directly but mediator at control
# Z also at control level (natural effects operate A→Y directly)
Y10 <- po(1, 0, 0)$Y  # Y(1, M(0), Z(0))
Y01 <- po(0, 1, 1)$Y  # Y(0, M(1), Z(1))

# Plot helper
plot_po <- function(Y00, Y11, Y10, Y01, title, subtitle = NULL) {
  df <- rbind(
    data.frame(Y = Y00, PO = "Y(0, M(0))"),
    data.frame(Y = Y11, PO = "Y(1, M(1))"),
    data.frame(Y = Y10, PO = "Y(1, M(0))"),
    data.frame(Y = Y01, PO = "Y(0, M(1))")
  )
  df$PO <- factor(df$PO, levels = c("Y(0, M(0))", "Y(1, M(0))", "Y(0, M(1))", "Y(1, M(1))"))

  ggplot(df, aes(x = Y, color = PO)) +
    geom_density(linewidth = 0.9) +
    scale_color_manual(values = c(
      "Y(0, M(0))" = "#377eb8", "Y(1, M(1))" = "#e41a1c",
      "Y(1, M(0))" = "#984ea3", "Y(0, M(1))" = "#ff7f00"
    )) +
    labs(title = title, subtitle = subtitle,
         x = "Outcome (Y)", y = "Density", color = "Potential Outcome") +
    theme_minimal() +
    theme(legend.position = "bottom")
}

decomp_table <- function(Y00, Y11, Y10, Y01) {
  data.frame(
    Effect = c("ATE", "NDE", "NIE"),
    Value = round(c(mean(Y11) - mean(Y00), mean(Y10) - mean(Y00), mean(Y11) - mean(Y10)), 3)
  )
}

Natural Effects (N)

Natural effects use cross-world counterfactuals: NDE sets \(A=1\) but keeps the mediator at its \(A=0\) level; NIE shifts the mediator while holding the direct treatment.

plot_po(Y00, Y11, Y10, Y01,
        title = "Natural Effects: Potential Outcome Distributions",
        subtitle = "Y = 0.8A + M² + 0.5Z + 0.4W + eps — M² creates right-skew shift")

decomp_table(Y00, Y11, Y10, Y01)
#>   Effect Value
#> 1    ATE 1.262
#> 2    NDE 0.800
#> 3    NIE 0.462

The gap between Y(0,M(0)) (blue) and Y(1,M(0)) (purple) is the NDE — a location shift from the direct \(0.8A\) term. The gap between Y(1,M(0)) (purple) and Y(1,M(1)) (red) is the NIE — a shape-changing shift because \(M^2\) introduces right-skew when \(M\) is larger under treatment.

Organic Effects (O)

Organic effects use the observed mediator distribution rather than a counterfactual one. The ODE compares treatment regimes holding \(M\) at the observed (mixture) distribution; the OIE captures the remaining pathway.

For visualization, the organic decomposition uses the same four potential outcomes, but the ODE/OIE partition differs conceptually — the ODE marginalizes over the observed \(M\) distribution rather than conditioning on \(M(0)\).

# Organic effects: compare treated vs control under observed M distribution
# ODE: E[Y(1,M)] - E[Y(0,M)] where M is from the observed distribution
# OIE: Total - ODE
# Simulate observed M (proper mixture: draw A ~ Bernoulli(0.5) per unit)
A_obs <- rbinom(n_oracle, 1, 0.5)
Z_obs <- 0.3 * A_obs + 0.4 * W_o + eps_Z
M_obs <- 0.5 * A_obs + 0.3 * W_o + 0.2 * Z_obs + eps_M

Y_1_Mobs <- 0.8 * 1 + M_obs^2 + 0.5 * Z_obs + 0.4 * W_o + eps_Y
Y_0_Mobs <- 0.8 * 0 + M_obs^2 + 0.5 * Z_obs + 0.4 * W_o + eps_Y

df_o <- rbind(
  data.frame(Y = Y00, PO = "Y(0, M(0)): Control"),
  data.frame(Y = Y11, PO = "Y(1, M(1)): Treated"),
  data.frame(Y = Y_1_Mobs, PO = "Y(1, M_obs): ODE reference"),
  data.frame(Y = Y_0_Mobs, PO = "Y(0, M_obs): ODE baseline")
)
df_o$PO <- factor(df_o$PO, levels = c("Y(0, M(0)): Control", "Y(0, M_obs): ODE baseline",
                                        "Y(1, M_obs): ODE reference", "Y(1, M(1)): Treated"))

ggplot(df_o, aes(x = Y, color = PO)) +
  geom_density(linewidth = 0.9) +
  scale_color_manual(values = c(
    "Y(0, M(0)): Control" = "#377eb8", "Y(1, M(1)): Treated" = "#e41a1c",
    "Y(1, M_obs): ODE reference" = "#984ea3", "Y(0, M_obs): ODE baseline" = "#ff7f00"
  )) +
  labs(title = "Organic Effects: Potential Outcome Distributions",
       subtitle = "ODE: gap between Y(1,M_obs) and Y(0,M_obs); OIE: remainder",
       x = "Outcome (Y)", y = "Density", color = "Potential Outcome") +
  theme_minimal() +
  theme(legend.position = "bottom")


data.frame(
  Effect = c("ODE", "OIE"),
  Value = round(c(mean(Y_1_Mobs) - mean(Y_0_Mobs),
                  (mean(Y11) - mean(Y00)) - (mean(Y_1_Mobs) - mean(Y_0_Mobs))), 3)
)
#>   Effect Value
#> 1    ODE 0.800
#> 2    OIE 0.462

Randomized Interventional Effects (RI)

RI effects handle mediator-outcome confounders (Z) by drawing the mediator from a randomized population distribution rather than an individual’s counterfactual. This avoids the cross-world assumption at the cost of a weaker decomposition.

# RI: draw M from population distribution under each treatment
# RDE: E[Y(1, M_pop(0))] - E[Y(0, M_pop(0))]
# RIE: E[Y(1, M_pop(1))] - E[Y(1, M_pop(0))]
# M_pop(a) is drawn independently from the marginal distribution of M(a)

# Shuffle mediator values to break individual-level pairing
set.seed(99)
M0_pop <- sample(po(0, 0, 0)$M)  # random draw from M(0) distribution
M1_pop <- sample(po(1, 1, 1)$M)  # random draw from M(1) distribution
Z0_pop <- sample(po(0, 0, 0)$Z)  # also shuffle Z to break M-Z confounding
Z1_pop <- sample(po(1, 1, 1)$Z)

Y_1_M0pop <- 0.8 * 1 + M0_pop^2 + 0.5 * Z0_pop + 0.4 * W_o + eps_Y
Y_0_M0pop <- 0.8 * 0 + M0_pop^2 + 0.5 * Z0_pop + 0.4 * W_o + eps_Y
Y_1_M1pop <- 0.8 * 1 + M1_pop^2 + 0.5 * Z1_pop + 0.4 * W_o + eps_Y

df_ri <- rbind(
  data.frame(Y = Y_0_M0pop, PO = "Y(0, M_pop(0))"),
  data.frame(Y = Y_1_M0pop, PO = "Y(1, M_pop(0))"),
  data.frame(Y = Y_1_M1pop, PO = "Y(1, M_pop(1))"),
  data.frame(Y = Y11, PO = "Y(1, M(1)): Individual")
)
df_ri$PO <- factor(df_ri$PO, levels = c("Y(0, M_pop(0))", "Y(1, M_pop(0))",
                                          "Y(1, M_pop(1))", "Y(1, M(1)): Individual"))

ggplot(df_ri, aes(x = Y, color = PO)) +
  geom_density(linewidth = 0.9) +
  scale_color_manual(values = c(
    "Y(0, M_pop(0))" = "#377eb8", "Y(1, M_pop(0))" = "#984ea3",
    "Y(1, M_pop(1))" = "#e41a1c", "Y(1, M(1)): Individual" = "#999999"
  )) +
  labs(title = "Randomized Interventional Effects: Potential Outcome Distributions",
       subtitle = "RDE: purple vs blue; RIE: red vs purple (gray = individual natural for reference)",
       x = "Outcome (Y)", y = "Density", color = "Potential Outcome") +
  theme_minimal() +
  theme(legend.position = "bottom")


data.frame(
  Effect = c("RDE", "RIE", "ATE"),
  Value = round(c(mean(Y_1_M0pop) - mean(Y_0_M0pop),
                  mean(Y_1_M1pop) - mean(Y_1_M0pop),
                  mean(Y_1_M1pop) - mean(Y_0_M0pop)), 3)
)
#>   Effect Value
#> 1    RDE 0.800
#> 2    RIE 0.462
#> 3    ATE 1.262

The randomized mediator (drawn from the population distribution) washes out individual-level confounding between M and Z, which is why the RI distributions (solid colors) differ slightly from the natural ones. The gray line shows the individual-level Y(1,M(1)) for comparison.

Recanting Twins (RT)

The RT decomposition separates the total effect into four causal pathways through all combinations of the mediator and mediator-outcome confounder channels:

# RT four pathways via (a_Y, a_Z, a_M) assignment:
# A→Y direct: vary A entering Y only
# A→Z→Y: vary A entering Z (and thus Y through Z), hold M at A=0
# A→M→Y: vary A entering M only, hold Z at A=0
# A→Z→M→Y: vary A entering both Z and M

Y_base <- po(0, 0, 0)$Y  # baseline: all at control

# A→Y: set a_Y=1, keep a_M=0, a_Z=0
Y_AY <- po(1, 0, 0)$Y
path_AY <- mean(Y_AY) - mean(Y_base)

# A→Z→Y: set a_Z=1 (Z shifts, which feeds into M and Y), keep a_Y=0, a_M=0
# But M also gets Z(1) input — we need to be careful about what "a_M=0" means here
# For RT, we construct: a_Y=0, a_Z=1, a_M=0
Z1_o <- 0.3 * 1 + 0.4 * W_o + eps_Z
M_Z1 <- 0.5 * 0 + 0.3 * W_o + 0.2 * Z1_o + eps_M  # M with Z(1) but a_M=0
Y_AZY <- 0.8 * 0 + M_Z1^2 + 0.5 * Z1_o + 0.4 * W_o + eps_Y
path_AZY <- mean(Y_AZY) - mean(Y_base)

# A→M→Y: set a_M=1, keep a_Y=0, a_Z=0
Z0_o <- 0.3 * 0 + 0.4 * W_o + eps_Z
M_A1 <- 0.5 * 1 + 0.3 * W_o + 0.2 * Z0_o + eps_M  # M with a_M=1 but Z(0)
Y_AMY <- 0.8 * 0 + M_A1^2 + 0.5 * Z0_o + 0.4 * W_o + eps_Y
path_AMY <- mean(Y_AMY) - mean(Y_base)

# A→Z→M→Y: set a_Z=1, a_M=1, keep a_Y=0
M_A1Z1 <- 0.5 * 1 + 0.3 * W_o + 0.2 * Z1_o + eps_M
Y_AZMY <- 0.8 * 0 + M_A1Z1^2 + 0.5 * Z1_o + 0.4 * W_o + eps_Y
path_AZMY <- mean(Y_AZMY) - mean(Y_base)

# Plot key pathway distributions
df_rt <- rbind(
  data.frame(Y = Y_base, PO = "Baseline: Y(0,M(0),Z(0))"),
  data.frame(Y = Y_AY, PO = "A->Y: Y(1,M(0),Z(0))"),
  data.frame(Y = Y_AMY, PO = "A->M->Y: Y(0,M(1),Z(0))"),
  data.frame(Y = Y11, PO = "Full treatment: Y(1,M(1),Z(1))")
)
df_rt$PO <- factor(df_rt$PO, levels = c("Baseline: Y(0,M(0),Z(0))",
                                          "A->Y: Y(1,M(0),Z(0))",
                                          "A->M->Y: Y(0,M(1),Z(0))",
                                          "Full treatment: Y(1,M(1),Z(1))"))

ggplot(df_rt, aes(x = Y, color = PO)) +
  geom_density(linewidth = 0.9) +
  scale_color_manual(values = c(
    "Baseline: Y(0,M(0),Z(0))" = "#377eb8",
    "A->Y: Y(1,M(0),Z(0))" = "#984ea3",
    "A->M->Y: Y(0,M(1),Z(0))" = "#ff7f00",
    "Full treatment: Y(1,M(1),Z(1))" = "#e41a1c"
  )) +
  labs(title = "Recanting Twins: Key Pathway Distributions",
       subtitle = "Four pathways decompose the gap between blue (baseline) and red (full treatment)",
       x = "Outcome (Y)", y = "Density", color = "Regime") +
  theme_minimal() +
  theme(legend.position = "bottom")


data.frame(
  Pathway = c("A->Y (direct)", "A->Z->Y (through Z only)",
              "A->M->Y (through M only)", "A->Z->M->Y (through Z then M)",
              "Sum of pathways", "ATE (total)"),
  Value = round(c(path_AY, path_AZY, path_AMY, path_AZMY,
                  path_AY + path_AZY + path_AMY + path_AZMY,
                  mean(Y11) - mean(Y_base)), 3)
)
#>                         Pathway Value
#> 1                 A->Y (direct) 0.800
#> 2      A->Z->Y (through Z only) 0.153
#> 3      A->M->Y (through M only) 0.249
#> 4 A->Z->M->Y (through Z then M) 0.462
#> 5               Sum of pathways 1.665
#> 6                   ATE (total) 1.262

The direct pathway (A→Y, purple vs blue) creates a location shift of \(\approx 0.8\). The mediator pathway (A→M→Y, orange) shifts the distribution and changes its shape because of the \(M^2\) non-linearity — this is precisely what distributional mediation reveals that mean-only decompositions miss.

Note that the RT pathway decomposition does not generally sum exactly to the ATE — the remainder is attributed to the intermediate confounding term, which captures the interaction between the Z and M channels.

References