DoubleML Pipeline

A robust Python library for automated causal inference utilizing Double Machine Learning (DoubleML) and FLAML AutoML.

By Google's APAC Marketing Effectiveness Taskforce

Note: doubleml-pipeline refers to and builds upon the underlying DoubleML Python package published by doubleml.org.

---

Why DoubleML Pipeline?

With DoubleML Pipeline, you can estimate the true incremental impact and ROI of your marketing interventions including pricing and promotions by controlling for complex, high-dimensional confounders.

Causal inference is challenging due to confounding bias, where standard regression models often fail to isolate true causal relationships. Furthermore, unlike standard supervised learning, Causal Machine Learning (CausalML) lacks observational ground truth data for unobserved counterfactuals, making verifiability a critical requirement for decision-making.

Double Machine Learning (DML) was selected as our base framework because it solves these challenges through rigorous theoretical and statistical guarantees—such as Neyman orthogonality and cross-fitting to eliminate bias—allowing causal estimates to be mathematically verified and trusted. DoubleML Pipeline operationalizes this framework by combining DML with automated machine learning (AutoML) to make causal estimation both scientifically robust and highly automated.

Key Features:

• Automated Nuisance Modeling: Automates nuisance parameter estimation using FLAML AutoML. It dynamically searches and tunes candidates (like LightGBM, XGBoost) under a specified time budget, removing the guesswork of model selection.
• Safe Business Unit Scaling: Normalizes and scales geographical and population data securely, ensuring that descaled results represent accurate real-world business metrics (KPI, ROI) without introducing mathematical singularities.
• Treatment Effect Heterogeneity (CATE): Automatically fits continuous splines or binary bases to estimate how treatment effects vary across different contexts (Conditional Average Treatment Effects).
• Sensitivity Analysis (OVB): Runs robust omitted variable bias (OVB) checks to evaluate how sensitive your causal estimates are to unobserved confounders. The bias bound is rho_U * scale * sqrt(cf_y * cf_d / (1 - cf_d)) with the residual-based scale = sqrt(sigma2 * nu2) (DoubleML's sensitivity_elements; recomputed from the out-of-fold residuals when unavailable; per-group moments for geo / item GATE rows) and the robustness value is RV = (sqrt(f^4 + 4 f^2) - f^2) / 2 with f = (|theta| - 1.96 se) / scale (Cinelli & Hazlett 2020), so the RV sits on the zero contour of the 95% lower bound; sensitivity_metrics.csv carries n_obs, f_q, sigma2, nu2, scale so the numbers can be verified. With a saturation family adopted in Phase 0 the analysis runs on the same PLR fitted with D = fbar(T; alpha_hat) (no extra fit), the benchmark covariate is chosen by correlation with fbar(T), and the bounds are conditional on alpha_hat (treatment_variable = "f_bar(T)" in the CSV).
• Bayesian MMM Prior Calibration: Extracts fitted distribution parameters (Normal, LogNormal, Beta) for treatment contribution (%) and ROI (e.g., discount) as prior knowledge. Because CATE estimated by DML fluctuates depending on conditions (e.g., seasonality), using only the mean of the distribution as input to Bayesian MMM priors is a consideration. Additionally, it is important to be aware that there may be a gap between ATE MMM estimates and CATE DML estimates due to differences in the estimands.
• Publication-Ready Visualizations: Generates comprehensive time series line plots, stacked contribution bar charts, cross-sectional ROI comparisons, and prior distribution shapes.

---

Getting Started

Prerequisites

Before running the pipeline, ensure you have:

• Python 3.10 or later.

Installation

You can install the library directly from PyPI:

pip install doubleml-pipeline

Alternatively, to install from source:

git clone https://github.com/google/doubleml-pipeline.git
cd doubleml-pipeline
pip install .

Quick Start Example

Here is how to initialize and run the causal pipeline in a standard Python environment, explicitly configuring all scaling, modeling, and financial parameters.

import pandas as pd
# Note: Module name uses underscores (_) instead of hyphens (-)
from doubleml_pipeline.preprocessing.scalers import CausalDataScaler
from doubleml_pipeline.workflow import CausalWorkflowOrchestrator

# 1. Load your dataset
df = pd.read_csv("path/to/your/marketing_data.csv")

# 2. Initialize the scaler with explicit scaling strategies per column
scaler = CausalDataScaler(
    standardize_cols=["holiday_index", "seasonality_factor"],
    population_standardize_cols=["competitor_spend"],
    population_median_normalize_cols=["sales_revenue", "tv_spend"],
    min_max_cols=["store_count"],
    median_normalize_cols=["discount_rate"],
    population_col="geo_population",
)

# 3. Initialize the orchestrator with all explicit configuration parameters
orchestrator = CausalWorkflowOrchestrator(
    # Primary target and treatments (single treatment to support optimization)
    y_col="sales_revenue",
    d_cols=["discount_rate"],
    x_cols_list=[
        ["holiday_index", "seasonality_factor", "competitor_spend"],
    ],
    # Nuisance estimation models & cross-validation setup
    ml_models_y=["lgbm", "xgboost"],
    ml_models_t=["lgbm", "xgboost"],
    n_folds_list=[5],
    output_dir="./causal_results",
    n_jobs=-1,
    # Conditional Average Treatment Effect (CATE) specification
    covariates_for_cate_list=[
        ["seasonality_factor"],
    ],
    cate_structure_list=[
        {"type": "auto_additive", "df": 3, "degree": 2, "include_intercept": False},
    ],
    time_budget=30,
    # Panel data dimensions
    date_col="date",
    geo_col="geo",
    item_col="item",
    geo_name_col="geo_name",
    item_name_col="item_name",
    # Sensitivity analysis (Omitted Variable Bias)
    run_sensitivity=True,
    sensitivity_scope=["total", "geo", "item"],
    # Ground truth validation (set True if ground truth effect columns exist)
    ground_truth_existence=False,
    ground_truth_effect_column=None,
    # Treatment type classification & monetary spend mapping
    treatment_types={
        "discount_rate": "percentage",
    },
    treatment_spend_cols={
        "discount_rate": "discount_cost",
    },
    sales_col="sales_revenue",
)

# 4. Execute the full 4-phase pipeline (Grid Search -> Shortlist -> Ensemble -> Consolidation)
orchestrator.run_full_pipeline(
    df=df,
    scaler=scaler,
    top_n=3,
    n_reps=10,
)

The orchestrator will execute the models in parallel, select the best estimators, run repetitions to calculate robust confidence bounds, perform sensitivity checks, and write all dataframes and plot artifacts to the specified ./causal_results directory.

Pipeline Execution Flow (4 Phases)

1. Grid Search Exploration: Fits nuisance models for outcome Y and treatment T across all candidate AutoML models and fold combinations in parallel.
2. Shortlisting: Ranks and selects top-performing models (top_n) based on combined error metrics.
3. Robust Ensemble Estimation: Re-estimates shortlisted models across n_reps iterations to build a weighted ensemble with empirical confidence bounds (2.5%–97.5%).
4. Consolidation & Reporting: Aggregates causal effects (CATE, incremental KPI, ROI), runs OVB sensitivity checks, and outputs charts.

Output Directory Structure

When the pipeline runs, it saves all intermediate evaluation metrics, candidate model repetitions, consolidated datasets, and visualization plots within the specified output_dir (e.g., ./causal_results/) as shown below:

Important Notes on Outputs & Optimization

1. ROI Charts: The ROI chart folders under roi_output/<treatment>/ (003_roi_by_month, 004_roi_by_geo, 005_roi_by_item) are conditionally generated depending on your configured treatment_types (default is "spend"), along with treatment_spend_cols and sales_col, as summarized below:

   treatment_types	treatment_spend_cols & sales_col	ROI Charts	Pipeline Behavior & Reason
   "spend" (default)	Not required	Not Generated	Not generated specifically because they are redundant with the output generated in Phase 4.
   "percentage", "impressions", or "price"	Provided	Generated	Explicitly generates the ROI charts (003_roi_by_month, 004_roi_by_geo, 005_roi_by_item).
   "percentage", "impressions", or "price"	Missing	Skipped (Warning)	Pipeline logs a warning and skips ROI charts because monetary spend mapping is missing.

2. Optimization: Optimization is limited to a single treatment. Two modes are available:

   optimize	simple_optimization	saturation block	Behavior
   True	True (derived when no family is set)	Not required	Linear (simple) optimization, unchanged from previous releases.
   True	False (derived when a family is set)	Required	Saturation-aware optimization (see below).
   True	False (passed explicitly)	Missing	ValueError at construction. Detailed optimization is only defined under a saturation model.
   True	True (forced for non-spend treatment types)	Optional	Simple optimization. The saturation-aware optimizer is only implemented for treatment_types = "spend".

3. Cluster DML keys (chosen automatically, every run): The cross-fitting folds and the cluster-robust standard errors use a 2-way Cluster DML whose keys are picked by the cluster-key determination rule — this runs whether or not saturation_config is given. The candidates are item_col, geo_col (each needs at least two distinct values) and a contiguous date block derived from date_col (never raw dates); the two candidates with the most clusters are adopted, one if only one is usable, none → observation-level inference. date_col, item_col and geo_col keep their visualization role in every chart and summary CSV regardless of which keys were adopted. The choice is logged (at the start of Phase 1 in a normal run, in the Phase 0 header otherwise):

   Cluster DML keys (determination rule applied automatically; no saturation_config needed): 2-way: cluster_date_block (7 clusters), geo (6 clusters); not adopted: item (2 clusters); date block = monthly blocks (weekly data)
   

   and recorded in metadata.json["cluster_dml"] (keys, n_way, cluster_counts, not_adopted, date_block_granularity, visualization_cols, keys_used, …); the HTML report shows the same information in its header. See docs/output_csv_schema.md.

4. Hyperparameter policy (hparam_policy, default "global"): Phases 1–3 fit a nuisance learner in every cross-fitting fold of every candidate and every repetition. With hparam_policy="global" the FLAML search (time_budget seconds, eval_method="auto", early_stop=True) runs once per (treatment, learner) at the start of Phase 1 — on the sample Phases 1–3 are fitted on (final_fit), with target Y for the Y-side learners and the PLR treatment D (f̄(T; α̂) when Phase 0 adopted a saturation family, T otherwise) for the T-side learners — and every fold then fits the winner's concrete estimator with those fixed hyperparameters (models.fixed_hparam_regressor.FixedHparamRegressor). The DML cross-fitting, CATE, repetitions, median aggregation and the combined error (now fed by the same validation RMSE measured inside each fold's training data with the method FLAML chose) are unchanged. The searches are persisted to step1_exploration_<d_col>/hparams_<d_col>_<learner>.json (one file per learner; a learner used on both sides has a Y and a T entry). hparam_policy="per_fold" restores the previous behaviour (a FLAML search inside every fold). Independently of the policy, Phase 1 reuses its results as repetition 0 of Phase 3, Phase 3 runs all (candidate, repetition) jobs as one parallel job list, and the parallel sections allocate workers = min(tasks, cores) processes with cores // workers learner threads each (logged as workers = W, threads/worker = h). The sensitivity analysis fits the benchmark's short model once per treatment and reports group-specific benchmark points (benchmark_cf_d, benchmark_cf_y, plus the new benchmark_rho and benchmark_scope columns of sensitivity_metrics.csv). Phase 0 is not affected by hparam_policy.

