Methods and literature

This page maps CausalTargeted APIs to the papers and books that define the estimands, identification conditions, and estimators. The implementations are Julia-native analogues of ideas developed in the LMTP / mediation literature (including the R lmtp and crumble packages); they are not line-for-line ports. Full bibliographic entries (with DOIs) are in References. Keys such as diaz2023lmtp match references.bib in the CDCS book. Engine naming (:lmtp / :mediation, not “crumble”) is summarised in NAMING.md.

Boundary with confounded offline reinforcement learning

Recent causal offline-RL work studies continuous-action policy learning when the behaviour policy uses latent state that is only partially recorded in the offline dataset. That problem is adjacent to, but not the same as, LMTP estimation: CausalTargeted estimates declared treatment-policy functionals, whereas offline-RL algorithms learn value functions and policies under sequential feedback. The package therefore does not expose flow-matching Q-learning. A future bridge may consume CausalDynamics observation and policy declarations, but robust lower-bound objectives and expressive continuous-action learners belong in an optional RL extension or separate package.

Modified treatment policies and LMTP

Scientific problem. Deterministic interventions that set a continuous exposure to a fixed value are often scientifically uninteresting and exacerbate positivity violations. Modified treatment policies (MTPs) shift or otherwise transform the natural value of treatment (e.g. raise exposure by one standard deviation, subject to clamps).

Point-treatment continuous MTP. run_lmtp_grid estimates a δ-indexed curve under a user-chosen ShiftPolicy, with cross-fitted outcome and treatment nuisances and optional TMLE fluctuation. Density-ratio options (gaussian, classification, hybrid) implement practical continuous-exposure clever covariates in the spirit of the LMTP literature.

The identifying DAG for the linear synthetic DGP (simulate_linear_mtp) is baseline confounding of a continuous exposure:

using CausalDynamics, Graphs, DAGMakie, CairoMakie

g = DiGraph(3)
add_edge!(g, 1, 2)  # W → A
add_edge!(g, 1, 3)  # W → Y
add_edge!(g, 2, 3)  # A → Y
fig = plot_with_adjustment_set(g, 2, 3, [1]; node_labels = ["W", "A", "Y"])
fig

Sequential / multi-time LMTP. SequentialPolicy / run_sequential_lmtp implement a practical recursive outcome regression with a TMLE-style correction at $t = 1$, following the sequential identification strategy of Díaz et al. (2023). Numeric A_t share an additive ShiftPolicy. Categorical A_t take per-time DiscreteTreatmentPolicy values in policies (dummy-coded Q; Díaz–Williams classification ratio at $t = 1$). Mixed continuous/discrete treatments are rejected. Pair with CausalDynamics TemporalEffectQuery + unroll_temporal_dag → identify, then sequential_spec_from_identification / plan_sequential (or sequential_identification_certificate) so estimation carries an explicit ID certificate. estimand_from_query on a TemporalEffectQuery still builds a numeric LongitudinalPolicy by default; pass nonempty policies and wide treatments (:A1, :A2, …) for a factor SequentialPolicy. Observational panels from CausalDynamics simulate_panel use the same wide layout (baseline / timed :a1,:a2 / terminal :y) that SequentialPolicy expects; execute_estimand merges the certificate before running sequential LMTP. Latent/filter outputs enter via CausalDynamics ObservationBridge / panel_from_latent_series before the same sequential path.

Discrete-time survival / event-time LMTP. SurvivalPolicy / run_survival_lmtp estimate event-free probability at a horizon under a common MTP shift on time-ordered treatments, with at-risk sequential regression and optional thin IPCW for censoring. Pair with CausalDynamics TemporalEffectQuery (outcome = event-free indicator at the horizon) and survival_identification_certificate. Competing risks remain out of scope. Synthetic gate: simulate_discrete_survival_mtp.

