QTE benchmark: py-qte vs R qte

# This document requires two prepared environments and fails loudly if either is
# missing:
#   - Python: `uv sync` (creates ../.venv)
#   - R:      `renv::restore()` inside docs/ (installs the pinned packages)
# Params files live in docs/: fast_params.yml for quick renders and
# prd_params.yml for the full run, e.g.:
#   quarto render docs --execute-params docs/fast_params.yml
#   quarto render docs --execute-params docs/prd_params.yml
# Override single values with -P, e.g. -P n_sim_confounded:1000000
docs_dir <- dirname(knitr::current_input(dir = TRUE))

# Force all libraries to be loaded from the renv environment.
options(renv.config.synchronized.check = FALSE)
renv::load(project = docs_dir, quiet = TRUE)

# Point reticulate at the Python environment one level above docs/.
library(reticulate)
venv <- file.path(docs_dir, "..", ".venv")
use_virtualenv(normalizePath(venv), required = TRUE)

# Report when every chunk tagged with `#| time_it: true` starts and how long it
# took to the console during compilation (not into the rendered document).
knitr::knit_hooks$set(time_it = local({
  start <- NULL
  function(before, options) {
    if (before) {
      start <<- Sys.time()
      message(sprintf("[%s] start %s", options$label, format(start, "%H:%M:%S")))
    } else {
      elapsed <- as.numeric(difftime(Sys.time(), start, units = "secs"))
      message(sprintf("[%s] time taken %.2fs", options$label, elapsed))
      invisible(NULL)
    }
  }
}))

library(qte)

# For parallel bootstrap in qte (r-package) we need BLAS/OpenMP to be single-threaded.
library(RhpcBLASctl)
blas_set_num_threads(1L)
omp_set_num_threads(1L)
pbapply::pboptions(type = "none")

N_ITER <- as.integer(params$n_iter)
N_ITER_SIMPLE <- as.integer(params$n_iter_simple)
N_WORKERS <- as.integer(params$n_workers)
SEED <- as.integer(params$seed)
N_SIM_UNCONFOUNDED <- as.integer(params$n_sim_unconfounded)
N_SIM_CONFOUNDED <- as.integer(params$n_sim_confounded)
N_SIM_PY_QTE_ONLY <- as.integer(params$n_sim_py_qte_only)
N_ITER_SIM_UNCONFOUNDED <- as.integer(params$n_iter_sim_unconfounded)
N_ITER_SIM_CONFOUNDED <- as.integer(params$n_iter_sim_confounded)
N_ITER_SIM_PY_QTE_ONLY <- as.integer(params$n_iter_sim_py_qte_only)

QS <- c(0.25, 0.5, 0.75)
QS_SIMPLE <- c(0.25, 0.5, 0.75, 0.9)
X_R <- ~ age + I(age^2) + education + black + hispanic + married + nodegree
X_PY <- "age + I(age**2) + education + black + hispanic + married + nodegree"

import time

import polars as pl
from great_tables import GT, loc, md, px, style

from qte.bootstrap import BootstrapConfig
from qte.cross_sectional import estimate_aipw_qte, estimate_simple_qte
from qte.datasets import load_lalonde

DATA = load_lalonde(controls_source="psid")
qs = [float(q) for q in r.QS]
qs_simple = [float(q) for q in r.QS_SIMPLE]

def make_source_note(n_obs: int, n_iter: int, n_workers: int) -> str:
    return (
        f"N = {int(n_obs):,}; {int(n_iter):,} bootstrap iterations; "
        f"{int(n_workers)} workers for parallel computations."
    )


def make_timing_table(ds: pl.DataFrame, estimation: str, source_note: str) -> GT:
    return (
        GT(ds.select("package", "computation", "seconds"))
        .cols_label(package="Package", computation="Computation", seconds="Seconds")
        .fmt_number(columns=["seconds"], decimals=1, use_seps=True)
        .tab_header(title=md("**Speed Comparison**"), subtitle=estimation)
        .tab_source_note(source_note)
    )


def make_comparison_data(qs, series) -> pl.DataFrame:
    """Tidy estimate/SE data for a package-vs-package comparison.

    ``series`` maps a ``(package, computation)`` label to an ``(effects, ses)``
    pair, each ordered as ``[QTE(q) for q in qs] + [ATT]``.
    """
    n_qte = len(qs)
    records = []
    for (package, computation), (effects, ses) in series.items():
        for i, q in enumerate(qs):  # quantiles
            records.append(
                {
                    "effect": "QTE",
                    "quantile": f"{q:.2f}",
                    "package": package,
                    "computation": computation,
                    "estimate": float(effects[i]),
                    "se": float(ses[i]),
                }
            )
        records.append(
            {
                "effect": "ATT",
                "quantile": "",
                "package": package,
                "computation": computation,
                "estimate": float(effects[n_qte]),
                "se": float(ses[n_qte]),
            }
        )
    return pl.DataFrame(records)