---

Saturation DML (optional)

By default the pipeline is unchanged: omit the saturation_config argument and Phases 1–4 behave exactly as before. Setting saturation_config['saturation']['family'] enables Phase 0, which estimates a diminishing-returns parameter alpha for each treatment and then feeds the transformed treatment f(T; alpha) into the existing DML pipeline.

Response families

family	f(T; alpha)	Interpretation of alpha
"tanh" (default)	tanh(alpha * T / 2)	Rate; T_1/2 = ln 3 / alpha is the half-saturation point.
"exp"	1 - exp(-alpha * T)	Rate; 1/alpha is the saturation scale.
"log"	log(1 + T / alpha)	Saturation scale (smaller = returns diminish sooner).
"power"	T ** alpha	Elasticity; must satisfy 0 < alpha < 1.

Configuration

Every saturation setting lives in a single saturation_config dict, and the presence of the saturation block is all that is required (family defaults to "tanh") — everything else defaults. Turning it on is therefore a one-key change:

orchestrator = workflow.CausalWorkflowOrchestrator(
    # ... all existing arguments ...
    optimize=True,
    saturation_config={"saturation": {"family": "power"}},
)

simple_optimization is derived: it becomes False automatically when a family is set, and True otherwise. Pass it explicitly to override.

The full set of keys, all shown at their defaults:

orchestrator = workflow.CausalWorkflowOrchestrator(
    # ... all existing arguments ...
    optimize=True,
    saturation_config={
        "saturation": {
            "family": "tanh",             # "tanh" (default) | "exp" | "log" | "power"
            "grid": {                     # alpha grid searched in Phase 0
                "mode": "auto",           # "auto" | "manual"
                "num": 25,
                "c_min": 0.25,            # auto: T_1/2 in [c_min*q50, c_max*q95]
                "c_max": 4.0,
                "max_extensions": 2,      # auto: x4 extension at an edge argmin
                "min": None,              # manual only (alpha units)
                "max": None,              # manual only (alpha units)
                "scale": "log",           # "log" | "linear" (manual spacing)
            },
            "basis_dim_T": "auto",        # "auto" | int (B-spline T-basis size J)
            "basis_tol": 1e-3,            # auto: rel. RMS approximation error
            "basis_dim_range": [6, 20],   # auto: search range for J
            "ridge": 1e-4,                # ridge penalty in the closed form (absolute)
            "t_ref_quantile": 0.95,       # T_ref = q95(T>0): all candidates pass through (T_ref, 1)
        },
        "alpha_selection": {
            "sel_fraction": 0.20,         # D_sel : D_inf split
            "n_seeds": 10,                # S: repetitions of R-loss profiling
            "k_folds": 5,                 # K for the alpha-selection stage
            "refine": "brent",            # "brent" | "none" (grid + Brent only)
            "linear_test_k": 1.0,         # adopt saturation iff d_bar < -crit * se (se = two-way cluster-robust se_obs)
            "linear_test_t_critical": False,  # crit = k (False) | Student-t value t_{G-1}^{-1}(Phi(k)) (True)
            "fallback_half_saturation_T": None,  # half-saturation point (units of T) used when Diag1/Diag2 fires
            "flatness_ratio": 10.0,       # Diag2: fires when max/min T_1/2 over the 1-SE set exceeds this
            "t_ref_tolerance": 1.0e-9,    # Diag6: rel. error tolerated between Phase 0 and Phase 1 T_ref
            "learners": {                 # Stage B nuisance learners (fixed)
                "learner_selection": False,   # True: legacy FLAML candidate search
                "y_model": {"class": "xgboost.XGBRegressor",
                            "params": {...}, "early_stopping_rounds": 30},
                "t_model": {"class": "xgboost.XGBRegressor",
                            "params": {"multi_strategy": "multi_output_tree",
                                       ...},
                            "early_stopping_rounds": 30},
                "small_sample_threshold": 2000,  # n_sel below -> RidgeCV
                "small_sample_model": "ridge_cv",
            },
            "output_dir": "step0_saturation",  # under output_dir
            # "time_budget": 30,          # learner_selection only
        },
        "final_fit": {"data": "d_inf"},   # "d_inf" (default) | "full_sample"
        # Saturation-aware optimization. Note this lives
        # under `optimization`, not `saturation`.
        "optimization": {
            "budget_mode": "fixed",       # "fixed" (budget K) | "roi_target" (no budget; ROI >= r)
            "budget_increment": 0.0,      # fixed only: K = (1 + x) * actual total of the target units
            "budget_reallocate": True,    # fixed only: redistribute what the ROI cap / negative units free up
            "cap_quantile": 1.0,          # support cap T_bar_t = max(Q_q(cell, T>0), T_actual); 1.0 = cell maximum
            "reallocate_within": "all",   # "all" (one KKT problem; spend moves across days) | "date" (one per date)
        },
    },
    # top-level optimizer scope (None = all rows of the column).
    opt_geos=None,
    opt_items=None,
    opt_threshold_roi=1.0,                # r >= 0 (default 1.0 across all optimizer paths); < 0 raises at initialization
)

