PHASE_18_SNMS.md

Phase 18 — G-estimation of Structural Nested Mean Models (SNMMs)

Status: COMPLETE — all chunks shipped (18a–18f, 18h–18k); 18g dropped (gesttools archived/broken, delicatessen used instead)

Depends on: Phase 2 (point infra), Phase 4 (treatment-model machinery), Phase 5 (longitudinal data shape), Phase 6 (effect-modification parser)

Motivation

SNMs are the third major approach to causal inference with time-varying confounding, alongside g-formula (gcomp, ICE) and IPW-MSM. Robins' original 1986 framework put them on equal footing with the other two; they have been less popular mainly because their outputs (per-stage blip parameters) are harder to explain than marginal counterfactual means, and because software support has been thin (gesttools is recent; geex is general-purpose but low-level).

causatr already covers the other two pillars of the Robins triangle (g-formula via gcomp/ICE; IPW-MSM via the self-contained density-ratio engine), so SNMs are the natural third leg. The reason to add them now — the reason I am reopening the "out of scope" call I made earlier — is that they are the only approach in causatr's target scope (GLM/GAM-based, sandwich-variance-based, no ML) that correctly handles effect modification by time-varying covariates. MSMs cannot (§ "Why SNMs, why now"); the g-formula can in principle but requires a correctly-specified high-dimensional outcome model; SNMs parameterize exactly the object of interest — the per-stage blip γ(a, l̄ₖ, āₖ₋₁; ψ) — and identify it via a moment condition that uses the treatment model as an "instrument."

Why SNMs, why now

Three applied scenarios causatr cannot answer today:

1. Time-varying effect modification. "Does the effect of a treatment increment at visit k depend on the patient's covariate trajectory up through k?" MSMs cannot condition on time-varying L in their target parameter (conditioning on post-treatment variables; Hernán & Robins Ch. 12.6, Robins 2000). G-formula can in principle but requires a correct high-dimensional outcome model. SNMs parameterize the blip directly and identify it cleanly under treatment-model correctness — this is the native use case documented in Vansteelandt & Joffe 2014 Section 5.
2. Near-positivity settings. When propensities approach 0 or 1, IPW weights explode and the MSM becomes unstable; the g-formula has to extrapolate the outcome model beyond covariate support. SNMs' estimating equations degrade more gracefully because the blip parameters are identified through covariation structure rather than via inverse weighting. This is the classical "structural nested models are the method of choice in pharmacoepi" argument (Robins & Hernán 2009).
3. Optimal dynamic treatment regimes as a side-product. Once the per-stage blip parameters are estimated, the optimal regime $d^*(l̄_k) = \arg\max_a γ(a, l̄_k; \hat\psi)$ falls out by inspection. IPW-MSM needs a whole extra estimation stage (value-search, Q-learning, etc.) for the same object. While DTR estimation is not a headline Phase 18 deliverable, the blip parameters support it downstream.

Scenario 1 is the headline. It also retroactively fixes a Phase 6 correctness gap (PHASE_6_INTERACTIONS.md did not say "modifier must be baseline" for IPW, which is now flagged by an explicit check and a doc note; see § "Relationship to Phase 6").

Scope

1. Point-treatment SNMM with linear blip γ(a, l; ψ) = a · (ψ₀ + Σⱼ ψⱼ · mⱼ), where mⱼ are user-supplied effect modifiers (baseline OR time-varying — both allowed, unlike IPW-MSM).
2. Longitudinal SNMM with per-stage blip γₖ(aₖ, l̄ₖ, āₖ₋₁; ψ); the stage-k blip may depend on time-varying covariates from l̄ₖ.
3. G-estimation via moment equations using the treatment model as an instrument.
4. Continuous and binary outcomes with an additive-linear blip (most studied and most robust case). Multiplicative / log-linear blips are deferred.
5. Sandwich variance via stacked EE: K propensity blocks + blip parameter block.
6. Bootstrap variance (straightforward — blip estimation wrapped in boot::boot()).
7. Truth-based simulation tests against analytical blip parameters on linear-Gaussian DGPs.
8. External cross-check against gesttools (Wallace et al.) on shared DGPs.