def make_estimate_comparison_table(data: pl.DataFrame, estimation: str, source_note: str) -> GT:
    """Render tidy comparison data (see ``make_comparison_data``).

    Rows are grouped under ``QTE``/``ATT``; package names become the top-level
    column spanners with the computation nested underneath. Estimate and
    standard-error columns get distinct backgrounds, and packages are separated
    by a vertical rule.
    """
    pairs = data.select("package", "computation").unique(maintain_order=True).rows()
    packages = list(dict.fromkeys(package for package, _ in pairs))
    computations = list(dict.fromkeys(computation for _, computation in pairs))

    def series(package, computation):
        # filter preserves the QTE...ATT order produced by make_comparison_data
        return data.filter((pl.col("package") == package) & (pl.col("computation") == computation))

    base = series(*pairs[0])
    frame = {"effect": base["effect"].to_list(), "quantile": base["quantile"].to_list()}
    leaf_labels = {}
    for package, computation in pairs:
        prefix = f"{package}.{computation}"
        sub = series(package, computation)
        frame[f"{prefix}.est"] = sub["estimate"].to_list()
        frame[f"{prefix}.se"] = sub["se"].to_list()
        leaf_labels[f"{prefix}.est"] = "Estimate"
        leaf_labels[f"{prefix}.se"] = md("Standard<br>Error")

    table = pl.DataFrame(frame)
    est_cols = [c for c in table.columns if c.endswith(".est")]
    se_cols = [c for c in table.columns if c.endswith(".se")]

    gt = (
        GT(table, groupname_col="effect")
        .row_group_order(["QTE", "ATT"])
        .cols_label(quantile="Quantile", **leaf_labels)
        .fmt_number(columns=est_cols + se_cols, decimals=1, use_seps=True)
        .tab_header(title=md("**Results Comparison**"), subtitle=estimation)
        .tab_source_note(source_note)
        .tab_style(style=style.fill(color="#E8F1FA"), locations=loc.body(columns=est_cols))
        .tab_style(style=style.fill(color="#FBF0E1"), locations=loc.body(columns=se_cols))
        .tab_style(style=style.fill(color="#D3E4F5"), locations=loc.column_labels(columns=est_cols))
        .tab_style(style=style.fill(color="#F6E0C4"), locations=loc.column_labels(columns=se_cols))
    )
    for package in packages:
        # narrower spanners first so the package spanner sits on top
        package_columns = []
        for computation in computations:
            est_col, se_col = f"{package}.{computation}.est", f"{package}.{computation}.se"
            if est_col not in frame:
                continue
            gt = gt.tab_spanner(
                computation, columns=[est_col, se_col], id=f"{package}|{computation}"
            )
            package_columns += [est_col, se_col]
        gt = gt.tab_spanner(package, columns=package_columns, id=f"{package}|pkg")
    if len(packages) > 1:
        divider = "#8A8A8A"
        second_package = packages[1]
        first_computation = next(
            computation for package, computation in pairs if package == second_package
        )
        first_of_second = f"{second_package}.{first_computation}.est"
        gt = gt.tab_style(
            style=style.borders(sides="left", color=divider, style="solid", weight=px(2)),
            locations=[
                loc.body(columns=first_of_second),
                loc.column_labels(columns=first_of_second),
                loc.spanner_labels(ids=[f"{second_package}|pkg"]),
                loc.spanner_labels(ids=[f"{second_package}|{first_computation}"]),
            ],
        )
    return gt

Summary

Here, we compare the performance of the py-qte and the qte package. We use a medium-sized machine (see the Environment section for details). The comparisons are based on the relatively small Lalonde dataset that comes with the qte package and on a medium-sized simulated dataset inspired by this section from the documentation of the DoubleML package.

In particular, estimators that involve an outcome model are relatively computationally intensive since they involve fitting a quantile regression model on a grid of quantiles in each bootstrap run. So speed matters here.

Even for the small Lalonde dataset (2,675 observations), we find sizeable differences between the qte and the py-qte packages in sequential computations. These can be closely tied to different implementation approaches of the respective estimators. Parallelized approaches benefit qte relatively more (in this setting). See Table 2. For estimation with confounders, the same holds true, but even after parallelization the gap remains sizeable, and the qte package takes about 3 times as long as the py-qte package (see Table 4).