The keys in detail:

Key	Default	Meaning
saturation.family	"tanh"	Used when the saturation block is present but family is omitted. family: None in the block still means a linear run.
saturation.grid.mode	"auto"	auto builds a dimensionless half-saturation grid T_1/2 in [c_min*q50(T>0), c_max*q95(T>0)], log-spaced, for tanh/exp/log. power always uses a linear alpha grid. manual uses min/max/scale in alpha units. A legacy grid with min/max but no mode is treated as manual.
saturation.grid.c_min, c_max	0.25, 4.0	Auto-grid range multipliers.
saturation.grid.max_extensions	2	When the argmin is at an auto-grid edge, extend that side x4 in T_1/2 (matrix algebra only, no ML refits).
saturation.basis_dim_T	"auto"	auto scans J in basis_dim_range and picks the smallest J whose worst-case relative RMS error of projecting fbar(.;alpha) is below basis_tol. The projection is a least-squares fit (SVD pseudo-inverse, never B'B) on a fixed grid over [0, T_max] made of 500 equally spaced points, 201 quantiles of T>0 and the Greville abscissae of the basis, so the basis matrix B has full column rank however skewed T is; T_max is the 99.9 % quantile of T>0 on the full fitting data (rows above it are clipped by the basis only, never by Phases 1-4). Bases with condition number cond(B) >= 1e8 are skipped; if no J meets basis_tol, the J with the smallest error is used (Diag1 warns). An integer pins J. A pinned or only-available basis with cond(B) >= 1e8 makes Stage B adopt the linear model (decision_reason = "basis_ill_conditioned") instead of trusting it.
saturation.basis_tol, basis_dim_range	1e-3, [6, 20]	See above.
saturation.t_ref_quantile	0.95	Quantile of T>0, computed once on the full fitting data before the split, that defines T_ref. Every candidate is normalized to fbar(T; alpha) = f(T; alpha) / f(T_ref; alpha), so all curves pass through (T_ref, 1) and the identity candidate is fbar_T = T / T_ref. This removes the alpha-dependent magnitude of f that let the ridge solve shrink near-identity candidates (large T_1/2) to zero. Always on.
alpha_selection.linear_test_k	1.0	Saturation is adopted only if d_bar < -crit*se, where d_s = L_s(alpha_hat) - L_s,T is the paired difference against the identity candidate fbar_T = T / T_ref and se is the cluster-robust se_obs of the per-observation differences d_i. With two adopted cluster keys (e.g. geo × date block) the se is two-way (Cameron–Gelbach–Miller: V = V_key1 + V_key2 - V_cell), which covers dependence along both keys; with a single key it is one-way on that key (there is no across-seed se). The one-way se on the fewer-cluster key is always reported next to it (se_one_way), and the decision line notes when it would have decided differently. crit = k by default (see the next row). Otherwise selected_family = None (linear). The adoption rate is tuned with k alone: k = 1.0 is the 1-SE rule, 2.0 a conservative choice, smaller values adopt saturation more readily. Independently of this test, an alpha_hat still on the upper grid edge after max_extensions extensions forces selected_family = None with decision_reason = "grid_edge_upper" (the data cannot distinguish the curve from a line); a lower-edge hit only warns.
alpha_selection.linear_test_t_critical	False	True replaces crit = k by the Student-t value with the same one-sided tail probability, crit = t_{G-1}^{-1}(Phi(k)), G = number of clusters (the smaller of the two counts for the two-way se). It compensates the noisy se of few clusters (k = 1: G = 30 -> 1.02, 10 -> 1.06, 4 -> 1.20) and can only make linear decisions more frequent; it is not applied to the Diag2 1-SE set. Off by default.
alpha_selection.fallback_half_saturation_T	None	A fixed half-saturation point in units of T (T_knee for log, alpha itself for power; same scale as half_saturation_T / T_ref in summary_alpha_selection.json, i.e. T / c). When Diag1 or Diag2 fires it replaces the estimate downstream and keeps the family; otherwise the 1-SE decision stands.
alpha_selection.flatness_ratio	10.0	Diag2 fires when T_1/2^max / T_1/2^min over the 1-SE set A_1SE = {alpha : L_bar(alpha) - L_bar(alpha_hat) <= se(alpha)} exceeds this (T_knee for log; for power the fixed rule alpha^max - alpha^min > 0.4). Warning only when a family is adopted, INFO otherwise.
alpha_selection.t_ref_tolerance	1e-9	Diag6: at the start of Phase 1 T_ref is recomputed from the treatment column Phase 1 receives; a relative error at or above this stops the run with an error (unit / hand-over bug).
optimization.budget_mode	"fixed"	fixed: a budget K is given and the allocation is solved in one KKT pass. roi_target: no budget; every unit receives min(T_bar_t, T_t^r) (requires opt_threshold_roi > 0).
optimization.budget_reallocate	True	fixed only. True (default): K' = K, solve once with the effective caps min(T_bar_t, T_t^r) (spend freed by negative-ROI units and the ROI cap is redistributed). False: K' = K_+ (negative-ROI spend removed), solve with the caps T_bar_t, then cut to T_t^r without redistribution (unspent_budget).
optimization.cap_quantile	1.0	q of the support cap T_bar_t = max(Q_q(cell_t), T_t^actual): Q_q is the q-quantile of the positive actual treatment within the unit's geo × item cell (consistent with T_ref = q95(T>0); zero-spend rows do not shift it), floored at the unit's own actual spend. 1.0 (default) is the cell maximum — today's behaviour, where one outlier day lifts the cap of every unit in its cell. With e.g. 0.95 the other units are capped at the 95 % point while the outlier day keeps its own spend as its cap: it is never increased, but may still be reduced by the KKT condition or the ROI cap like any unit. Must be in (0, 1].
optimization.reallocate_within	"all"	Scope of the KKT problem. "all": one problem over every allocation unit (spend moves across days, geos and items; today's behaviour). "date": one problem per date with the per-date budget K_d = (1 + budget_increment) * sum_{t in d} T_t^actual — no cross-day movement; mu* becomes per date (median with marginal_roi_at_optimum_min / _max in optimization_metrics.csv). See the no-carryover note under "Saturation-aware optimization".
opt_threshold_roi (top level)	1.0	ROI floor r >= 0, applied per row across both simple_optimization=True and False (including budget_mode='roi_target').
opt_geos, opt_items (top level)	None	Restrict the allocation units of the saturation-aware optimizer to these geos / items (None = all rows of the column). Rows outside keep their actual values in the optimized_* columns.
alpha_selection.learners	see above	Fixed Stage B learners. learner_selection: True restores the legacy FLAML search. Without XGBoost >= 2.0 a warning is emitted and scikit-learn stand-ins are used: HistGradientBoostingRegressor for Y and one native multi-output RandomForestRegressor for T.
alpha_selection.output_dir	"step0_saturation"	Sub-directory of output_dir that receives the Phase 0 files.

What Phase 0 adds

1. Top-level split into D_sel (alpha selection) and D_inf (estimation). With two cluster keys, the split unit is their intersection cell (e.g. geo × date block). No cell straddles the two halves, and every level of each key appears in both halves and, when it has at least K cells in D_sel, in every Stage B fold (others_split_coverage.csv).
2. Alpha selection by minimizing the profiled R-loss over the half-saturation grid, with optional Brent refinement. The same CATE basis p(X) is used as in Phases 1–3. Every nuisance prediction is out-of-fold (cross-fitting: each (seed, fold) model predicts only its held-out rows) and the inner CATE coefficients of fold k are fitted on the out-of-fold residuals of the other folds, so no in-sample residual enters the loss. The intercept-only curve is drawn for visualization only.
3. Identity-candidate comparison. All candidates are normalized through the reference point (T_ref, 1), T_ref = q95(T>0): fbar(T; alpha) = f(T; alpha) / f(T_ref; alpha). The saturated fit must beat the identity candidate fbar_T = T / T_ref by more than crit = linear_test_k cluster-robust standard errors (two-way on the two adopted keys when both are known; optionally a Student-t critical value with linear_test_t_critical). Otherwise selected_family = None. An alpha_hat left on the upper grid edge after all extensions also forces None (decision_reason = "grid_edge_upper"). 3b. Diagnostics Diag1–Diag6. Diag1: basis approximation error eps(J) on the final grid (warn if >= basis_tol; J is fixed after the initial grid) and basis conditioning cond(B) (warn if >= 1e8; the basis is then not trusted and the linear model is adopted, decision_reason = "basis_ill_conditioned"); the single Diag1 line also shows T_max (99.9 % quantile of T>0) and the number of D_sel rows the basis clips at it. Diag2: 1-SE flatness, the ratio of the largest to the smallest half-saturation point in the 1-SE set (warn if > flatness_ratio and a family is adopted; INFO otherwise). Diag3: r_T = sd(T~)/sd(T/T_ref) and Diag4: the shape quantiles fbar(q_p(T>0); alpha_hat) with their linear references q_p / T_ref (supporting values, no rule). Diag5: kappa(alpha_hat) = f(T_ref; alpha_hat) for tanh / exp (no rule). Diag6: unit consistency of T_ref at the start of Phase 1 (error stop). Diag1 / Diag2 warn and continue; the fixed fallback_half_saturation_T, if set, is used downstream when they fire.
4. If saturation is adopted, alpha is treated as a constant and Phases 1–3 run unmodified on fbar(T; alpha) (so the estimated theta is the effect of one T_ref worth of treatment). Otherwise they run on the raw treatment, and Phase 4 falls back to the linear optimizer.
5. Every row is scored. With the default final_fit.data = "d_inf", Phases 1–3 are fitted on D_inf only, but the fitted nuisance models and the CATE function are then applied to the held-out D_sel rows as well. Everything under step4_consolidated_results/ — the consolidated CSV, the CATE / incremental-KPI / ROI summaries, the charts and the Phase 4 optimization — therefore covers all rows, not just D_inf. The consolidated CSV carries an estimation_sample column (D_inf = fitted row / D_sel = scored row), the 1c metrics report number_of_fitted_rows and number_of_scored_rows, and the log states e.g. Consolidated frame covers all 180 rows: 142 fitted (D_inf) + 38 scored with the fitted CATE (D_sel). The OVB sensitivity analysis keeps using the fitted rows only.

Phase 0 outputs

All written to output_dir/step0_saturation/. With several treatments each one gets its own step0_saturation/<treatment>/ directory. Every file name carries a prefix that says what it belongs to: diag1_ … diag6_ for the six identification diagnostics, summary_ for the result of the alpha selection, others_ for supporting material without a diagnostic rule. Every PNG ends with a footnote that explains its annotations, states the rule / status of its diagnostic for this run and names the CSV columns behind it.

File	Contents
summary_alpha_selection.json	family, selected_family, alpha_hat, half_saturation_T (T_knee for log), alpha_grid_argmin, rloss_min (L_bar(alpha_hat)), refine_method, d_bar, se (cluster-robust se_obs), linear_test_k, decision, decision_reason (1se_pass / 1se_fail / grid_edge_upper), diag1_fired, J, basis_approx_error, diag2_fired, diag2_ratio, r_T, fbar_q25 … fbar_q90, linear_ref_q25 … linear_ref_q90, saturation_at_T_ref (null for log / power), grid_mode, grid_extensions, grid_edge_hit (none / upper / lower), T_ref, T_max, cluster_keys, cluster_fallback_tier, n_sel, n_inf, n_seeds, k_folds, learners, plus extras (for example phase4_optimization, diag6_status, flatness_ratio, basis_tol, fallback_half_saturation_T_used, t_ref_quantile, T_ref_raw, grid_variable, and the divisor of the normalized treatment: normalization, divisor_constant, divisor_population_col, divisor_formula, divisor_description). NaN is written as null. Phases 1–3 and Phase 4 read T_ref from this file.
summary_alpha_grid.csv	One row per grid point: grid_index, half_saturation_T, alpha, mean_r_loss (L_bar(alpha)), se_alpha (the Diag2 se(alpha)), in_1se_set (membership in A_1SE), fbar_q25, fbar_q50, fbar_q75, fbar_q90 (phi_p(alpha) = fbar(q_p(T>0); alpha)), basis_approx_error_alpha (the Diag1 eps(J) per alpha). This is the complete data behind the R-loss profile charts (it replaces the former alpha_profile_<treatment>.csv).
summary_rloss_curves.csv	seed, grid_index, alpha, rloss_cate, rloss_ate, rloss_identity.
summary_alpha_profile_<treatment>.png	The one-page Phase 0 summary, 2 + 2 + 1 panels from top to bottom: Grid construction (histogram of T>0 with q50, q95 = T_ref, the automatic grid range and every grid point) and Seed-wise alpha distribution (per-seed alpha_hat_s histogram); Profiled R-loss curve with A_1SE and Fitted saturation curve (normalized through (T_ref, 1)) — fbar(T; alpha_hat), the identity candidate T / T_ref, (T_ref, 1), the shapes of the smallest / median / largest grid alpha and the curve Phases 1–4 actually use labelled FINAL; and the full-width Key stats & process metrics table whose diagnostics column tags each row with diag1 … diag6. The footnotes define T_ref, q05 / q50 / q95, the grid ends 0.25*q50 (fastest saturation) / 4*q95 (mildest, practically linear) and list every file of the folder.
summary_alpha_seed-wise_alpha_distribution.csv	The data behind the seed-wise alpha distribution panel: seed, alpha_hat_s, rloss_at_alpha_hat_s, rloss_identity_s, d_s (visualization only, no rule).
summary_alpha_selection.csv, summary_phase0_diagnostics.json	Run-level: one row per treatment (decision, alpha_hat, Diag1–Diag5 flags, compute accounting) and the Stage A split / cluster-key plan with every Phase 0 warning.
diag1_basis_approx_error.png	Diag1: eps(J)(alpha) per grid point against basis_tol, with J, max eps(J) and the fired status in the title.
diag2_rloss_profile.png	Diag2: R-loss vs T_1/2: CATE curve (±1 se), intercept curve (visualization only), identity candidate L_bar_T (±k se), the 1-SE set A_1SE, alpha_hat, decision; the footnote states the flatness ratio against flatness_ratio and whether Diag2 fired.
diag3_identification_diagnostics.csv	Diag3: r_T, sd_T, sd_T_resid, share_T_zero, q50_T_pos, q95_T_pos, T_max, T_ref.
diag4_f_transform_quantiles.png	Diag4: the four shape quantiles phi_25, phi_50, phi_75, phi_90 over alpha with the linear references q_p / T_ref.
diag5_saturation_at_T_ref.png	Diag5 (tanh / exp only): kappa(alpha) = f(T_ref; alpha), the share of the ceiling reached at T_ref, over the grid with kappa(alpha_hat) marked.
diag6_t_ref_consistency.csv	Diag6, written at the start of Phase 1: treatment, status (pass / fail / skipped when selected_family = None), T_ref_phase0, T_ref_phase1, rel_error, t_ref_tolerance, t_ref_quantile, n_rows_phase1, note.
others_split_coverage.csv, others_split_coverage_detail.csv, others_split_coverage.png	Stage A / Stage B bookkeeping: per key level n_cells_sel, n_cells_inf, covered, n_folds_covered_min, fold_covered (the summary CSV, drawn as the bar chart of n_cells_sel / n_cells_inf per level with uncovered levels in red); and the per-cell detail seed, cluster_id, key1_value, key2_value, fold (-1 for D_inf), split (sel / inf).

The HTML report gains a 000_saturation tab, which embeds diag2_rloss_profile.png and summary_alpha_profile_<treatment>.png and shows decision, decision_reason, alpha_hat, half_saturation_T, T_ref, d_bar, se, diag2_ratio, saturation_at_T_ref and r_T.

Saturation-aware optimization

With optimize=True, simple_optimization=False, one treatment of type spend and selected_family != None, Phase 4 solves

max sum_t theta_t fbar(T_t) s.t. sum_t T_t <= K', 0 <= T_t <= T_bar_t^eff

in one KKT pass: one allocation unit is one row (item × geo × date) inside opt_periods ∩ opt_geos × opt_items; theta_t = p(X_t)' beta_hat; the support cap T_bar_t = max(Q_q(cell_t), T_t^actual) is the optimization.cap_quantile-quantile q of the positive actual treatment within the unit's geo × item cell, floored at the unit's own actual spend (the default q = 1.0 is the cell maximum, i.e. the historical cap; with q = 0.95 an outlier day keeps its own actual spend as its cap — it can be cut by the KKT condition or the ROI cap but no unit is ever pushed above its cap — while every other unit is capped at the cell's 95 % point); the ROI threshold r = opt_threshold_roi is a per-unit cap T_t^r (the root of theta_t fbar(T)/T = r; 0 when theta_t fbar'(0) <= r; inf when r = 0), T_bar_t^eff = min(T_bar_t, T_t^r); units with theta_t <= 0 receive 0. T_t*(mu) = clip((fbar')^-1(mu / theta_t), 0, T_bar_t^eff) with the closed forms of all four families and a single bisection on mu so that sum_t T_t*(mu) = K' (mu* = 0 when the caps absorb the budget). The former "exclude below the threshold → re-solve → iterate" loop no longer exists. budget_mode = roi_target skips the budget and sets T_t* = T_bar_t^eff. Any other combination (several treatments or selected_family = None) falls back to the existing optimizer with a warning; a non-spend treatment type forces simple_optimization=True at construction (see above); simple_optimization=False is otherwise never overridden when saturation_config is given, and opt_threshold_roi < 0 raises at initialization.