Non-scope

• Survival SNMs (SNFTMs, SNCFTMs). Accelerated-failure-time and cumulative-failure-time structural nested models are a distinct literature with their own parameterizations and estimators (Robins 1992; Joffe et al. 2012). Out of scope; these belong in the separate survival package (etverse/survatr).
• Multiplicative / log-linear blips for binary outcomes with causal OR / RR. Requires a different g-estimating equation and awkward finite-sample behavior. Deferred.
• Optimal dynamic regimes as a first-class output. The blip parameters already encode the optimal regime, but constructing a causatr_regime object with inference is a separate methodological question (value-function inference, DTR confidence sets). Out of scope.
• Rank-preserving SNMs with heterogeneous treatment effects at the individual level. Phase 18 assumes "no effect modification by unobserved factors beyond what the blip captures" (Robins 1994; Vansteelandt & Joffe 2014) — the semiparametric version. The stronger rank-preservation assumption is not needed for the additive-linear mean case.
• Multivariate treatment SNMs. Out of scope; compose with Phase 8 if and when it ships.
• Matching + SNMs. Architecturally incompatible (SNMs are a moment-based estimator, not a reweighting or matching one).

Design

Point-treatment SNMM

For a single-time-point treatment with linear blip $\gamma(a, l; \psi) = a \cdot (\psi_0 + \sum_j \psi_j \cdot m_j)$, the g-estimating equation is a residualised-treatment moment condition: $$ E_n\Big[ \big(A_i - \hat{E}[A \mid L_i]\big) \cdot \mathbf{m}_i \cdot \big(Y_i - \gamma(A_i, L_i; \psi)\big) \Big] = 0, $$ where $\mathbf{m}i = (1, m{1,i}, \dots)^\top$ is the modifier vector. Under correct specification of $\hat{E}[A \mid L]$ (the treatment model), the residual $A_i - \hat{E}[A \mid L_i]$ is orthogonal to $L_i$ by construction, so the moment condition identifies $\psi$ without requiring the outcome-model component to be correctly specified. This is exactly the 2SLS-with-residualized-treatment structure, specialised to a causal-effect-modifier parameterisation.

Solution. For linear blip this has a closed form: $$ \hat\psi = \left(\sum_i \mathbf{m}_i^{\otimes 2} \cdot (A_i - \hat{E}[A \mid L_i]) \cdot A_i\right)^{-1} \sum_i \mathbf{m}_i \cdot (A_i - \hat{E}[A \mid L_i]) \cdot Y_i, $$ a single matrix inversion. For nonlinear blips (not in Phase 18 scope) the moment equation becomes nonlinear in $\psi$ and rootSolve::multiroot() or equivalent is required.

What's "new" vs IPW: only the estimating equation. The treatment model $\hat{E}[A \mid L]$ is fit by the same fit_treatment_model() machinery Phase 4 already uses, so for SNMs under binary / continuous / categorical / count treatment, Phase 4's propensity-family dispatch carries over unchanged.

Longitudinal SNMM

Per-stage blip $\gamma_k(a_k, \bar{l}k, \bar{a}{k-1}; \psi)$. Robins' g-estimating equation at stage $k$ uses the transformed outcome $$ H(\psi)(k) = Y - \sum_{j \geq k} \gamma_j(A_j, \bar{L}j, \bar{A}{j-1}; \psi), $$ which is the hypothetical outcome an individual would have had if treatment from $k$ onward had been set to 0. Under sequential exchangeability and correct specification of the treatment models $f(A_k \mid \bar{A}_{k-1}, \bar{L}k)$, $H(\psi^*)(k) \perp A_k \mid \bar{A}{k-1}, \bar{L}k$ at the truth $\psi^*$ (Robins 1994). The g-estimating equation that exploits this is $$ \sum_k E_n\Big[ \big(A_k - \hat{E}[A_k \mid \bar{A}{k-1}, \bar{L}k]\big) \cdot q_k(\bar{A}{k-1}, \bar{L}_k) \cdot H(\psi)(k) \Big] = 0, $$ summed over all periods $k \in {0, \dots, K}$, where $q_k$ is a user-chosen (for efficiency) vector of functions of history — typically $q_k = \mathbf{m}_k$, the modifier vector at stage $k$, which may include time-varying components of $\bar{L}_k$. For linear blips the system remains linear in $\psi$ and solves by matrix inversion.