For a larger simulated dataset (250,000 observations), we find that without confounders py-qte is orders of magnitude as fast even after parallelization. See Table 6. For the AIPW estimator we find that the parallelized py-qte version on a medium-sized dataset (10,000 observations and 5 covariates) with a reasonable number of bootstrap iterations is more than 20 times faster than the qte version. See Table 9.

Lastly, the py-qte version can provide results for 100,000 observations and 200 bootstrap iterations in about 260 seconds. This is about half the time the parallelized qte version took for a dataset that has only 20% of the rows. Interestingly, the parallelized py-qte version with 10 workers is about 10 times faster than the sequential version, which suggests that the chosen parallelization strategy based on threads is close to optimal for larger datasets. See Table 11.

Overall, these results stress the importance of algorithmic choices and parallelization strategies for the bootstrap iterations.

Cross-Sectional Data

We compare cross-sectional estimators on the Lalonde dataset. The Lalonde dataset is a relatively small dataset with only 2,675 observations.

Estimation under Unconfoundedness

# Unconditional QTE with no covariates == R's IPW branch with x = NULL.
run_r_simple <- function(cl) {
  set.seed(SEED)
  t0 <- Sys.time()
  res <- unc_qte(
    yname = "re78", dname = "treat", data = lalonde.psid, xformla = NULL,
    probs = QS_SIMPLE, est_method = "ipw", target = "qtt",
    cl = cl, biters = N_ITER_SIMPLE
  )
  list(
    seconds = as.numeric(Sys.time() - t0, units = "secs"),
    vals = list(
      qte = as.numeric(res$qte), se = as.numeric(res$qte.se),
      ate = as.numeric(res$ate), ate_se = as.numeric(res$ate.se)
    )
  )
}

r_simple_seq <- run_r_simple(1)
r_simple_par <- run_r_simple(N_WORKERS)
r_simple_seq_seconds <- r_simple_seq$seconds
r_simple_par_seconds <- r_simple_par$seconds
r_simple_seq_vals <- r_simple_seq$vals
r_simple_par_vals <- r_simple_par$vals

def run_py_simple(n_workers):
    cfg = BootstrapConfig(n_iter=int(r.N_ITER_SIMPLE), seed=int(r.SEED), n_workers=n_workers)
    t0 = time.perf_counter()
    res = estimate_simple_qte(DATA, "re78", "treat", qs=qs_simple, bootstrap_config=cfg)
    return time.perf_counter() - t0, res


py_simple_seq_seconds, res_py_simple_seq = run_py_simple(1)
py_simple_par_seconds, res_py_simple_par = run_py_simple(int(r.N_WORKERS))

simple_order = [f"QTE({q:.2f})" for q in qs_simple] + ["ATE"]


def py_simple_cols(res):
    rows = {
        f"QTE({row['q']:.2f})": (row["effect"], row["se"]) for row in res.qtt.iter_rows(named=True)
    }
    rows["ATE"] = (res.att["effect"][0], res.att["se"][0])
    return [rows[k][0] for k in simple_order], [rows[k][1] for k in simple_order]


def r_simple_cols(vals):
    eff = [float(x) for x in vals["qte"]] + [float(vals["ate"])]
    err = [float(x) for x in vals["se"]] + [float(vals["ate_se"])]
    return eff, err


simple_source_note = make_source_note(DATA.height, int(r.N_ITER_SIMPLE), int(r.N_WORKERS))
simple_table_html = make_estimate_comparison_table(
    make_comparison_data(
        qs_simple,
        {
            ("py-qte", "Sequential"): py_simple_cols(res_py_simple_seq),
            ("py-qte", "Parallel"): py_simple_cols(res_py_simple_par),
            ("qte", "Sequential"): r_simple_cols(r.r_simple_seq_vals),
            ("qte", "Parallel"): r_simple_cols(r.r_simple_par_vals),
        },
    ),
    "Estimation under Unconfoundedness",
    simple_source_note,
).as_raw_html()

data = pl.DataFrame(
    {
        "language": ["Python", "Python", "R", "R"],
        "package": ["py-qte", "py-qte", "qte", "qte"],
        "computation": ["Sequential", "Parallel", "Sequential", "Parallel"],
        "workers": [1, int(r.N_WORKERS), 1, int(r.N_WORKERS)],
        "iterations": [int(r.N_ITER_SIMPLE)] * 4,
        "seconds": [
            py_simple_seq_seconds,
            py_simple_par_seconds,
            float(r.r_simple_seq_seconds),
            float(r.r_simple_par_seconds),
        ],
    }
)

simple_timing_html = make_timing_table(
    data, "Estimation under Unconfoundedness", simple_source_note
).as_raw_html()