The saturation optimizer additionally writes optimized_incremental_KPI_<t> and optimized_roi_<t> into the consolidated CSV and into 007_optimization/optimization_metrics.csv, where every geo / item / all row also carries budget_mode, budget_reallocate, reallocate_within, cap_quantile, assumptions (no_carryover), K, K_out, unspent_budget, marginal_roi_at_optimum (mu*; with marginal_roi_at_optimum_min / _max, equal to it under reallocate_within = "all"), n_units, n_units_excluded_negative_roi, n_units_capped_by_roi and n_units_at_support_cap (K, K_out, unspent_budget and the counts are aggregated over the units of that row, so a geo row can show a negative unspent_budget when budget was reallocated into it; budget_mode, the scope and mu* describe the whole problem). The same values appear in the HTML descriptive block of the 007_optimization tab. When treatment_type="spend" it also writes, per geo, per item and for the total, the two-line chart <prefix>_original_vs_optimal_treatment.png (purple original vs green optimal treatment); the incremental effect is covered by the before / after chart below.

Both optimizers write, per geo, per item and for the total, the outcome chart <prefix>_<outcome>.png: the red observed outcome (Original KPI) against the blue projected outcome under the optimized spend (Estimated total KPI with optimal treatment, Y + (KPI_opt - KPI) summed per date — for the simple optimizer the KPI of the cut rows leaves the total). The HTML tabs 007_optimization and roi_output_007_optimization (<treatment>) show it in the same position whatever simple_optimization is: right after the optimizer's own chart and above the before / after charts.