Relationship to Phase 6 (effect modification)

The Phase 6 parser parse_effect_mod() (R/effect_modification.R) already detects A:modifier terms in confounders = ~ L + sex + A:sex. Phase 18 reuses the same parser to build the blip parameterisation: every A:modifier (including time-varying modifiers from the ICE lag-expansion machinery, lag1_A:L_k etc.) becomes a blip component with its own $\psi$ parameter. The difference from Phase 6 under IPW is that Phase 18's SNM identifies the blip under treatment-model correctness alone, without the baseline-only modifier restriction that IPW-MSM requires. This closes the Phase 6 gap documented in PHASE_6_INTERACTIONS.md: users who want time-varying effect modification under IPW-style identification should switch to estimator = "snm".

Stacked sandwich variance

The estimating-equation system is $$ \omega(L, A, Y; \alpha, \psi) = \begin{pmatrix} s_\alpha^{\text{treatment}}(L, A) \ \omega_\psi^{\text{g-est}}(L, A, Y; \alpha, \psi) \end{pmatrix}, $$ with $s_\alpha$ stacked across $K$ time points for the longitudinal case. The bread is block-triangular (the treatment score does not depend on $\psi$); the $\omega_\psi$ row has cross-derivatives $\partial \omega_\psi / \partial \alpha$ because the residual $A_k - \hat{E}[A_k \mid \cdot]$ depends on $\hat\alpha$. The sandwich-variance primitives handle this with the same pattern used for IPW: prepare_model_if() on each treatment model, numDeriv::jacobian on the g-estimating-equation closure for the cross-derivative, vcov_from_if() to aggregate the per-individual IF on $\hat\psi$.

Clustering. For repeated measures within individual (typical longitudinal data), the per-individual IF is aggregated across periods within individual before the crossprod, matching the cluster-robust convention established in variance_if_matching() (aggregate score within cluster, then square). This is the natural cluster structure for longitudinal SNMs.

Contrasts

contrast() on an SNM fit returns the blip parameters $\hat\psi$ directly, with standard errors and CIs. For the user-facing "effect" interpretation, two default summaries:

• Average blip effect in the study population: $(1/n) \sum_i \gamma(a, L_i; \hat\psi) - \gamma(0, L_i; \hat\psi)$ averaged over a user-supplied treatment value $a$ (or $a = 1$ for binary). This is an ATE-like scalar.
• Per-modifier-stratum blip effect: for each level of a categorical modifier, the blip value. This is the time-varying effect modification result.

The contrast() output object is causatr_result as usual; the fit_type = "snm" attribute dispatches to an SNM-aware print.causatr_result() and plot.causatr_result() that foreground the blip-parameter table and a modifier-stratum effect plot.

Treatment types

SNMs in Phase 18 support the same treatment types as Phase 4 IPW, with the same family dispatch:

Treatment type	Supported?	Notes
binary	✓	$\hat{E}[A \mid L]$ via logit GLM.
continuous	✓	$\hat{E}[A \mid L]$ via gaussian GLM; this is the canonical SNMM case.
categorical ($k > 2$)	✓	Residualisation via multinomial nnet::multinom; blip becomes level-specific $\gamma(a_k, l; \psi_k)$. Adds complexity; may be chunked separately.
count	✓	Poisson or NB treatment model via propensity_family.
multivariate	⛔	Out of scope; Phase 8 composition deferred.

Composition with other phases