Repeated outcomes under a static treatment. When treatment is fixed once and the same response is measured at several times (wide Y1…YT), run_repeated_outcome_msm estimates the unstructured profile $τ(t)=E[Y_t\mid do(A=1)]-E[Y_t\mid do(A=0)]$ with a joint influence-function covariance. Shared cross-fit propensity and per-time outcome regressions yield $\widehat\Sigma$ so contrasts such as $τ(t_3)-τ(t_2)$ use msm_contrast. Default estimator=:tmle fluctuates $Q_t$ on the shared clever covariate $A/g-(1-A)/(1-g)$; :eif is the untargeted one-step. Long (id, time, Y) tables pivot with unstack_repeated_outcomes. Missingness is a complete-profile policy (time_indexed=true); $Y_t$ are never imputed. This is the point-treatment / multi-outcome setting (Rosenblum & van der Laan 2010; R tmle::tmleMSM as a conceptual reference), not sequential LMTP (time-varying $A_t$). Pair with identify_repeated_outcomes (CausalTargeted helper over CausalDynamics TotalEffectQuery) when each $Y_t$ shares a backdoor set. Parametric MSMs project that profile (or stacked means $μ(t,a)$) onto a design with GLS: run_parametric_repeated_msm / ParametricRepeatedOutcomeMSM (:constant, :linear_time, :factor_time, custom matrices, or :mean_treatment_time with target=:mean). Pass cluster=:cluster (or a length-n id vector) for a cluster-robust sandwich on $\\widehat{\\Sigma}$; point estimates are unchanged (sampling hierarchy, not BLUP). Synthetic gates: simulate_repeated_outcome_ate, simulate_mean_treatment_time_msm. Generative nested $U$ DGPs live in CausalDynamics (RandomEffectSpec).

Gaussian MMRM (optional MixedModels extension). For a static treatment and repeated Gaussian outcomes in long (id, time, Y) form, fit_mmrm / run_mmrm fit outcome ~ treatment * time + baseline + (1 | id) (default) or an :unstructured random-effects approximation (1 + visit | id) with an internal categorical visit factor. mixed_g_computation supplies visit-specific contrasts through the common parametric standardisation layer (default random_effects=:zero). Each visit is standardised over its observed rows: with unbalanced follow-up, visit means may represent different empirical populations and use row weights rather than equal subject or cluster weights. This is a parametric trial-style reference beside LMTP/MSM, not a substitute for discrete longitudinal LMTP. The covariance choice in fit_mmrm describes the fitted random-effects structure; it is separate from g-computation's coefficient uncertainty. Requires the optional extension dependencies. Stress: docs/stress/mmrm_stress.qmd.

Transport weights. domain_transport_weights / transport_weighted_mean provide marginal IPTW-style domain reweighting after a CausalDynamics TransportQuery certificate. Decision. choose_policy evaluates labelled Estimands with execute_estimand and returns a PolicyChoice (max/min scalar TE).

Targeted learning, Super Learner, and cross-fitting

Adjustment covariates are encoded before they enter Super Learner. An internal fitted schema uses StatsModels default DummyCoding for string and categorical columns and preserves its levels, dummy-column meanings, and column order across all folds. Numeric values (including ordinary integer columns) remain one continuous column, while Bool becomes numeric 0/1. Missing-data handling (:drop, :impute, :ipcw, or :ipcw_impute) occurs before schema fitting; the schema itself does not impute or standardise.

Finite-support treatments use DiscreteTreatmentPolicy: T=1 via run_discrete_lmtp, multi-time via SequentialPolicy policies. Integer columns remain numeric on the continuous MTP path unless factorised. Random Forest, EvoTrees, and XGBoost receive numeric dummy columns rather than native categorical features; Random Forest therefore computes mtry from the encoded feature count.

Missing data. Full policy catalogue: Missingness. In brief: incompleteness is stratum × PCH rung (Structural $R$ in CausalDynamics; Dynamical sequential / survival gaps; Observable policies here). handle_missing_data supports :drop, :impute, :ipcw, and :ipcw_impute before schema fitting; IPCW weights enter IF summaries so :drop and :ipcw need not coincide. Survival censoring IPCW is not MAR missing $S_T$. Opt-in impute_posterior + run_lmtp_grid(...; imputation=) pools under a MAR certificate. Stress: missingness grid.