Both optimizers (simple and saturation-aware) write, per geo, per item and for the total, three before / after charts — two line charts stacked vertically on an identical vertical scale, "before optimization" on top and "after optimization" below, every panel on exactly the width and x-axis layout of <prefix>_timeseries_spend_vs_est.png so the dates line up: <prefix>_timeseries_spend_before_after.png (purple actual spend (treatment_spend_cols) / dark green optimized spend), <prefix>_timeseries_incremental_effect_before_after.png (red incremental effect of the actual spend / blue incremental effect under the optimized spend) and <prefix>_timeseries_roi_before_after.png (red ROI per date = incremental KPI / spend, 0 without spend / blue ROI under the optimized spend, with the ROI threshold as a dotted reference line) — every panel of the three charts is exactly half the height of spend_vs_est at the same vertical scale, so each before / after pair fits one screen — and optimization_yearly_summary.csv (per dimension and calendar year plus Total: spend_before, incremental_kpi_before, roi_before, spend_after, incremental_kpi_after, roi_after). The HTML tabs 007_optimization and roi_output_007_optimization (<treatment>) show the three charts in that order under the existing optimization chart and render the CSV as a matrix below them that follows the dimension dropdown: one row per metric (Spend, Incremental KPI, ROI) under "Before optimization" / "After optimization", one column per year plus Total, numbers only in the cells; the matrix always fits the width of the tab (its font size adapts to the widest number, down to 9 px, after which the label columns narrow and, as a last resort, the matrix scrolls inside the tab rather than cutting digits).