Table 1: Unconditional QTE Estimates on the Lalonde Data.

Results Comparison
Estimation under Unconfoundedness
	py-qte				qte
Quantile	Sequential		Parallel		Sequential		Parallel
	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error
QTE
0.25	−11,041.0	586.1	−11,041.0	586.1	−11,041.0	617.0	−11,041.0	604.8
0.50	−16,455.9	743.7	−16,455.9	743.7	−16,455.9	727.9	−16,455.9	723.1
0.75	−19,911.5	814.7	−19,911.5	814.7	−19,911.5	787.0	−19,911.5	767.2
0.90	−23,867.8	1,948.6	−23,867.8	1,948.6	−23,839.0	1,904.7	−23,839.0	1,916.2
ATT
	−15,204.8	688.6	−15,204.8	688.6	−15,204.8	654.2	−15,204.8	651.2
N = 2,675; 1,000 bootstrap iterations; 10 workers for parallel computations.

Table 2: Unconditional QTE Runtimes on the Lalonde Data.

Speed Comparison
Estimation under Unconfoundedness
Package	Computation	Seconds
py-qte	Sequential	0.7
py-qte	Parallel	0.4
qte	Sequential	3.3
qte	Parallel	0.5
N = 2,675; 1,000 bootstrap iterations; 10 workers for parallel computations.

Notice first that py-qte replicates the point estimates of the qte package and produces similar standard errors. Notice that py-qte produces consistent standard error estimates across the sequential and the parallelized runs, as it should.

We see that py-qte is in a sequential scenario much faster than qte. This is probably mostly because np.quantile provides an efficient implementation of computing weighted quantiles whereas the qte package runs a solver for each quantile. The speed-up from parallelization is very large for qte and sizeable but relatively low for py-qte.

AIPW Estimator

set.seed(SEED)
t0 <- Sys.time()
res_r_seq <- unc_qte(
  yname = "re78", dname = "treat", data = lalonde.psid, xformla = X_R,
  probs = QS, est_method = "aipw", target = "qte", cl = 1, biters = N_ITER
)
r_seq_seconds <- as.numeric(Sys.time() - t0, units = "secs")

r_seq_vals <- list(
  qte = as.numeric(res_r_seq$qte), se = as.numeric(res_r_seq$qte.se),
  ate = as.numeric(res_r_seq$ate), ate_se = as.numeric(res_r_seq$ate.se)
)

set.seed(SEED)
t0 <- Sys.time()
res_r_par <- unc_qte(
  yname = "re78", dname = "treat", data = lalonde.psid, xformla = X_R,
  probs = QS, est_method = "aipw", target = "qte", cl = N_WORKERS, biters = N_ITER
)
r_par_seconds <- as.numeric(Sys.time() - t0, units = "secs")

r_par_vals <- list(
  qte = as.numeric(res_r_par$qte), se = as.numeric(res_r_par$qte.se),
  ate = as.numeric(res_r_par$ate), ate_se = as.numeric(res_r_par$ate.se)
)

cfg = BootstrapConfig(n_iter=int(r.N_ITER), seed=int(r.SEED), n_workers=1)
t0 = time.perf_counter()
res_py_seq = estimate_aipw_qte(
    DATA,
    "re78",
    "treat",
    qs=qs,
    ps_x_formular=r.X_PY,
    or_x_formular=r.X_PY,
    bootstrap_config=cfg,
)
py_seq_seconds = time.perf_counter() - t0

cfg = BootstrapConfig(n_iter=int(r.N_ITER), seed=int(r.SEED), n_workers=int(r.N_WORKERS))
t0 = time.perf_counter()
res_py_par = estimate_aipw_qte(
    DATA,
    "re78",
    "treat",
    qs=qs,
    ps_x_formular=r.X_PY,
    or_x_formular=r.X_PY,
    bootstrap_config=cfg,
)
py_par_seconds = time.perf_counter() - t0

order = [f"QTE({q:.2f})" for q in qs] + ["ATE"]


def py_cols(res):
    rows = {
        f"QTE({row['q']:.2f})": (row["effect"], row["se"]) for row in res.qtt.iter_rows(named=True)
    }
    rows["ATE"] = (res.att["effect"][0], res.att["se"][0])
    return [rows[k][0] for k in order], [rows[k][1] for k in order]


def r_cols(vals):
    eff = [float(x) for x in vals["qte"]] + [float(vals["ate"])]
    err = [float(x) for x in vals["se"]] + [float(vals["ate_se"])]
    return eff, err