• Phase 10 (longitudinal IPW) + SNM triangulation. Two parallel approaches to the same longitudinal problem. Under correct specification both are consistent; disagreement is a diagnostic red flag. An explicit triangulation test (SNM blip × IPW-MSM marginal mean on a simulated DGP) is the scientific payoff of having both in the same package.
• Phase 14 (IPCW) × SNM. Censoring weights can be incorporated into the SNM g-estimating equation (Yiu & Su 2022; Boatman & Vock 2019): weight the stage-k moment condition by cumulative IPCW up to $k$. Subsection to be added to PHASE_14_IPCW.md once Phase 18 ships.
• Phase 9 (inference infrastructure). Survey weights and cluster-robust SE at the individual level compose straightforwardly with the SNM stacked sandwich. The longitudinal-cluster aggregation (§ "Clustering") is a prerequisite; full survey-weighted SNMs wait for Phase 9.
• Survival. Survival SNMs (SNFTMs / SNCFTMs) are owned by the separate survival package (etverse/survatr), not by Phase 18.

Chunks

Chunk	Scope	Status	Depends on
18a	Route estimator = "snm" to fit_snm(); validate linear-blip specification; reject non-linear blips with informative error	✅	Phase 2
18b	Point-treatment SNMM: compute_snm_blip_point() solves the linear g-estimating equation; compute_snm_contrast() returns blip param table or averaged blip effect; variance_if_snm() stacked EE sandwich; validated against delicatessen + DTRreg	✅	18a, Phase 4
18b½	Treatment-free outcome model for efficiency augmentation. causat(..., treatment_free = ~ L) enables joint estimation of (β, ψ) following Vansteelandt & Joffe (2014) / DTRreg's tf.mod. The treatment-free model absorbs L→Y variance, reducing SEs by 30–45%. Stacked sandwich: (α_trt, θ_joint = (β, ψ)). Validated against delicatessen (3 DGPs) and DTRreg (EM case: exact match on point + SE).	✅	18b
18c	Phase-6 parser integration: parse_effect_mod() already produces the modifier list; wire its output into the blip parameterisation; both baseline and time-varying modifiers accepted (point case)	✅	18b, Phase 6
18d	Longitudinal SNMM: fit_snm() fits $K$ per-period treatment models; compute_snm_blip_longitudinal() uses backward sequential g-estimation (Robins 1994) with per-stage blip parameters $\psi_k$; cluster-aggregated sandwich via stacked EE with cross-stage derivatives and per-individual IF aggregation; validated against DTRreg	✅	18b, Phase 5
18e	Time-varying EM truth-based test: 2-period DGP with $M_0 = 1{L_0 > 0}$ (baseline) and $M_1 = 1{L_1 > 0}$ (post-treatment); both stages have modifier in blip; truth $(\psi_{00}, \psi_{0M}, \psi_{10}, \psi_{1M}) = (1.15, 2, 2, 2)$; DTRreg cross-check with history=0; IPW-MSM bias demonstration	✅	18d
18f	Triangulation test: SNM longitudinal blip-averaged effect vs Phase 10 IPW-MSM + ICE g-comp on a shared DGP; three pairwise checks + collider-bias negative control; delicatessen cross-check via shared fixture data; history=0 (no-lag Markov model) support for all longitudinal estimators	✅	18d, Phase 10
18g	External cross-check against gesttools::gestMultiple() — DROPPED: gesttools archived from CRAN (2025-04-26), broken NAMESPACE on GitHub, DataCombine dependency also archived. delicatessen cross-checks in 18b/18f serve as the primary external oracle instead	⛔ dropped	18d
18h	Categorical / count treatment extensions (residualisation via multinomial / Poisson / NB treatment models)	✅	18b
18i	Bootstrap variance: blip estimation wrapped in boot::boot(); per-individual cluster aggregation for longitudinal	✅	18b, 18d
18j	S3 dispatch: print.causatr_result() / plot.causatr_result() / summary / tidy / coef on estimator = "snm" — foreground blip-parameter table + per-modifier-stratum pointrange plot	✅	18b
18k	Documentation, vignette (snm.qmd — categorical/count/bootstrap/S3 sections), FEATURE_COVERAGE_MATRIX.md update, CLAUDE.md update	✅	18a–18j