HTML report

When generate_html=True the interactive report is written to the root of output_dir as <treatment>_<outcome>_<family>.html, for example tv_spend-digital_spend_sales_power.html. Multiple treatments are joined with a hyphen, and the _<family> segment is omitted when saturation is not enabled (tv_spend_sales.html).

---

Analyzing Results with AI (Browser Chat or ADK Agent)

There are two ways to explore the causal estimates with Gemini:

1. Directly in report.html (no server required): Open <output_dir>/report.html in your browser, select the ✨ AI Assistant (Gemini) tab, and enter a Gemini API Key to chat interactively with the precomputed CausalResultsReader outputs and bundled retail-analytics skill.
2. Local ADK Agent (generate_agent=True): By default, the pipeline also writes a runnable ADK agent to <output_dir>/agents/causal_ml_insight_agent/. The agent answers business questions such as "Which promotion pays back, where, and when?" from the pipeline outputs. Its tools read the results and do the arithmetic deterministically, so the LLM interprets the DoubleML estimates instead of computing causal effects itself. The domain knowledge comes from Agent Skills bundled in skills/, starting with retail-analytics. Coding agents can use these skills directly too.

cp causal_results/agents/causal_ml_insight_agent/.env.example \
   causal_results/agents/causal_ml_insight_agent/.env  # Set your project.