py_seq_e, py_seq_s = py_cols(res_py_seq)
py_par_e, py_par_s = py_cols(res_py_par)
r_seq_e, r_seq_s = r_cols(r.r_seq_vals)
r_par_e, r_par_s = r_cols(r.r_par_vals)

aipw_data = make_comparison_data(
    qs,
    {
        ("py-qte", "Sequential"): (py_seq_e, py_seq_s),
        ("py-qte", "Parallel"): (py_par_e, py_par_s),
        ("qte", "Sequential"): (r_seq_e, r_seq_s),
        ("qte", "Parallel"): (r_par_e, r_par_s),
    },
)

aipw_source_note = make_source_note(DATA.height, int(r.N_ITER), int(r.N_WORKERS))
table_html = make_estimate_comparison_table(
    aipw_data, "Estimation via AIPW", aipw_source_note
).as_raw_html()

timing = pl.DataFrame(
    {
        "language": ["Python", "Python", "R", "R"],
        "package": ["py-qte", "py-qte", "qte", "qte"],
        "computation": ["Sequential", "Parallel", "Sequential", "Parallel"],
        "workers": [1, int(r.N_WORKERS), 1, int(r.N_WORKERS)],
        "iterations": [int(r.N_ITER)] * 4,
        "seconds": [py_seq_seconds, py_par_seconds, float(r.r_seq_seconds), float(r.r_par_seconds)],
    }
)

timing_html = make_timing_table(timing, "Estimation via AIPW", aipw_source_note).as_raw_html()

Table 3: AIPW Comparison of py-qte and qte.

Results Comparison
Estimation via AIPW
	py-qte				qte
Quantile	Sequential		Parallel		Sequential		Parallel
	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error
QTE
0.25	−7,647.4	2,089.0	−7,647.4	2,089.0	−7,646.8	1,899.1	−7,646.8	1,919.6
0.50	−12,682.7	2,806.1	−12,682.7	2,806.1	−12,684.5	2,779.0	−12,684.5	2,837.9
0.75	−16,522.7	1,788.6	−16,522.7	1,788.6	−16,522.6	1,883.4	−16,522.6	1,814.2
ATT
	−12,535.5	1,771.9	−12,535.5	1,771.9	−13,194.8	1,582.2	−13,194.8	1,668.3
N = 2,675; 100 bootstrap iterations; 10 workers for parallel computations.

In terms of results, we see in Table 3 that the point estimates are close between py-qte and qte. The standard errors are once more consistent across the sequential and parallelized runs.

Table 4: AIPW Runtimes Comparison of py-qte and qte.

Speed Comparison
Estimation via AIPW
Package	Computation	Seconds
py-qte	Sequential	18.7
py-qte	Parallel	3.4
qte	Sequential	73.2
qte	Parallel	10.3
N = 2,675; 100 bootstrap iterations; 10 workers for parallel computations.

We further see in Table 4 that the sequential version of py-qte is much faster than the qte version. This is probably largely due to the fact that py-qte computes \(F_{AIPW}(y)\) on the y-grid in a single (vectorized) pass whereas the qte version computes each point on the y-grid separately in a loop. The parallelized version of py-qte retains a sizeable advantage over the qte version although the gap is considerably narrower than in the sequential case.

Simulated Data

The sections above use the Lalonde dataset, which has only a few thousand observations. To see how well the py-qte and the qte package do on large dataset we resort to simulated data.

from qte.cross_sectional import simulate_covariate_data, simulate_simple_data


def add_correct_features(ds):
    return ds.with_columns(
        ind1=(pl.col("x1") > 0.1).cast(pl.Float64),
        ind2=((pl.col("x0") * pl.col("x2")) > 0).cast(pl.Float64),
        sq=(2.5 + 0.3 * pl.col("x1")).sqrt(),
    )


EFFECT_SCALE = 1000.0
sim_simple = simulate_simple_data(
    int(r.N_SIM_UNCONFOUNDED),
    treatment_effect=EFFECT_SCALE,
    error_var=EFFECT_SCALE**2,  # keep the residual sd on the same scale as the effect
    seed=int(r.SEED),
)
sim_cov = add_correct_features(
    simulate_covariate_data(
        int(r.N_SIM_CONFOUNDED), effect_scale=EFFECT_SCALE, seed=int(r.SEED) + 1
    )
)

# Plain dicts so reticulate hands R a data.frame (polars objects are not converted).
SIM_SIMPLE_R = sim_simple.to_dict(as_series=False)
SIM_COV_R = sim_cov.to_dict(as_series=False)

sim_simple_df <- as.data.frame(py$SIM_SIMPLE_R)
sim_cov_df <- as.data.frame(py$SIM_COV_R)