Supported combination matrix

Full accounting of which dimensions interact with SNM and their status. This matrix governs what Phase 18 ships vs defers vs rejects.

Treatment type

Treatment	SNM support	Chunk	Notes
Binary	✅ shipped	18b	Logistic propensity via fit_treatment_model()
Continuous	✅ shipped	18b	Gaussian propensity; canonical SNMM case
Categorical ($k > 2$)	✅ shipped	18h	Multinomial residualisation via nnet::multinom; per-level blip parameters
Count (Poisson/NB)	✅ shipped	18h	Via propensity_family dispatch; scalar residual R = A − λ̂
Multivariate	⛔ rejected	—	Out of scope; Phase 8 composition deferred

Treatment timing

Timing	SNM support	Chunk	Notes
Point	✅ shipped	18b	Closed-form g-estimation
Longitudinal	✅ shipped	18d	Per-stage blip via backward sequential g-estimation; cluster-robust sandwich

Outcome family

All outcome families work with the additive-linear blip (the blip enters as $Y - \gamma$ in the moment condition; no link-function interaction). The blip interpretation changes: for non-gaussian outcomes with additive blip, $\psi$ is an additive mean difference on the response scale, not a log-odds or log-rate ratio. Multiplicative/log-linear blips for binary/count outcomes are deferred (see Non-scope).

Family	SNM support	Chunk	Notes
Gaussian	✅ truth-tested	18b	Design doc DGP; delicatessen cross-check
Binomial	✅ truth-tested	18b-ext	Additive blip = risk difference; linear probability DGP at n=3000
Poisson	✅ truth-tested	18b-ext	Additive blip on count scale; linear-in-A DGP
Gamma	✅ truth-tested	18b-ext	Additive blip on mean scale
Negative binomial	✅ truth-tested	18b-ext	Via MASS::negative.binomial(theta) family
Quasibinomial	✅ truth-tested	18b-ext	Same DGP as binomial; overdispersion parameter irrelevant for blip EE
Betareg	✅ truth-tested	18b-ext	Via betareg::betareg; additive blip on (0,1) scale

Estimand

Estimand	SNM support	Notes
ATE	✅ shipped	Default: blip parameter table or averaged blip effect
ATT/ATC	⛔ rejected	SNMs estimate the blip (conditional on covariates), not marginal means under subpopulations. ATT/ATC do not have natural blip analogs under additive linear parameterisation
by-stratified	✅ shipped	Per-stratum averaged blip: global ψ̂ with stratum-specific colMeans(M_s); delta-method SE from global vcov

Contrast type

Contrast	SNM support	Notes
Difference	✅ shipped	treatment_values = c(0, 1) gives $(a_1 - a_0) \cdot (\hat\psi_0 + \sum \hat\psi_j \bar{m}_j)$
Ratio	⛔ rejected	Additive-linear blip parameterises a difference, not a ratio. Multiplicative blips deferred
OR	⛔ rejected	Same reason as ratio

Variance method

Method	SNM support	Chunk	Notes
Sandwich (analytic)	✅ shipped	18b	Stacked EE, block-triangular bread, cross-derivative via numDeriv::jacobian
Sandwich (numeric fallback)	✅ shipped	18b-ext	Defensive tryCatch in variance_if_snm(): catches models without model.matrix/coef/residuals support; re-throws as causatr_snm_no_estfun with pointer to bootstrap
Bootstrap	✅ shipped	18i	Blip re-estimation in boot::boot(); point row-resampling + longitudinal ID-cluster resampling
Cluster-robust	✅ shipped	18b-ext	Explicit cluster = threading via vcov_from_if(); fit-time propagation; dimension-safe subsetting to fit rows

Model class (treatment model)

Model	SNM support	Notes
GLM (stats::glm)	✅ shipped	Default via fit_treatment_model()
GAM (mgcv::gam)	✅ smoke-tested	Via propensity_model_fn = mgcv::gam; bread uses $Vp; as.numeric() wrap for predict.gam() 1-D array bug