adk web causal_results/agents

To regenerate the agent with your own skills or model, or to add one to the results of an earlier run:

from doubleml_pipeline.agent import generate_agent

generate_agent("./causal_results", skills_dir="./my_skills")

If you run one orchestrator per channel, one agent can cover all of the runs. It compares the channels side by side but never adds their effects up, since each run is estimated separately:

generate_agent(
    ["results/Spend_TV", "results/Spend_Digital", "results/Spend_Flyer"],
    output_dir="results",
)
# adk web results/agents

User Guide

For a comprehensive step-by-step walkthrough of pipeline settings, modeling workflows, and visual slide guides, refer to the User Guide.

---

Engaging with the Project

Documentation & Examples

• User Guide: For a step-by-step walkthrough of pipeline settings and visual slide guides, refer to the User Guide.
• API Reference: For a comprehensive list of modules, classes, and function signatures, refer to the API Reference.
• Output CSV Schema: For detailed specifications of all output CSV files, including column names, data types, and descriptions, refer to the Output CSV Schema.
• Interactive Samples: Check out the sample Jupyter notebooks (IN_skincare and JP_retailer), which walk through loading realistic simulated panel marketing data, engineering lagged features, scaling, running the DML pipeline, and plotting the results.

Contributing & Feedback

We welcome contributions!

• Reporting Issues: If you encounter bugs or want to request features, please open an issue in the GitHub repository issue tracker.
• Submitting Changes: Please see CONTRIBUTING.md for our guidelines on submitting pull requests.

---

Acknowledgments & Dataset Credits

We gratefully acknowledge our partners for their collaboration on the simulation datasets included under examples/:

• examples/IN_skincare/: Simulation data and example materials provided in collaboration with Analytic Edge.
• examples/JP_retailer/: Simulation data designed and provided by Hakuhodo Inc., informed by real-world Japanese marketing dynamics.

---

References

1. Debiased Machine Learning method Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, James Robins. Double/debiased machine learning for treatment and structural parameters. arXiv:1608.00060

2. Hyperparameter Tuning for Causal Inference with Double Machine Learning: A Simulation Study Martin Spindler, et al. arXiv:2402.04674

3. Debiased Machine Learning of Conditional Average Treatment Effects and Other Causal Functions Vira Semenova, Victor Chernozhukov. arXiv:1702.06240

4. Multiway Cluster Robust Double/Debiased Machine Learning Harold D. Chiang, Kengo Kato, Yukitoshi Matsushita, Takuya Ishihara. arXiv:1909.03489