run_r_sim <- function(data, xformla, est_method, target, cl, biters) {
  set.seed(SEED)
  t0 <- Sys.time()
  res <- unc_qte(
    yname = "y", dname = "treat", data = data, xformla = xformla,
    probs = QS, est_method = est_method, target = target,
    cl = cl, biters = biters
  )
  list(
    seconds = as.numeric(Sys.time() - t0, units = "secs"),
    vals = list(
      qte = as.numeric(res$qte), se = as.numeric(res$qte.se),
      ate = as.numeric(res$ate), ate_se = as.numeric(res$ate.se)
    )
  )
}

r_sim_simple_par <- run_r_sim(sim_simple_df, NULL, "ipw", "qtt", N_WORKERS, N_ITER_SIM_UNCONFOUNDED)
r_sim_cov_par <- run_r_sim(sim_cov_df, ~ x3 + x4 + ind1 + ind2 + sq, "aipw", "qte", N_WORKERS, N_ITER_SIM_CONFOUNDED)

def run_py_sim_simple(n_workers):
    cfg = BootstrapConfig(
        n_iter=int(r.N_ITER_SIM_UNCONFOUNDED), seed=int(r.SEED), n_workers=n_workers
    )
    t0 = time.perf_counter()
    res = estimate_simple_qte(sim_simple, "y", "treat", qs=qs, bootstrap_config=cfg)
    return time.perf_counter() - t0, res


def run_py_sim_cov(n_workers):
    cfg = BootstrapConfig(
        n_iter=int(r.N_ITER_SIM_CONFOUNDED), seed=int(r.SEED), n_workers=n_workers
    )
    t0 = time.perf_counter()
    res = estimate_aipw_qte(
        sim_cov,
        "y",
        "treat",
        qs=qs,
        ps_x_formular="x1 + x3 + I(x0 > 0)",
        or_x_formular="x3 + x4 + ind1 + ind2 + sq",
        bootstrap_config=cfg,
    )
    return time.perf_counter() - t0, res


py_sim_simple_seq_seconds, res_py_sim_simple_seq = run_py_sim_simple(1)
py_sim_simple_par_seconds, res_py_sim_simple_par = run_py_sim_simple(int(r.N_WORKERS))
py_sim_cov_seq_seconds, res_py_sim_cov_seq = run_py_sim_cov(1)
py_sim_cov_par_seconds, res_py_sim_cov_par = run_py_sim_cov(int(r.N_WORKERS))

def make_sim_tables(title, py_seq, py_par, r_par, source_note):
    py_seq_seconds, py_seq_res = py_seq
    py_par_seconds, py_par_res = py_par
    series = {
        ("py-qte", "Sequential"): py_cols(py_seq_res),
        ("py-qte", "Parallel"): py_cols(py_par_res),
    }
    rows = {
        "language": ["Python", "Python"],
        "package": ["py-qte", "py-qte"],
        "computation": ["Sequential", "Parallel"],
        "workers": [1, int(r.N_WORKERS)],
        "seconds": [py_seq_seconds, py_par_seconds],
    }
    if r_par is not None:
        series[("qte", "Parallel")] = r_cols(r_par["vals"])
        rows["language"].append("R")
        rows["package"].append("qte")
        rows["computation"].append("Parallel")
        rows["workers"].append(int(r.N_WORKERS))
        rows["seconds"].append(float(r_par["seconds"]))
    table_html = make_estimate_comparison_table(
        make_comparison_data(qs, series), title, source_note
    ).as_raw_html()
    timing_html = make_timing_table(pl.DataFrame(rows), title, source_note).as_raw_html()
    return table_html, timing_html


sim_simple_table_html, sim_simple_timing_html = make_sim_tables(
    "Simulated data without covariates",
    (py_sim_simple_seq_seconds, res_py_sim_simple_seq),
    (py_sim_simple_par_seconds, res_py_sim_simple_par),
    r.r_sim_simple_par,
    make_source_note(r.N_SIM_UNCONFOUNDED, r.N_ITER_SIM_UNCONFOUNDED, r.N_WORKERS),
)
sim_cov_table_html, sim_cov_timing_html = make_sim_tables(
    "Simulated data with covariates",
    (py_sim_cov_seq_seconds, res_py_sim_cov_seq),
    (py_sim_cov_par_seconds, res_py_sim_cov_par),
    r.r_sim_cov_par,
    make_source_note(r.N_SIM_CONFOUNDED, r.N_ITER_SIM_CONFOUNDED, r.N_WORKERS),
)

Estimation under Unconfoundedness

Table 5: QTE Estimates on the Simulated Unconfounded Data.