Transportability / generalizability

Feature	SNM support	Notes
target = (sampling model)	❌ deferred	Requires fitting treatment model on study rows only, averaging blip over target covariate distribution. Conceptually straightforward but needs design work. See Phase 17 cross-reference

Censoring / IPCW

Feature	SNM support	Notes
censoring = row filter	✅ inherited	Outcome NAs excluded via get_fit_rows() (same as gcomp)
ipcw = TRUE	❌ deferred	Weight moment condition by cumulative IPCW (Yiu & Su 2022). See Phase 14 cross-reference

External / survey weights

Feature	SNM support	Notes
weights = external	✅ shipped	Propagated to both treatment model and blip EE (RM_w = RM * w); sandwich bread and score also weighted; matches manual WLS formula sum(wRY)/sum(wRA)
Survey design	❌ deferred	Full survey-weighted SNMs wait for Phase 9 infrastructure

Effect modification

Feature	SNM support	Chunk	Notes
Baseline modifier	✅ shipped	18b	Via Phase 6 parse_effect_mod()
Time-varying modifier (point)	✅ shipped	18c	Parser wiring; headline feature. Truth-based tests + delicatessen + DTRreg cross-checks; IPW bias demonstration
Time-varying modifier (long.)	✅ shipped	18d–18e	Per-stage blip with $\bar{L}_k$ terms (18d); truth test with $M_0, M_1$ (18e); DTRreg cross-check; IPW bias demonstration
Multiple modifiers	✅ shipped	18b	Validated with 2-modifier DGP + delicatessen

Stabilised weights

Not applicable — SNMs are moment-based, not reweighting estimators.

Missing data

Pattern	SNM support	Notes
Outcome NA (MCAR)	✅ handled	get_fit_rows() excludes NAs; same as gcomp
Treatment NA	🟡 inherited	Via treatment model's na.action; not SNM-specific
Covariate NA	🟡 inherited	Via prepare_data() upstream; not SNM-specific
MAR (multiple imputation)	❌ deferred	Phase 21 causat_mice() composes via Rubin pooling

Intervention types

SNMs do not use interventions = — the blip parameter is the estimand directly. The treatment_values = c(a0, a1) path computes averaged blip effects. This is a structural difference from gcomp/IPW/AIPW, not a missing feature.

Intervention	SNM support	Notes
static()	⛔ rejected	contrast(fit, interventions = ...) → causatr_snm_no_interventions
shift()	⛔ rejected	Same
treatment_values =	✅ shipped	Averaged blip effect with delta-method SE

Diagnostics

Feature	SNM support	Notes
diagnose() dispatch	❌ deferred	SNM diagnostics differ from gcomp/IPW: treatment-residual orthogonality, moment-condition diagnostics, per-modifier blip plots. Phase 11 did not add SNM branch; future chunk

Invariants

• estimator = "snm" requires at least one A:modifier term in confounders for the blip to be non-trivial; otherwise the SNMM reduces to a simple two-stage residualisation estimator of the constant ATE, which is fine but warrants an rlang::inform() note ("no effect modifiers specified — blip reduces to a single ATE parameter").
• The treatment model under SNM must use the same machinery as IPW (fit_treatment_model()), so whatever Phase 4 accepts on the propensity side, Phase 18 accepts on the treatment-residualisation side. If IPW rejects a treatment/intervention combination, so does SNM — but the failure modes differ: IPW rejects static-on-continuous because of Dirac issues; SNM does not have the equivalent issue because it never evaluates counterfactual densities, only residual moments.
• interventions = is not used for SNM estimation — the blip parameter is the estimand directly; there is no need to specify a counterfactual treatment to recover $\hat\psi$. contrast() on an SNM fit therefore takes a different argument signature: contrast(fit, treatment_values = c(0, 1)) returns the average blip effect at the supplied treatment values. Users who try to pass interventions = get a targeted error pointing to treatment_values =.
• Under correct model specification, the SNM point estimate on a DGP must agree with the Phase 10 longitudinal-IPW point estimate up to Monte Carlo noise (chunk 18f). Disagreement indicates a bug or misspecification, not a method difference.
• The Phase 6 doc gap closure (PHASE_6_INTERACTIONS.md "baseline-only modifier under IPW" note) cross-links to this phase as the recommended alternative for users whose modifier is time-varying. Keep the cross-link intact when PHASE_6 or PHASE_18 is edited.