G-computation SEs. run_gcomp percentile intervals use a refitting bootstrap (redraw rows, recompute cross-fitted $Q$). Missing-data strategies also refit their response weights and covariate preprocessing on each resample of the original rows; bootstrap_refits_missingness records this behaviour. With estimator=:plugin, n_boot=0 returns the point estimate and unavailable inferential SE/limits (NaN). Its prediction_se measures fitted-contrast variation only. For a complete-data binary ATE with independent rows, estimator=:onestep adds the cross-fitted AIPW correction and uses the matching influence curve. It requires one outcome fit and, unless propensity is known, one propensity fit per fold, using the package's usual cross-validated binomial stack rather than inverse training-MSE weighting. A known scalar or row-aligned probability vector can describe randomised treatment assignment. Both avoid hundreds of bootstrap outcome-model refits. This is a distinct estimator; it does not silently replace plug-in g-computation.

Covariance and sampling units. Discrete static contrasts share outcome fits and calculate variance from paired influence contributions. Cluster IDs remain aligned after missingness handling; IPCW weighting precedes cluster aggregation. Cluster-aware folds keep dependent rows together. Bootstrap cluster copies receive distinct IDs and unsuccessful replicates are counted. A result with too few successful draws has status=:insufficient_success and unavailable limits, rather than a substituted point-fit interval. Staggered DiD aggregation retains shared-control and repeated-unit covariance; cell-only external tables explicitly report unavailable covariance.

Policy regularity. run_lmtp_grid retains its clamped-policy approximation, with inference_status=:approximate_clamped_policy on nonzero shifts. A clamp maps continuous tails to boundary masses, so its score limits are not ordinary nonparametric EIF intervals. Hard truncation is prespecified by default (cv_trunc=false); the opt-in variance-minimising heuristic is not cross-validation and does not account for truncation bias. The multiplier helper centres columns internally, but resampling cannot repair a biased or incomplete score. Continuous point-intervention contrasts return unavailable inferential SEs. Discrete sequential LMTP uses the full cumulative-ratio one-step score; continuous sequential/survival routes preserve point fits but report unavailable uncertainty until their policy and full score contract are supported. Mean-imputed covariates and missingness indicators are preprocessing, rather than a general causal recovery method. Grid fits that actually impute covariates retain point estimates but report :unsupported_covariate_imputation and unavailable inferential limits; this is distinct from posterior draws for missing outcomes under a declared Gaussian MAR model.

This fitted-schema path covers CausalTargeted's g-computation, LMTP Gaussian and classification/hybrid density-ratio nuisances, fold caching, sequential LMTP, survival LMTP, and missing-data nuisance models. CausalMediation reuses fit_covariate_schema / design_matrix for fold-stable string and categorical covariates when running run_mediation_grid.

Parametric regression standardisation. fit_parametric_gcomp fits an outcome model; gcomp_mean, gcomp_contrast, and gcomp_interaction select target rows, apply interventions to copies, predict on the response scale, and average over those rows. The default backend=:glm retains Gaussian/identity, binomial/logit, Gamma/log, and fixed or estimated NB2/log fits. The optional backend=:mixed uses supported Gaussian LMM and fixed- or estimated-shape NB2 random-intercept fits. run_parametric_gcomp fits and contrasts in one call. These are parametric references separate from cross-fitted run_gcomp and LMTP.

Training and target tables may have different rows, sizes, and order. The target need not contain an outcome column. Target selection precedes intervention; gcomp_interaction compares treatment effects within the observed modifier subgroups. Each selected target predictor must be complete-case, without silent row dropping or imputation. Interventions rebuild the fitted formula schema, including categorical contrasts, transformations, and interactions; unseen fixed-effect levels are rejected. Mixed grouping identifiers cannot be intervened on. A grouping factor such as site need not be a treatment-assignment unit: treatment may vary between individuals within a site. The static-treatment restriction remains in the repeated-outcome convenience interfaces.