Results Comparison
Simulated data without covariates
	py-qte				qte
Quantile	Sequential		Parallel
	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error
QTE
0.25	997.1	5.7	997.1	5.7	997.0	6.3
0.50	1,001.1	4.8	1,001.1	4.8	1,001.1	5.6
0.75	1,000.2	5.4	1,000.2	5.4	1,000.3	5.7
ATT
	997.5	4.1	997.5	4.1	997.5	4.3
N = 250,000; 200 bootstrap iterations; 10 workers for parallel computations.

Table 6: QTE Runtimes on the Simulated Unconfounded Data.

Speed Comparison
Simulated data without covariates
Package	Computation	Seconds
py-qte	Sequential	0.8
py-qte	Parallel	0.1
qte	Parallel	10.2
N = 250,000; 200 bootstrap iterations; 10 workers for parallel computations.

As in the Lalonde comparison, the two packages agree on the point estimates and standard errors. We omit the sequential qte run here because it is very slow on larger data; even against the parallel qte run py-qte is faster. Because the simulated design is randomized, the treatment effect is identified by the unadjusted comparison; simulate_simple_data generates it with the treatment effect set by EFFECT_SCALE.

AIPW Estimator

The simulated dataset has five covariates: the outcome model uses five transformed features of them, whereas the propensity score uses three of the raw covariates. Note that qte::unc_qte takes a single xformla used for both the propensity score and the outcome regression, so the R run below fits both models on the same five features while the py-qte run uses separate formulas. The population (“true”) effects of this design are reported in Table 7. They are obtained by generating a large dataset at the true parameter values and subsequently computing the quantities of interest from that dataset.

import gc

import numpy as np

N_TRUE = int(10e6)
true_ds = simulate_covariate_data(N_TRUE, effect_scale=EFFECT_SCALE, seed=int(r.SEED))
true_y0 = true_ds["y_0"].to_numpy()
true_y1 = true_ds["y_1"].to_numpy()
true_ate = float(np.mean(true_y1 - true_y0))
true_qte = np.quantile(true_y1, qs) - np.quantile(true_y0, qs)

# Drop the large simulated population before the later runs.
del true_ds, true_y0, true_y1
_ = gc.collect()

true_rows = pl.concat(
    [
        pl.DataFrame({"Quantity": ["ATE"], "Estimate": [true_ate]}),
        pl.DataFrame({"Quantity": [f"QTE({q:.2f})" for q in qs], "Estimate": true_qte}),
    ]
)
true_table_html = (
    GT(true_rows)
    .fmt_number(columns=["Estimate"], decimals=1, use_seps=True)
    .tab_header(
        title=md("**True Effects**"),
        subtitle="Confounded simulation (population QTE)",
    )
    .tab_source_note(f"True values from a simulation with N = {N_TRUE:,} units.")
    .as_raw_html()
)

Table 7: True Mean and Quantile Treatment Effects of the Confounded Simulation.

True Effects
Confounded simulation (population QTE)
Quantity	Estimate
ATE	499.5
QTE(0.25)	−498.7
QTE(0.50)	498.7
QTE(0.75)	1,497.6
True values from a simulation with N = 10,000,000 units.

Table 8: AIPW QTE Estimates on the Simulated Confounded Data.

Results Comparison
Simulated data with covariates
	py-qte				qte
Quantile	Sequential		Parallel
	Estimate	Standard Error	Estimate	Standard Error	Estimate	Standard Error
QTE
0.25	−435.1	15.6	−435.1	15.6	−431.9	15.2
0.50	518.2	19.5	518.2	19.5	518.0	19.3
0.75	1,492.7	17.9	1,492.7	17.9	1,494.1	16.7
ATT
	519.6	11.6	519.6	11.6	519.2	10.7
N = 10,000; 200 bootstrap iterations; 10 workers for parallel computations.

Table 9: AIPW QTE Runtimes on the Simulated Confounded Data.

Speed Comparison
Simulated data with covariates
Package	Computation	Seconds
py-qte	Sequential	173.9
py-qte	Parallel	21.6
qte	Parallel	493.9
N = 10,000; 200 bootstrap iterations; 10 workers for parallel computations.