DGP for truth-based tests

Point-treatment with effect modification (chunk 18b)

L ~ N(0, 1),   M = L > 0
A | L ~ N(0.5 · L, 1)                 (continuous treatment)
Y | A, L, M = 2 + 3·A + 1.5·L + 2·A·M + ε,  ε ~ N(0, 1)

Linear blip: γ(a, l, m; ψ) = a · (ψ_0 + ψ_M · m)
True ψ_0 = 3, ψ_M = 2

Longitudinal with time-varying effect modification (chunk 18e — the headline test)

L_0 ~ N(0, 1)
A_0 | L_0 ~ Bernoulli(expit(0.3 · L_0))
L_1 = 0.5 · L_0 + 0.3 · A_0 + ε_L,  ε_L ~ N(0, 0.5)
M_1 = L_1 > 0                               (time-varying modifier — post-treatment!)
A_1 | L_1, A_0 ~ Bernoulli(expit(0.3 · L_1 + 0.2 · A_0))
Y | A_0, A_1, L_0, L_1, M_1 = 2 + 1·A_0 + 2·A_1 + 2·A_1·M_1 + 1.5·L_0 + 0.5·L_1 + ε_Y

Stage-0 blip: γ_0(a_0, l_0; ψ) = a_0 · ψ_{00} = 1·a_0     (no time-varying EM at k=0)
Stage-1 blip: γ_1(a_1, l̄_1, a_0; ψ) = a_1 · (ψ_{10} + ψ_{1M} · M_1) = a_1 · (2 + 2·M_1)

Truth: ψ_{00} = 1, ψ_{10} = 2, ψ_{1M} = 2

The Phase 18 longitudinal SNM must recover (ψ_{00}, ψ_{10}, ψ_{1M}) = (1, 2, 2). The Phase 6 IPW-MSM on the same DGP with confounders = ~ L_0 + L_1 + A_1:M_1 (which is the correct specification if M_1 were baseline) is biased because M_1 is post-treatment — the bias should be visible in the same test.

References

• Robins JM (1986). A new approach to causal inference in mortality studies with a sustained exposure period. Math Modelling 7:1393–1512. (foundational)
• Robins JM (1994). Correcting for non-compliance in randomized trials using structural nested mean models. Comm Stat Theory Methods 23:2379–2412.
• Robins JM (2000). Marginal structural models versus structural nested models as tools for causal inference. In Statistical Models in Epidemiology: The Environment and Clinical Trials, Halloran & Berry (eds.), Springer. (direct comparison, relevant to the time-varying-EM argument)
• Vansteelandt S, Joffe MM (2014). Structural nested models and g-estimation: The partially realized promise. Stat Sci 29:707–731. (modern review — primary teaching reference)
• Naimi AI, Moodie EEM, Auger N, Kaufman JS (2017). Constructing inverse probability weights for continuous exposures: a comparison of methods. Epidemiology 28:709–717. (applied-epi primer for the residualisation argument)
• Moodie EEM, Chakraborty B, Kramer MS (2012). Q-learning for estimating optimal dynamic treatment rules from observational data. Can J Statist 40:629–645. (DTR link)
• Wallace MP, Moodie EEM, Stephens DA (2017). Dynamic treatment regimen estimation via regression-based techniques: Introducing R package DTRreg. J Stat Software 80:1–20.
• Yiu A, Su L (2022). Joint estimation of treatment and covariate-censoring models with applications to structural nested models. Statistica Sinica. (IPCW + SNM composition)
• Boatman JA, Vock DM (2019). Estimating and testing a mean for correlated failure time data with informative censoring. Biostatistics. (Phase 14 × Phase 18 reference)