For mixed predictions, random_effects=:zero sets latent effects to zero; :marginal integrates the response mean over the supported fitted Gaussian random-effects distribution. These agree for Gaussian identity-link models, including supported random coefficients. For NB2/log random-intercept models, :marginal multiplies exp(Xβ) by exp(σ_b²/2). It does not average fitted subject-specific effects. State the mode explicitly for new nonlinear mixed calls. GLM inference keeps covariance=:hc3 by default, with :model and :none available; mixed fixed-effect delta inference uses the fitted fixed-effect covariance. Ratios use log-scale intervals and multiple effects use their joint coefficient covariance. Delta inference holds fitted shape, dispersion, and variance components fixed. Marginalising predictions does not propagate their estimation uncertainty. The existing estimated-shape NB2 mixed_g_computation parametric bootstrap refits the hierarchical model; bootstrap_gcomp_interaction remains the GLM row-refit bootstrap. Unsupported mixed bootstrap requests fail explicitly. These procedures do not alter the cross-fitted inference of run_gcomp.

The same GLM workflow remains valid:

using CausalTargeted, DataFrames, StatsModels
training = DataFrame(Y=[1.0, 2.0, 2.5, 3.5, 1.3, 2.2],
                     A=[0, 1, 0, 1, 0, 1], X=[0.0, 0.0, 1.0, 1.0, 2.0, 2.0])
fit = fit_parametric_gcomp(@formula(Y ~ A + X), training; family=:gaussian)
gcomp_contrast(fit, select(training, Not(:Y)); treatment=:A,
               reference=0, comparison=1)

For a non-longitudinal mixed model, individuals can have different treatments within a site. Load the optional extension stack before fitting:

using CausalTargeted, MixedModels, FastGaussQuadrature, NLopt, SpecialFunctions
using DataFrames, StatsModels, Random, Distributions
rng = MersenneTwister(30)
training = DataFrame(site=repeat(1:16, inner=6),
                     A=repeat([0.0, 1.0, 0.0, 1.0, 0.0, 1.0], 16),
                     Group=repeat(["low", "high"], inner=3, outer=16),
                     X=randn(rng, 96))
site_effect = 0.4 .* randn(rng, 16)
training.Y = 1 .+ 0.6 .* training.A .+ 0.3 .* training.X .+
             0.2 .* (training.Group .== "high") .+
             site_effect[training.site] .+ 0.1 .* randn(rng, 96)
fit = fit_parametric_gcomp(@formula(Y ~ A * Group + X + (1 | site)), training;
                           backend=:mixed, family=:gaussian, id=:site)
target = select(training, Not(:Y))
gcomp_contrast(fit, target; treatment=:A, reference=0.0,
               comparison=1.0, random_effects=:zero)
gcomp_interaction(fit, target; treatment=:A, reference=0.0,
                  comparison=1.0, modifier=:Group,
                  modifier_reference="low", modifier_comparison="high",
                  random_effects=:zero)

For an NB2/log mixed random-intercept model, the two modes have different means. theta supplies a fixed shape; omitting it selects the dedicated estimated-shape fitter:

μ = exp.(0.2 .+ 0.4 .* training.A .+ 0.1 .* training.X .+
         site_effect[training.site])
training.Count = [rand(rng, NegativeBinomial(2.0, 2.0/(2.0+m))) for m in μ]
nbfit = fit_parametric_gcomp(@formula(Count ~ A + X + (1 | site)), training;
                             backend=:mixed, family=:negbin, theta=2.0, id=:site)
gcomp_contrast(nbfit, select(training, Not(:Count)); treatment=:A,
               reference=0.0, comparison=1.0, scale=:ratio,
               random_effects=:marginal)

The trial-style repeated-outcome front end remains available. Here arm is static within each subject; each observed visit contributes its actual rows:

panel = DataFrame(subject=repeat(1:24, inner=3), visit=repeat(0:2, 24),
                  arm=repeat(Float64.(isodd.(1:24)), inner=3),
                  W=repeat(randn(rng, 24), inner=3))
panel.Y = 1 .+ 0.5 .* panel.arm .+ 0.2 .* panel.visit .+
          0.1 .* panel.W .+ 0.3 .* randn(rng, 72)
run_mmrm(panel; outcome=:Y, treatment=:arm, time=:visit,
         id=:subject, baseline=[:W], random_effects=:zero)

Neither random effects nor parametric fitting identify a causal effect by themselves. Interpretation still requires an appropriate treatment and target population definition, exchangeability, positivity, consistency, and a suitable outcome model.

At **small n***, rich libraries overfit. recommend_run_options / adaptive_learners prefer lean GLM/mean stacks when n < 80, consistent with the Super Learner principle that the library must be *estimable at the sample size at hand. For binary propensity or outcome nuisances where probability calibration matters, fit_super_learner(...; family=:binomial) defaults to :nnloglik (R method.NNloglik). The discrete Super Learner is metalearner=:cv_selector (Phillips dSL / sl3 Lrnr_cv_selector). Pass a nested_sl_candidate in the library to nest an ensemble Super Learner inside that selector (eSL-inside-dSL; opt-in, not a new default). LMTP density-ratio classifiers stay on :invmse. Categorical outcomes use family=:multinomial (convex combination of class-probability matrices; sl3 loss_loglik_multinomial). Categorical treatments in LMTP use run_discrete_lmtp with Díaz–Williams classification density ratios, not a multinomial propensity. Arm contrasts (e.g. SS vs R) can use run_discrete_lmtp_contrast, which fits two static policies and differences the estimates (independent SE approximation). For binary presence outcomes (Y ∈ {0,1}), pass family_outcome = :binomial so outcome nuisances use logistic Super Learner fits rather than the default Gaussian location model:

run_discrete_lmtp_contrast(
    df, :arm, :infected;
    arm_hi = "SS", arm_ref = "R", levels = ["R", "SS", "SC"],
    baseline = [:weight],
    family_outcome = :binomial,
    learners_outcome = (:glm, :mean),
)