if int(r.N_SIM_PY_QTE_ONLY) > 0:
    sim_py_only = add_correct_features(
        simulate_covariate_data(
            int(r.N_SIM_PY_QTE_ONLY), effect_scale=EFFECT_SCALE, seed=int(r.SEED) + 2
        )
    )
    sim_cfg = BootstrapConfig(
        n_iter=int(r.N_ITER_SIM_PY_QTE_ONLY), seed=int(r.SEED), n_workers=int(r.N_WORKERS)
    )
    t0 = time.perf_counter()
    res_py_only_seq = estimate_aipw_qte(
        sim_py_only,
        "y",
        "treat",
        qs=qs,
        ps_x_formular="x1 + x3 + I(x0 > 0)",
        or_x_formular="x3 + x4 + ind1 + ind2 + sq",
        bootstrap_config=sim_cfg._replace(n_workers=1),
    )
    py_only_seq_seconds = time.perf_counter() - t0
    t0 = time.perf_counter()
    res_py_only_par = estimate_aipw_qte(
        sim_py_only,
        "y",
        "treat",
        qs=qs,
        ps_x_formular="x1 + x3 + I(x0 > 0)",
        or_x_formular="x3 + x4 + ind1 + ind2 + sq",
        bootstrap_config=sim_cfg,
    )
    py_only_par_seconds = time.perf_counter() - t0
    py_only_table_html, py_only_timing_html = make_sim_tables(
        "py-qte only (no R)",
        (py_only_seq_seconds, res_py_only_seq),
        (py_only_par_seconds, res_py_only_par),
        None,
        make_source_note(r.N_SIM_PY_QTE_ONLY, r.N_ITER_SIM_PY_QTE_ONLY, r.N_WORKERS),
    )
else:
    py_only_table_html = ""
    py_only_timing_html = ""

py-qte only

With n_sim_py_qte_only set we run a separate confounded simulation of that size with only py-qte. The tables therefore contain only py-qte columns.

Table 10: AIPW QTE Estimates for the py-qte-only Run.

Results Comparison
py-qte only (no R)
	py-qte
Quantile	Sequential		Parallel
	Estimate	Standard Error	Estimate	Standard Error
QTE
0.25	−494.1	5.1	−494.1	5.1
0.50	500.9	6.1	500.9	6.1
0.75	1,494.7	5.9	1,494.7	5.9
ATT
	499.8	3.7	499.8	3.7
N = 100,000; 200 bootstrap iterations; 10 workers for parallel computations.

Table 11: AIPW QTE Runtimes for the py-qte-only Run.

Speed Comparison
py-qte only (no R)
Package	Computation	Seconds
py-qte	Sequential	2,101.2
py-qte	Parallel	265.1
N = 100,000; 200 bootstrap iterations; 10 workers for parallel computations.

Environment

system_ <- Sys.info()[["sysname"]]


get_linux_os <- function(key) {
  if (!file.exists("/etc/os-release")) {
    return(NA_character_)
  }
  x <- grep(paste0("^", key, "="), readLines("/etc/os-release"), value = TRUE)[1]
  if (is.na(x)) {
    return(NA_character_)
  }
  value <- gsub('^"|"$', "", sub(paste0("^", "PRETTY_NAME", "="), "", x))
  if (grepl("Void Linux", value)) value <- paste0(value, " btw.")
  value
}
proc_value <- function(path, key) {
  if (!file.exists(path)) {
    return(NA_character_)
  }
  x <- grep(key, readLines(path), value = TRUE)[1]
  if (is.na(x)) {
    return(NA_character_)
  }
  trimws(sub("^[^:]*:", "", x))
}

if (system_ == "Darwin") {
  cpu_name <- system("sysctl -n machdep.cpu.brand_string", intern = TRUE)
} else {
  cpu_name <- proc_value("/proc/cpuinfo", "model name")
}

library(ps)
env_stats <- data.frame(
  item = c(
    "OS", "Kernel", "CPU", "Cores (physical)", "Cores (logical)",
    "Memory", "R", "BLAS", "Python"
  ),
  value = c(
    if (system_ == "Linux") get_linux_os("PRETTY_NAME") else system_,
    paste(Sys.info()[["sysname"]], Sys.info()[["release"]], Sys.info()[["machine"]]),
    cpu_name,
    parallel::detectCores(logical = FALSE),
    parallel::detectCores(logical = TRUE),
    sprintf("%.1f GiB", ps_system_memory()$total / 1e9),
    R.version.string,
    unname(extSoftVersion()["BLAS"]),
    sub(" .*", "", reticulate::py_config()$version_string)
  )
)
knitr::kable(env_stats, col.names = c("", ""), row.names = FALSE)

OS	Darwin
Kernel	Darwin 27.0.0 arm64
CPU	Apple M4 Pro
Cores (physical)	14
Cores (logical)	14
Memory	51.5 GiB
R	R version 4.6.1 (2026-06-24)
BLAS	/opt/homebrew/Cellar/openblas/0.3.34/lib/libopenblasp-r0.3.34.dylib
Python	3.13.11