For nonnegative count outcomes (parasite totals, colony counts), use family_outcome = :poisson or :negbin. Nuisance models use log-link GLMs; profiled NB shape per cross-fit fold; TMLE targeting uses a log-link fluctuation on the fitted mean. The reported estimand is $E[Y \mid \mathrm{do}(A)]$ on the count scale (#36):

run_discrete_lmtp_contrast(
    df, :arm, :count;
    arm_hi = 1, arm_ref = 0, levels = [0, 1],
    baseline = [:weight],
    family_outcome = :negbin,
    learners_outcome = (:glm_nb, :mean),
)

When excess zeros are scientifically co-primary, prefer the two-part hurdle path below rather than a single count model.

The same family_outcome keyword is forwarded by run_lmtp_grid and run_discrete_lmtp. Multi-time factor recodes use the same policies on run_sequential_lmtp. Optional MLJ / MLP candidates are opt-in: they can improve recovery on some DGPs in a single synthetic draw while diluting others (overfitting vs generalisation). Prefer repeated Monte Carlo and library ablations before changing defaults.

Panel bridge (CausalDynamics → CausalTargeted)

Temporal identification in CausalDynamics yields adjustment certificates and plan_targeted_estimation selects a runner symbol and wide column names. CausalTargeted executes the plan with run_estimation_plan.

using CausalDynamics, CausalTargeted, DataFrames

spec = TemporalDAGSpec(
    nodes = [
        TemporalNodeSpec(
            :grid_type;
            temporal_support = FromOnsetSupport(0),
            value_representation = :attribute,
            causal_role = :assigned,
        ),
        TemporalNodeSpec(:fec; value_representation = :state),
    ],
    edges = [LaggedEdge(:grid_type, :fec, 0)],
)
u = unroll_temporal_dag(spec, 4)
query = TemporalEffectQuery(:grid_type, :fec, 2, 2)
plan = plan_targeted_estimation(
    u, query, names(wide);
    data = wide,              # optional: sets plan.estimability
    outcome_specs = Dict(     # optional: hurdle → :two_part_discrete_lmtp
        :fec => NodeOutcomeSpec(OutcomeKind.hurdle, :fec_bin, :fec_intensity),
    ),
)
res = run_estimation_plan(
    wide, plan;
    arm_hi = "SS", arm_ref = "R", levels = ["R", "SS", "SC"],
    learners_outcome = recommend_run_options(nrow(wide)).learners_outcome,
)

For dynamical queries (t_treat < t_outcome on the same variable), the planner selects :sequential_lmtp or :lmtp_grid; see run_lmtp_grid kwargs (lower_q, upper_q, shift is not a grid keyword). Stress notebook: apodemus_panel_stress.qmd is maintained outside the public repository because the cohort data are private.

Super Learner candidate roles

• :glm, :glm_interact, and :glm_quad provide ordinary, interaction-expanded, and quadratic-expanded GLMs (logistic under family=:binomial).
• :glmnet and its lasso/ridge aliases provide regularised linear candidates through MLJLinearModels.
• :randomforest provides a bagged Random Forest through MLJ and DecisionTree.jl. Its conservative defaults support stable nuisance estimation at modest sample sizes.
• :evotree and :evotree_deep provide Julia gradient boosting through EvoTrees.jl.
• :xgboost provides gradient boosting through MLJ and XGBoost.jl.
• :mlj_mlp and :mlj_nn_binary provide optional MLJFlux neural candidates.
• :mean is the intercept-only reference and fallback learner. It is not a standalone treatment-contrast engine: run_gcomp, run_lmtp_grid, and related entry points reject learners=(:mean,) with an ArgumentError. Always pair :mean with at least one treatment-dependent candidate (typically :glm).

The base DEFAULT_SL_LEARNERS remains (:glm, :mean). RICH_SL_LEARNERS adds :randomforest for bagging/model-class diversity, while :xgboost remains explicit opt-in because the rich library already contains EvoTrees boosting. SMALL_N_SL_LEARNERS and adaptive_learners remain conservative and do not add either flexible tree learner automatically.

using CausalTargeted, MLJ, MLJLinearModels
using MLJDecisionTreeInterface, MLJXGBoostInterface

fit = fit_super_learner(
    X,
    y;
    learners = (:glm, :glmnet, :randomforest, :xgboost, :mean),
)

Current Random Forest configuration:

Current XGBoost defaults:

Interventional mediation grids

Natural direct/indirect effects (Pearl, 2001; Robins & Greenland, 1992; VanderWeele, 2015) require cross-world assumptions that fail under intermediate confounding. Interventional (randomised interventional) effects (Vansteelandt & Daniel, 2017) and stochastic intervention mediation (Díaz & Hejazi, 2020; Hejazi et al., 2023) weaken those assumptions. Liu, Williams, Rudolph & Díaz (2024) unify modern mediation estimands with MTPs; the R package crumble (Liu et al., 2025 tutorial) is a software companion—Julia APIs use mediation names (run_mediation_grid, engine :mediation), with soft-deprecated run_crumble_* / :crumble aliases.

Implementation note. run_mediation_grid estimates TE / NDE / NIE under continuous MTP shifts via nested Monte Carlo and cross-fitted nuisances. Nested-MC variability is first-class: mediation_n_mc_sweep and mediation_stability_summary quantify SE and sign stability across n_mc (essential at small n).

A minimal mediation DAG (A → M → Y, A → Y) for interpreting those contrasts:

using DAGMakie, CairoMakie

fig, _ax, _p = dagplot_mediation(["A", "M", "Y"])
fig

Fold/δ cache. build_mediation_fold_cache (and the LMTP analogue) reuse outcome / mediator / exposure fits across δ within folds—same statistical estimand, lower wall time.

Positivity and support

Positivity (overlap) is necessary for identification of interventional means (Hernán & Robins, 2020; Petersen et al., 2012). MTPs are often designed so that shifted exposures remain in the support of the observed treatment law (Díaz et al., 2023).

Sensitivity to unmeasured confounding

Even with correct adjustment sets, estimates can tip under omitted confounding. CausalTargeted exposes diagnostic tipping-point and partial-R² calibrations inspired by Cinelli & Hazlett (2020); complementary classical tools include VanderWeele & Ding (2017) E-values and Rosenbaum (2002) sensitivity models.

Never silently replace a user DAG with a discovery graph in production defaults.

Outcome paths and model choice

Biological and trial endpoints differ in support, time structure, and whether zero inflation is scientifically primary. The table below maps common patterns to CausalTargeted engines. All targeted paths share cross-fitted nuisances, optional handle_missing, and cluster= sandwich SEs where noted. None of these replace identification: pair with CausalDynamics certificates and the query ladder (association → intervention → counterfactual).

Two-part hurdle LMTP. Build presence ($I(Y>0)$) and intensity (e.g. log(y) on positives, with missing when $Y=0$). Presence uses binomial outcome nuisances; intensity fits on the positive subsample only, so the reported contrast is conditional on observed infection, not the marginal $E[Y \mid \mathrm{do}(A)]$ unless you define a composite estimand separately.

res = run_two_part_discrete_lmtp_contrast(
    df, :grid_type;
    presence = :fec_bin,
    intensity = :fec_logpos,
    arm_hi = "SS", arm_ref = "R", levels = ["R", "SS", "SC"],
    baseline = [:weight],
    family_presence = :binomial,
    cluster = :mouse_id,
)
res.presence.estimate    # co-primary: infection risk difference
res.intensity.estimate   # secondary: log-burden among positives

Use suggest_family_outcome or recommend_run_options(n; outcome=col) to set family_outcome for single-outcome runs. For DAG faithfulness on hurdle nodes, see the next section; for repeated Gaussian visits, compare MSM (IF) with MMRM (parametric) on the same panel before interpreting biology.

Hurdle CI testing

Zero-inflated biological outcomes are often operationalised as a two-part hurdle: presence $I(Y > 0)$ and log intensity among positives. Gaussian partial correlation on $\log(1 + Y)$ is misspecified for faithfulness checks on those nodes.

Use partial correlation for continuous block variables. Use test_implied_hurdle_independences when a statement involves a hurdle-split node (binomial presence GLM; Gaussian intensity GLM among positives). String / categorical predictors (e.g. grid_type) use StatsModels DummyCoding, aligned with discrete LMTP (#34) and two-part LMTP (#35).

results = test_implied_hurdle_independences(
    statements, df,
    Dict(:fec => (:fec_bin, :fec_intensity), :eimeria => (:eimeria_bin, :eimeria_intensity)),
)

Identification certificates (CausalDynamics bridge)

Estimation attaches provenance via identification_certificate / attach_run_metadata!. Upstream ID uses Pearl’s do-calculus toolkit (Pearl, 2009; Shpitser & Pearl, 2006) and g-methods (Robins, 1986; Robins, 2000). Time-indexed queries use TemporalDAGSpec / unroll_temporal_dag / TemporalEffectQuery in CausalDynamics (see that package’s References).

Small-n profile

Conservation biology, ecology, and early trials often have tens to low hundreds of units. recommend_folds, SMALL_N_SL_LEARNERS, and recommend_run_options encode memory-safe, positivity-aware defaults (parallel=false, higher mediation n_mc when n is small). See the Small-n checklist.

What we deliberately do not claim

• Full parity with every option in R lmtp / crumble (GPU Riesz nets, all mediation estimand flavours, competing-risks survival LMTP).
• That tipping-point / partial-R² helpers replace design-based identification.
• That Super Learner at n ≈ 30 recovers oracle rates—diagnostics exist precisely because they often do not.