The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
81,867 characters
Calibration Strategies for Robust Causal Estimation: Theoretical and Empirical Insights on Propensity Score-Based Estimators
\if11
{ \singlespacing
\title{\bf \hspace{.2cm}
Calibration Strategies for Robust Causal Estimation: Theoretical and Empirical Insights on Propensity Score-Based Estimators
\\}
\author[a,b]{Sven Klaassen}
\author[c,a,$\dagger$]{Jan Rabenseifner}
\author[d]{Jannis Kueck}
\author[b]{Philipp Bach}
\affil[a]{\it \small Economic AI, Nürnberger Str. 262 A, 93059 Regensburg, Bayern, Germany}
\affil[b]{\it \small Chair of Statistics with Application in Business Analytics,
University of Hamburg Business School, Moorweidenstr. 18, 20148 Hamburg, Hamburg, Germany}
\affil[c]{\it \small Chair of Mathematics and Statistics in Economics, University of Hamburg Business School, Von-Melle-Park 5, 20146 Hamburg, Hamburg, Germany}
\affil[d]{\it \small Düsseldorf Institute for Competition Economics, Heinrich Heine University Düsseldorf, Universitätstr. 1, 20225 Düsseldorf, Nordrhein-Westphalen, Germany
} \date{}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Calibration Approaches for Causal Estimation: Theoretical and Empirical Insights on Propensity Score-Based Estimators}
\end{center}
} \fi
\singlespacing{
\begin{abstract}
\footnotesize
The partitioning of data for estimation and calibration critically impacts the performance of propensity score based estimators like inverse probability weighting (IPW) and double/debiased machine learning (DML) frameworks. We extend recent advances in calibration techniques for propensity score estimation, improving the robustness of propensity scores in challenging settings such as limited overlap, small sample sizes, or unbalanced data. Our contributions are twofold: First, we provide a theoretical analysis of the properties of calibrated estimators in the context of DML. To this end, we refine existing calibration frameworks for propensity score models, with a particular emphasis on the role of sample-splitting schemes in ensuring valid causal inference. Second, through extensive simulations, we show that calibration reduces variance of inverse-based propensity score estimators while also mitigating bias in IPW, even in small-sample regimes. Notably, calibration improves stability for flexible learners (e.g., gradient boosting) while preserving the doubly robust properties of DML. A key insight is that, even when methods perform well without calibration, incorporating a calibration step does not degrade performance, provided that an appropriate sample-splitting approach is chosen.
\end{abstract}
\bigskip
\footnotesize
{\it Key words:} Balancing, Calibration, Double Machine Learning, IPW, Sample Splitting
\bigskip
\baselineskip 24pt
}
\begingroup
\footnotetext{$\dagger$ Corresponding author email: [email removed]}
\addtocounter{footnote}{-0}
\endgroup
\doublespacing
\setlength{\abovedisplayskip}{3pt}
\setlength{\belowdisplayskip}{3pt}
\setlength{\abovedisplayshortskip}{3pt}
\setlength{\belowdisplayshortskip}{3pt}
\section{Introduction}\label{sec:intro}
\subsection{Motivation}
In many settings of the causal inference literature,
researchers are interested in the effect of a (binary) treatment $D\in\{0,1\}$ on an outcome $Y\in \mathbb{R}$. If the treatment is not assigned randomly, a common assumption is the so-called unconfoundedness assumption $ Y(0), Y(1) \perp D \mid X$ assuming that the potential outcomes $Y(1)$ and $Y(0)$ are independent of the actual treatment status $D$ conditional on control variables $X$. Let
\begin{align}\label{def::prop_score}
m_0(x) := P(D=1|X=x) = \mathbb{E}[D|X=x].
\end{align}
be the propensity score. As famously shown in \citet{rosenbaum1983central} conditioning on the propensity score is sufficient to effectively account for the confounding through $X$
\begin{equation*}
Y(0), Y(1) \perp D \mid m_0(X).
\end{equation*}
Consequently, propensity scores are a cornerstone of modern causal inference for addressing confounding in observational studies. They enable balancing treatment and control groups, allowing for unbiased estimation of treatment effects under unconfoundedness. Commonly, propensity scores are used in methods such as inverse probability weighting (IPW), matching, stratification, Bayesian causal inference and more recently double machine learning (DML). Moreover, effective propensity adjustment requires sufficient overlap for propensity score-based estimators. When overlap between treatment and control groups is limited, or treatment assignment is unbalanced, propensity scores can become extreme (i.e. close to 0 or 1) leading to instability of the causal estimates. In such cases, estimators can suffer from inflated variance as extreme weights in IPW disproportionately amplify small errors in propensity score estimation. Similarly, matching algorithms may struggle to find suitable matches, resulting in biased estimates. These challenges highlight the importance of robust and well-calibrated propensity score models to maintain the reliability of causal estimates.
To mitigate instability, researchers often enforce common support by trimming or bounding propensity scores. For example, observations with propensity scores below a certain threshold or outside a predefined range are excluded from analysis. While this approach can reduce variance, it does so at the expense of bias, as it discards valuable information and reduces sample size. Such trade-offs are particularly problematic in small-sample settings, where the exclusion of even a few observations can significantly impact the precision and validity of treatment effect estimates.
To use the propensity score for balancing, the property
\begin{align}\label{eq::calibration_prop_score}
m_0(X) = \mathbb{E}[D|m_0(X)]
\end{align}
is crucial, e.g. see Theorem 2 in \citet{rosenbaum1983central}. In the classification literature, Equation (\ref{eq::calibration_prop_score}) is known as the so-called calibration property.\footnote{To be precise, this notion of calibration is also referred to as conditional calibration, which is equivalent to probabilistic calibration for binary outcomes \citep{Gneiting2011probcalib}.} Intuitively, a (binary) classifier $\hat{m}(\cdot)$ is well calibrated if the percentage of positive labels ($D=1$) is approximately $m$ for all instances with $\hat{m}(X)\approx m$. As the true propensity score is the conditional expectation it is calibrated
\begin{equation*}
m_0(X) = \mathbb{E}[D|X] = \mathbb{E}[\mathbb{E}[D|X]|\mathbb{E}[D|X]] = \mathbb{E}[D|m_0(X)] .\quad P\text{-a.s.}
\end{equation*}
In most settings the true propensity score is not known such that it is typically estimated via some classification algorithm such as logistic regression, random forest, boosting methods or even deep neural networks. Consequently, the resulting estimator $\hat{m}(\cdot)$ of $m_0(\cdot)$ might not be calibrated, e.g. the percentage of treated units with $\hat{m}(X)\approx m$ might differ substantially from $m$ for certain values of $m\in (0,1)$. Accurately estimating treatment probabilities is crucial for valid causal inference. Inverse propensity score estimators aim to balance covariates between treatment and control groups. This balancing ensures more reliable and unbiased treatment effect estimates by aligning the covariate distributions, improving the comparability of the treatment groups in terms of observable characteristics. Whereas methods like isotonic regression and Platt scaling refer to the calibration of predictions, critical questions remain about the optimal integration of these calibration techniques in causal estimation workflows.
\subsection{Related Literature}
Recent advances in calibration for propensity score estimation and causal inference emphasize three interconnected themes: the adaptation of machine learning calibration techniques to causal settings, stabilization strategies for inverse probability weighting (IPW), and theoretical insights into finite-sample performance. Initial studies by \citet{deshpande2023calibrated} showed significant improvements using single split calibration on data with deterministic treatment assignments and complex settings with hidden confounders. They specifically highlight variance reduction properties and analyze the regret of recalibration of propensity scores. Their proofs are independent of the sample size \(N\), as their calibration framework consists of splitting the data into training and calibration sets. This distinction is crucial for our work, as we focus on sample size dependent calibration algorithms.
\citet{Gutman2024platt} demonstrated that post-processing propensity scores using methods such as Platt scaling, referred to as post-calibration, improves treatment effect estimation by correcting propensity score distortions beyond covariate balancing. Calibration significantly benefits miscalibrated models like tree-based methods (e.g. gradient boosting), while logistic regression, contrary to the findings of \citet{deshpande2023calibrated}, shows minimal gains. This suggests that post-calibration allows flexible learners to retain their predictive accuracy while enhancing robustness to data challenges such as small sample sizes, model misspecification, class imbalance, or limited overlap, thus challenging the conventional trade-off between model complexity and calibration. \citet{ballinari2025improving} use nested cross-fitting, reserving distinct data folds for propensity estimation and calibration, which is a straightforward application of the double machine learning theory of \citet{Chernozhukov2018dml}. They reported instability in small samples, a limitation attributed to reduced effective sample sizes for both steps. These studies collectively identify a tension between calibration flexibility and stability, particularly in finite-sample regimes. Efforts to address these challenges have taken different paths. For instance, \citet{vanderlaan2024stabilized} proposed a stratified calibration for treated and control groups to stabilize IPW weights. In another work, they extended calibration frameworks to estimate heterogeneous treatment effects \citep{van2023causal}. Meanwhile, \citet{vanderlaan2024automaticdoublyrobustinference} extend automatic debiased machine learning (autoDML) \citep{chernozhukov2022autoDML, chernozhukov2024automaticdebiasedmachinelearning} by calibrating both outcome regression and Riesz representer. They demonstrate that the calibration step yields doubly robust asymptotically linear estimators, reducing nuisance rate requirements.
Theoretical work by \citet{gamarnik1998efficient}, \citet{mammen2007additive} and \citet{wuthrich2023isotonic} established consistency guarantees for isotonic regression, but practical implementations often struggle with convergence rates in small samples, as shown by \citet{yang2019contraction}. Furthermore, \citet{yang2019contraction} established that isotonic projection is non-contractive under the $\ell_{\infty}$ norm, exacerbating edge instability in propensity score estimates. This theoretical insight justifies the empirical requirement of clipping extreme probabilities, which explains the unstable treatment effect estimates in \citet{ballinari2025improving} using isotonic regression under limited overlap.
This paper systematically evaluates how calibration performance depends on data partitioning for propensity estimation and calibration. While existing studies fix specific splitting strategies (e.g., single-split or nested cross-fitting), we show that the choice of partitioning interacts critically with sample size, clipping thresholds, and complexity of the data generating process. For instance, \citet{vanderlaan2024stabilized} suggest calibration on the full-sample, which avoids reserving data exclusively for calibration. This approach can mitigate instability without sacrificing theoretical guarantees, a hypothesis we test across multiple data-generating processes (DGP). Similarly, we reconcile the debates about stratified versus pooled calibration by demonstrating that efficient reuse of cross-fitted propensity scores obviates the need for group-specific adjustments in many settings. Further, beyond the work of \citet{ballinari2025improving}, we provide a theoretical extension of the double machine learning theory to allow for different sample-splitting schemes. By synthesizing these insights, our work clarifies when and how calibration improves ATE estimation, providing a bridge between theoretical calibration properties and practical implementation challenges.
\textbf{Plan of the Paper.}
The rest of the paper is organized as follows. Section \ref{sec:prop_score} introduces propensity score calibration, highlighting different approaches and providing details on the properties of isotonic regression. Section \ref{sec:dml} proposes calibration algorithms for estimating the average treatment effect and establishes theoretical guarantees under double machine learning (DML), including convergence rates and asymptotic normality. Section \ref{sec:sim} demonstrates robustness through simulations across diverse and challenging data-generating processes. Details on the calibration of partially linear regression models, along with proofs of the theorems, are provided in Section 1 of the Supplementary Material. Implementation details for reproducibility, an analysis of how calibration affects the normalization of propensity score weights and covariate balance, as well as extended simulation results and sensitivity analyses, can be found in Sections 2 and 3 of the Supplement.
\section{Propensity Score Calibration}\label{sec:prop_score}
Let $D\in\{0,1\}$ be a binary treatment variable with covariates $X\in\mathcal{X}\subseteq\mathbb{R}^d$. The propensity score is defined as $m_0(x) := P(D = 1 \mid X = x) = \mathbb{E}[D \mid X = x].$ Since the propensity score represents a conditional expectation, it is calibrated such that $m_0(X) = \mathbb{E}[D \mid m_0(X)].$ The goal is to achieve a similar balancing property of an estimated version of $m_0(X)$. Given an estimate $\hat{m}(X)$ of $m_0(X)$, we consider popular calibration methods such as isotonic regression. Generally, we consider calibration procedures based on the pseudo-sample
$((D_1, \hat{m}(X_1)), \dots,$ $(D_N, \hat{m}(X_N)))$. The calibration algorithm approximates $\mathbb{E}[D|\hat{m}(X)]$ which typically differs from $\mathbb{E}[D|m_0(X)]$. In the following, we denote the calibrated propensity score by $\tilde{m}: \mathcal{X} \to [0,1], \quad x\mapsto \tilde{m}(x).$
\subsection{Rate Comparison $\hat{m}(\cdot)$ and $\tilde{m}(\cdot)$}
For any estimate $\hat{m}(\cdot)$ of $m_0(\cdot)$ the mean-squared-error decomposes as
\begin{align}\label{eq:rmse_decomposition}
\|\hat{m}(X) - m_0(X)\|_{P,2} &= \left(\mathbb{E}\left[\textrm{Var}\big(m_0(X)|\hat{m}(X)\big)\right] + \|\mathbb{E}[m_0(X)|\hat{m}(X)] - \hat{m}(X)\|_{P,2}^2 \right)^{1/2}.
\end{align}
The first term denotes the expected precision of $\hat{m}(\cdot)$, while the second term is the calibration error. Calibration procedures minimizing mean square error in the pseudo-sample approximate $\mathbb{E}[D|\hat{m}(X)]$, reducing this error.
\begin{assumption}\label{assumption::calibration_rate}
Let $\tilde{m}(X)$ be an estimator of $\tilde{m}_0(X):=\mathbb{E}[D|\hat{m}(X)]$, with $\|\tilde{m}(X) - \tilde{m}_0(X)\|_{P,2}\le \tilde{\varepsilon}_N.$
\end{assumption}
\begin{lemma}\label{lemma::rmse_decomposition}
Under Assumption \ref{assumption::calibration_rate}:
\begin{align} \label{eq:calibration_rate}
\|\tilde{m}(X) - m_0(X)\|_{P,2} &\le \left(\mathbb{E}\left[\textrm{Var}\big(m_0(X)|\hat{m}(X)\big)\right]\right)^{1/2} + \tilde{\varepsilon}_N.
\end{align}
\end{lemma}
Comparing \eqref{eq:rmse_decomposition} and \eqref{eq:calibration_rate} shows rate differences when the calibration error $\|\mathbb{E}[m_0(X)|\hat{m}(X)] - \hat{m}(X)\|_{P,2}^2$ dominates. Improvements occur if $\tilde{\varepsilon}_N = o(\|\mathbb{E}[m_0(X)|\hat{m}(X)] - \hat{m}(X)\|_{P,2})$. This is relevant in double ML settings requiring $\|\hat{m}(X) - m_0(X)\|_{P,2} = o(N^{-1/4})$: If $\tilde{\varepsilon}_N = o(N^{-1/4})$, $\tilde{m}(\cdot)$ satisfies the requirement with potential rate improvements from faster convergence.
\subsection{Calibration Methods}\label{sec:calib_methods}
\label{ch:isotonic_theory}
Model calibration has roots in meteorological forecasting to address reliability in probabilistic weather predictions \citep{toth2006weathercalib, dawid2014probforecast, gneiting2014weathercalib}.
Calibration is often achieved via post-hoc point calibrators that transform \(\hat{m}(X)\) into $f(\hat{m}(X))$, balancing two aims: (1) calibration validity, and (2) preservation of predictive sharpness (i.e., $f\circ\hat{m}(X)$ approximates \(\hat{m}\)’s discriminative power) \citep{gneiting2007calibrationsharpness, gupta2020binarycalibration}.
Common approaches include parametric methods like Platt scaling, which fits a logistic sigmoid to \(f(X)\) \citep{platt1999plattscale, cox1958sigmoid}; non-parametric methods such as histogram binning, partitioning predictions into fixed intervals \citep{zadrozny2001calibiontrees, gupta2021binaryhistogram}, and isotonic regression, learning a monotonic transform via empirical risk minimization \citep{zadrozny2002multiclass, barlow1972isotonic}; as well as conformal methods like Venn-Abers predictors, refining calibration through cross-conformal inference \citep{vovk2014vennaberspredictors}.
Isotonic calibration, while distribution-free and tuning parameter-free, achieves asymptotic guarantees with \(O_P(N^{-1/3})\) convergence \citep{zhang2002isotonicbounds, van2023causal}. In contrast, histogram binning requires explicit bin specification and trades flexibility for finite-sample validity \citep{gupta2021binaryhistogram}. We focus on three common calibration methods, noting that calibration is an active area of research with potential for improvements and new proposals that require further simulation testing\footnote{The description of Platt scaling and Venn-Abers calibration can be found in Supplement \ref{S-appendix:calib_method}. As our theoretical results are built around isotonic regression, we introduce it here in more detail.}.
\subsubsection*{Isotonic Regression}
\begin{definition}
Given an estimate $\hat{m}(\cdot)$ of the propensity score, perform an isotonic regression as$f = \arg\min_{f\in\mathcal{F}_{\text{iso}}}\sum_{i=1}^N (D_i - (f\circ\hat{m})(X_i))^2$ with $\mathcal{F}_{\text{iso}}$ being the set of non-decreasing functions. The calibrated propensity score is then given by $\tilde{m} = f\circ\hat{m}.$
\end{definition}
Especially, the in-sample calibration property $\mathbb{E}_n[D|\tilde{m}(X_i)] = \tilde{m}(X_i)$
for $i\in\{1,\dots,N\},$ seems to be desirable (e.g. \citet{wuthrich2023isotonic}). We consider an estimated propensity score $\hat{m}(\cdot)$ based on a separate sample such that $\hat{m}(\cdot)$ can be considered a fixed function. Let $U:=\hat{m}(X)$ and define the pseudo sample as $Z:=(D,U)$, where $(Z_i)_{i=1}^N$ are iid. copies of $Z$. Consider the following assumption
\begin{assumption} \label{assumption::monotonicity}
The regression function $\tilde{m}_0(u):= \mathbb{E}[D|U = u]$ is monotone.
\end{assumption}
\begin{remark}
Assumption \ref{assumption::monotonicity} is substantially weaker than requiring monotonicity in the original covariates $X$, intuitively only requiring the preliminary propensity score $\hat{m}(\cdot)$
to be monotone "on average". When $\hat{m}(\cdot)$ is not injective, distinct values of $X$ may map to the same $U$, allowing $\tilde{m}_0(U)$ to average over these $X$ values. This aggregation can smooth out non-monotonic behavior in $X$, making monotonicity in $U$ easier to satisfy. However, if $\hat{m}(\cdot)$ is bijective, each $U$ uniquely determines $X$, and monotonicity of $\tilde{m}_0(U)$ becomes equivalent to monotonicity of $\mathbb{E}[D|X = x]$ using the reparameterization $u = \hat{m}(x)$.
\end{remark}
Define $\tilde{m}$ as the estimator obtained by isotonic regression on the sample $(Z_i)_{i=1}^N$.
\begin{lemma}[Convergence of Isotonic Regression]
\label{lemma::convergence_isotonic_regression}
Under Assumption \ref{assumption::monotonicity}, the isotonic regression estimator $\tilde{m}$ satisfies: $\|\tilde{m}(U) - \tilde{m}_0(U)\|_{P,2} = O_P(N^{-1/3})$.
\end{lemma}
The rate follows from the bracketing entropy bound in Theorem 2.7.5 of \citet{vaart2023empirical} for monotonic functions. Relying on Lemma 3.4.3 of \citet{vaart2023empirical} gives the desired rate as mentioned in the corresponding Section 3.4.3.2. Earlier work by \citet{BirmanSolomjak1967} and \cite{geer2000empirical} established foundational approximation and entropy arguments, extended to additive isotonic regression in \citet{mammen2007additive}. For detailed convergence properties see e.g. \citet{zhang2002isotonicbounds}.
Notably, similar convergence rates hold in settings where the data is not strictly i.i.d., such as when using sample splitting or cross-fitting to estimate \(\hat{m}(\cdot)\). Theorem 2 in \cite{vanderlaan2024stabilized} demonstrates that under cross-fitting regimes, where \(\hat{m}(\cdot)\) is trained on an independent sample and applied to the estimation sample, the pseudo-sample $Z$ is still sufficiently weakly dependent for the \(O_P(N^{-1/3})\) rate to hold. This aligns with the mean squared error (MSE) decomposition in \eqref{eq:rmse_decomposition}, where the \(O_P(N^{-1/3})\) \(L_2(P)\) convergence rate of isotonic regression (Lemma \ref{lemma::convergence_isotonic_regression}) ensures that the calibration error term \(\|\mathbb{E}[m_0(X)|\hat{m}(X)] - \hat{m}(X)\|_{P,2}^2\) decays as \(O_P(N^{-2/3})\). Theorem 1 in \cite{vanderlaan2024stabilized} further establishes distribution-free calibration guarantees, bounding the calibration error by \(O_P(N^{-2/3})\), irrespective of the smoothness of the inverse propensity score or the dimension of the covariates. Together, these results accommodate the calibration algorithms introduced in Section \ref{sec:dml}.
\section{Calibration for Double Machine Learning}\label{sec:dml}
The following section combines calibration and double machine learning. The first part focuses on high-level conditions
for different double machine learning algorithms, whereas the second part states explicit conditions for particular double machine learning models with isotonic regression.
\subsection{Double Machine Learning Theory and Algorithms}
In this section, we state conditions which enable a re-estimation step for nuisance estimators in the double machine learning framework if the complexity of the re-estimation procedure is not too large.
As in \citet{Chernozhukov2018dml} we denote $\theta_0\in \Theta\subset \mathbb{R}$ the parameter of interest. The leading example is the average treatment effect (ATE) $\theta_0 = \mathbb{E} [ Y(1) - Y(0) ]$.
Further, we assume that $\theta_0$ satisfies the moment condition,
\begin{equation}\label{assumption::moment_condition}
\mathbb{E}_P[\psi(W;\theta_0, \eta_0)]=0,
\end{equation}
where $\psi$ is a known score function, the data $W$ is a random element in $(\mathcal{W}, \mathcal{A}_\mathcal{W})$ with probability measure $P\in\mathcal{P}_N$ and $\eta_0$ is the true value of the nuisance parameter $\eta\in T$, where $T$ is the convex subset of a normed vector space with norm $\|\cdot\|_T$.\\
The previous setting describes the standard double machine learning framework introduced in \cite{Chernozhukov2018dml}. For simplicity, we restrict ourselves to the case of linear score functions, that is
\begin{align}\label{assumption::linear_score}
\psi(w;\theta, \eta) = \psi^a(w;\eta)\theta + \psi^b(w;, \eta),
\end{align}
for all $w\in \mathcal{W}$, $\theta\in\Theta$ and $\eta\in T$.
Further, since we would like to consider a scenario with a re-estimation or calibration step, which might not affect all nuisance parameters, we define $\eta_0 = (\eta_0^{(1)}, \eta_0^{(2)})$,
where $\eta_0^{(2)}$ should be re-estimated as for example when $\eta_0^{(2)}$ is a propensity score to be calibrated. Correspondingly define $T = T^{(1)}\times T^{(2)}$.
Algorithm \ref{alg:dml2:standard} recaps the standard version of the double machine learning algorithm based on cross-fitting (cf. Definition 3.2 in \citet{Chernozhukov2018dml}).\\
Let $(W_i)_{i=1}^N$ be iid. copies of $W$ with probability measure $P$. To simplify notation, assume that $N$ is divisible by $K$.
\begin{algorithm}[!ht]
\caption{(uncalibrated) DML 2 Algorithm}
\setstretch{0.9}
\label{alg:dml2:standard}
{\footnotesize
\begin{algorithmic}[1]
\State \hspace*{0.02in} \textbf{Input:} Data $(W_i)_{i=1}^N$. A $K$-fold random partition $(I_k)_{k=1}^K$ of $[N]=\{1,\dots,N\}$ such that each fold $I_k$ is of size $n=N/K$. For each $k\in [K]=\{1,\dots,K\}$, define $I_k^c:=\{1,\dots,N\}\setminus I_k$.
\State For each $k\in[K]$, fit a machine learning estimator
$$\hat{\eta}_{0,k} = \hat{\eta}_{0}\big((W_i)_{i\in I_k^c}\big)$$
of $\eta_0$, where $\hat{\eta}_{0,k}$ is a random element in $T$, where the randomness only depends on the $(W_i)_{i\in I_k^c}$.
\State Construct the estimator $\hat{\theta}_0$ as the solution to
$$\frac{1}{K}\sum_{k=1}^K \mathbb{E}_{n,k}[\psi(W;\hat{\theta}_0, \hat{\eta}_{0,k})]=0,$$
where $\mathbb{E}_{n,k}[\psi(W)]=n^{-1}\sum_{i\in I_k} \psi(W_i)$ is the empirical expectation over $(W_i)_{i\in I_k}$.
\end{algorithmic}
}
\end{algorithm}
The standard DML 2 algorithm employs cross-fitting to handle the complexity of the estimated nuisance elements $\hat{\eta}_{0,k}$.
\begin{remark}\label{theorem::dml_asympt_dist}
Theorem 3.1 in \citet{Chernozhukov2018dml} shows that the estimator according to Algorithm \ref{alg:dml2:standard} is asymptotically normally distributed. More specifically, it holds
\begin{align}\label{eq:dml_normal_dist}
\sqrt{N}\sigma^{-1} (\hat{\theta}_0 - \theta_0) = \frac{1}{N}\sum_{i=1}^N \bar{\psi}(W_i) + O_P(\rho_N) \rightsquigarrow \mathcal{N}(0,1)
\end{align}
uniformly over $P\in\mathcal{P}_N$, where the size of the remainder term obeys $\rho_N:=N^{-1/2} + r_N + r'_N + N^{1/2}\lambda_N + N^{1/2}\lambda'_N \lesssim \delta_N$, with $\delta_N \ge N^{-1/2}$. Here, $\bar{\psi}(\cdot):=\sigma^{-1}J_0^{-1} \psi(\cdot; \theta_0, \eta_0)$ is the influence function and the approximate variance is $\sigma^2:= J_0^{-2}\mathbb{E}_P[\psi(W; \theta_0, \eta_0)^2]$.
\end{remark}
In Remark \ref{theorem::dml_asympt_dist} it is assumed that $\hat{\eta}_{0,k}= (\hat{\eta}_{0,k}^{(1)}, \hat{\eta}_{0,k}^{(2)}) \in \mathcal{T}_N$ with probability $1-o(1)$, where $\mathcal{T}_N$ is a suitable nuisance realization set.
To enable the re-estimation of nuisance elements $\hat{\eta}_{0,k}^{(2)}$ the algorithm and the nuisance realization set $\mathcal{T}_N$ has to be slightly adapted. A simple and straightforward adaption is presented in Algorithm \ref{alg:dml2:naive}.
\begin{algorithm}[!htb]
\caption{(nested $K$-fold cross-fitting calibration) DML 2 Algorithm}
\label{alg:dml2:naive}
{\footnotesize
\begin{algorithmic}[1]
\State \textbf{Input:} Data $(W_i)_{i=1}^N$. A $K$-fold random partition $(I_k)_{k=1}^K$ of $[N]=\{1,\dots,N\}$ such that each fold $I_k$ is of size $n=N/K$. For each $k\in [K]=\{1,\dots,K\}$ define $I_k^c:=\{1,\dots,N\}\setminus I_k$.
\State For each $k\in[K]$, fit a machine learning estimator
$$\hat{\eta}^{(1)}_{0,k} = \hat{\eta}^{(1)}_{0}\big((W_i)_{i\in I_k^c}\big)$$
of $\eta_0^{(1)}$, where $\hat{\eta}^{(1)}_{0,k}$ is a random element in $T^{(1)}$, where the randomness only depends on the $(W_i)_{i\in I_k^c}$.
\State For each $k\in[K]$, split the training partition $I_k^c$ into two disjoint samples $I_{k,1}^c$ and $I_{k,2}^c$.
\State Use the first subset $I_{k,1}^c$ to fit a machine learning estimator
$$\hat{\eta}^{(2)}_{0,k} = \hat{\eta}^{(2)}_{0}\big((W_i)_{i\in I_{k,1}^c}\big)$$
of $\eta^{(2)}_0$, where $\hat{\eta}^{(2)}_{0,k}$ is a random element in $T^{(2)}$, where the randomness only depends on the $(W_i)_{i\in I_{k,1}^c}$.
\State Use the second subset $I_{k,2}^c$ and the estimated nuisance element $\hat{\eta}_{0,k}^{(2)}$ to fit a re-estimation procedure
$$\tilde{\eta}^{(2)}_{0,k} =\tilde{\eta}^{(2)}_{0}\big((W_i)_{i\in I_{k,2}^c}, \hat{\eta}^{(2)}_{0,k}\big)$$
of $\eta_0^{(2)}$, where $\tilde{\eta}^{(2)}_{0,k}$ is a random element in $T^{(2)}$, where the randomness only depends on the $(W_i)_{i\in I_{k}^c}$.
\State Construct the estimator $\tilde{\theta}_0$ as the solution to
$$\frac{1}{K}\sum_{k=1}^K \mathbb{E}_{n,k}\big[\psi\big(W;\tilde{\theta}_0, (\hat{\eta}^{(1)}_{0,k}, \tilde{\eta}^{(2)}_{0,k})\big)\big]=0,$$
where $\mathbb{E}_{n,k}[\psi(W)]=n^{-1}\sum_{i\in I_k} \psi(W_i)$ is the empirical expectation over $(W_i)_{i\in I_k}$.
\end{algorithmic}
}
\end{algorithm}
As already mentioned, the standard double machine learning procedure in Algorithm \ref{alg:dml2:standard} uses cross-fitting to handle the complexity of estimated nuisance elements $\hat{\eta}_{0,k}$. Algorithm \ref{alg:dml2:naive} is a straightforward extension, which leaves the cross-fitting unchanged. The approach just employs a nested sample splitting procedure, such that the calibrated nuisance elements $\tilde{\eta}^{(2)}_{0,k}$ still depend only on the observations of the training sample $(W_i)_{i\in I_{k}^c}$. As a consequence, only the predictive performance of the calibrated nuisance estimators must be ensured.
\begin{remark}\label{theorem::dml_asympt_dist_naive}
Let $\tilde{\theta}_0$ be the estimator according to Algorithm \ref{alg:dml2:naive}.
Under the assumptions in Theorem 3.1 in \citet{Chernozhukov2018dml}, Equation \eqref{eq:dml_normal_dist} in Remark \ref{theorem::dml_asympt_dist} holds analogously for $\tilde{\theta}_0$.
\end{remark}
Typically, Assumption 3.2 in Theorem 3.1 in \citet{Chernozhukov2018dml} requires high-quality nuisance estimators. In particular, the re-estimation procedure (calibration) has to converge sufficiently fast, i.e. $
\|\tilde{\eta}^{(2)}_{0,k} -\eta^{(2)}_{0}\|_{P,2} \lesssim \varepsilon_N=o(N^{-1/4}).$
Since calibration properties are most important for small samples, further splitting of the training sample for calibration might not be desirable.
Therefore, we state Assumptions 3.1 and 3.2 of \citet{Chernozhukov2018dml} for an adapted nuisance realization set $\widetilde{\mathcal{T}}_N$. It is worth noting that in the following assumption the calibrated nuisance elements $\tilde{\eta}_{0}^{(2)}((W_i)_{i\in [N]})) )$ may depend on the full data.
This allows us to introduce new estimation algorithms that rely on more sophisticated splitting rules for calibration.
\begin{assumption}\label{assumption::dml}
Let \( c_0 > 0 \), \( c_1 > 0 \), and \( q \geq 2 \) be some finite constants such that \( c_0 \leq c_1\), and let \( \{\delta_N\}_{N \geq 1} \) and \( \{\Delta_N\}_{N \geq 1} \) be some sequences of positive constants converging to zero such that \( \delta_N \geq N^{-1/2} \). Also, let \( K \geq 2 \) be some fixed integer, and let \( \{\mathcal{P}_N\}_{N \geq 1} \) be some sequence of sets of probability distributions \( P \) of \(\mathcal{W}\) on \( W \).
\textbf{Assumption 3.1 (Linear scores with approximate Neyman orthogonality)} For all \( N \geq 3 \) and \( P \in \mathcal{P}_N \), the following conditions hold: \textup{(\romannumeral 1)} The true parameter value \( \theta_0 \) obeys (\ref{assumption::moment_condition}). \textup{(\romannumeral 2)} The score \( \psi \) is linear in the sense of (\ref{assumption::linear_score}). \textup{(\romannumeral 3)} The map \( \eta \mapsto E_P[\psi(W;\theta,\eta)] \) is twice continuously Gateaux-differentiable on \( T \). \textup{(\romannumeral 4)} The score \( \psi \) obeys the Neyman orthogonality or, more generally, the Neyman \( \lambda_N \) near-orthogonality condition at \( (\theta_0, \eta_0) \) with respect to the nuisance realization set \( \widetilde{\mathcal{T}}_N \subset T \) for $
\lambda_N := \sup_{\eta \in \mathcal{T}_N} \left| \partial_\eta \mathbb{E}_P[\psi(W;\theta_0,\eta_0)][\eta - \eta_0]\right| \leq \delta_N N^{-1/2}.$
\textup{(\romannumeral 5)} The identification condition holds; namely, the singular values of the matrix $J_0 := \mathbb{E}_P[\psi^a(W; \eta_0)]$ are between \( c_0 \) and \( c_1 \).
\bigskip
\textbf{Assumption 3.2 (Score regularity and quality of nuisance parameter estimators)} For all \( N \geq 3 \) and \( P \in \mathcal{P}_N \), the following conditions hold:
\begin{itemize}
\item[(a)] Given a random subset \( I \) of \( [N] \) of size \( n = N/K \), the nuisance parameter estimator \( \tilde{\eta}_0 := (\hat{\eta}_{0}^{(1)}((W_i)_{i\in I^c}), \tilde{\eta}_{0}^{(2)}((W_i)_{i\in [N]})) \) belongs to the realization set \( \widetilde{\mathcal{T}}_N \) with probability at least \( 1 - \Delta_N \), where \( \widetilde{\mathcal{T}}_N = \widetilde{\mathcal{T}}^{(1)}_N \times \widetilde{\mathcal{T}}^{(2)}_N \) contains \( \eta_0 \) and is constrained by the following conditions.
\item[(b)] The moment conditions hold:
\[
m_N := \sup_{\eta \in \widetilde{\mathcal{T}}_N} \mathbb{E}_P[|\psi(W; \theta_0, \eta)|^{q}]^{1/q} \leq c_1;
\quad
m_N' := \sup_{\eta \in \widetilde{\mathcal{T}}_N} \mathbb{E}_P[|\psi^a(W; \eta)|^{q}]^{1/q} \leq c_1.
\]
\item[(c)] The following conditions on the statistical rates \( r_N \), \( r_N' \), and \( \lambda_N' \) hold:
\[
r_N := \sup_{\eta \in \widetilde{\mathcal{T}}_N} \left|\mathbb{E}_P[\psi^a(W; \eta)] - \mathbb{E}_P[\psi^a(W; \eta_0)]\right| \leq \delta_N,
\]
\[
r_N' := \sup_{\eta^{(1)} \in \widetilde{\mathcal{T}}^{(1)}_N} \left(\mathbb{E}_P\left[\left(\psi\big(W; \theta_0, (\eta^{(1)}, \eta_0^{(2)})\big) - \psi\big(W; \theta_0, (\eta_0^{(1)}, \eta_0^{(2)})\big)\right)^2\right]\right)^{1/2} \leq \delta_N ,
\]
\[
\lambda_N' := \sup_{r\in(0,1), \eta \in \widetilde{\mathcal{T}}_N} \left|\partial^2_r \mathbb{E}_P[\psi(W; \theta_0, \eta_0 + r (\eta - \eta_0))]\right| \leq \delta_N / \sqrt{N}.
\]
\item[(d)] The variance of the score \( \psi \) is non-degenerate: $c_0\le \mathbb{E}_P[\psi(W; \theta_0, \eta_0)^2].$
\end{itemize}
\end{assumption}
Indeed, if the calibration method is not too complex (see Assumption \ref{assumption::dml_calibration_complexity}), for example when isotonic regression is used for calibration, the additional sample split in algorithm \ref{alg:dml2:naive} can be avoided by calibrating the predictions on each ``test''-fold $I_k$ which are used to estimate the target parameter $\theta_0$. This procedure is described in Algorithm \ref{alg:dml2:calibration}.
\begin{algorithm}[!htb]
\caption{($k$-fold cross-fitting calibration) DML 2 Algorithm}
\label{alg:dml2:calibration}
{\footnotesize
\begin{algorithmic}[1]
\State \textbf{Input:} Data $(W_i)_{i=1}^N$. A $K$-fold random partition $(I_k)_{k=1}^K$ of $[N]=\{1,\dots,N\}$ such that each fold $I_k$ is of size $n=N/K$. For each $k\in [K]=\{1,\dots,K\}$ define $I_k^c:=\{1,\dots,N\}\setminus I_k$.
\State For each $k\in[K]$, fit a machine learning estimator
$$\hat{\eta}_{0,k} = \hat{\eta}_{0}\big((W_i)_{i\in I_k^c}\big)$$
of $\eta_0$, where $\hat{\eta}_{0,k}$ is a random element in $T$, where the randomness only depends on the $(W_i)_{i\in I_k^c}$.
\State For each $k\in[K]$, rely on estimated nuisance element $\hat{\eta}_{0,k}^{(2)}$ to fit a re-estimation procedure
$$\tilde{\eta}^{(2)}_{0,k} =\tilde{\eta}^{(2)}_{0}\big((W_i)_{i\in I_{k}}, \hat{\eta}^{(2)}_{0,k}\big)$$
of $\eta_0^{(2)}$, where $\tilde{\eta}^{(2)}_{0,k}$ is a random element in $T^{(2)}$.
\State Construct the estimator $\tilde{\theta}_0$ as the solution to
$$\frac{1}{K}\sum_{k=1}^K \mathbb{E}_{n,k}\big[\psi\big(W;\tilde{\theta}_0, (\hat{\eta}^{(1)}_{0,k}, \tilde{\eta}^{(2)}_{0,k})\big)\big]=0,$$
where $\mathbb{E}_{n,k}[\psi(W)]=n^{-1}\sum_{i\in I_k} \psi(W_i)$ is the empirical expectation over $(W_i)_{i\in I_k}$.
\end{algorithmic}
}
\end{algorithm}
Although the calibrated nuisance estimator depends on the full data in this case, we will show an analog result as in Theorem 3.1 in \citet{Chernozhukov2018dml} under the Assumptions \ref{assumption::dml} and \ref{assumption::dml_calibration_complexity} in Theorem \ref{theorem::dml_asympt_dist_calibration}. As a slight modification of Algorithm \ref{alg:dml2:calibration} one can use different $K$-fold cross-fitting procedures for the estimated nuisance elements. For example, $2$-fold cross-fitting as described in Algorithm \ref{alg:dml2:calibration_split}, uses half of the data for nuisance estimation and the other half for calibration. Consequently, the calibration step might be more stable.
\begin{algorithm}[!htb]
\caption{(single split cross-fitting calibration) DML 2 Algorithm}
\label{alg:dml2:calibration_split}
{\footnotesize
\begin{algorithmic}[1]
\State \textbf{Input:} Data $(W_i)_{i=1}^N$. A $K$-fold random partition $(I_k)_{k=1}^K$ of $[N]=\{1,\dots,N\}$ such that each fold $I_k$ is of size $n=N/K$. A $2$-fold random partition $\tilde{I}_1\cup \tilde{I}_2 = \{1,\dots,N\}$ and $\tilde{I}_1\cap \tilde{I}_2 = \emptyset$. For each $k\in [K]=\{1,\dots,K\}$ define $I_k^c:=\{1,\dots,N\}\setminus I_k$ and for $j\in [2]=\{1,2\}$ define $\tilde{I}_j^c:=\{1,\dots,N\}\setminus \tilde{I}_j$.
\State For each $k\in[K]$, fit a machine learning estimator
$$\hat{\eta}^{(1)}_{0,k} = \hat{\eta}^{(1)}_{0}\big((W_i)_{i\in I_k^c}\big)$$
of $\eta_0^{(1)}$, where $\hat{\eta}^{(1)}_{0,k}$ is a random element in $T^{(1)}$, where the randomness only depends on the $(W_i)_{i\in I_k^c}$.
\State For each $j\in[2]$, fit a machine learning estimator
$$\hat{\eta}^{(2)}_{0,j} = \hat{\eta}^{(2)}_{0}\big((W_i)_{i\in \tilde{I}_j^c}\big)$$
of $\eta_0^{(2)}$, where $\hat{\eta}^{(2)}_{0,j}$ is a random element in $T^{(2)}$, where the randomness only depends on the $(W_i)_{i\in \tilde{I}_j^c}$.
\State For each $j\in[2]$, rely on estimated nuisance element $\hat{\eta}_{0,j}^{(2)}$ to fit a re-estimation procedure
$$\tilde{\eta}^{(2)}_{0,j} =\tilde{\eta}^{(2)}_{0}\big((W_i)_{i\in \tilde{I}_{j}}, \hat{\eta}^{(2)}_{0,j}\big)$$
of $\eta_0^{(2)}$, where $\tilde{\eta}^{(2)}_{0,j}$ is a random element in $T^{(2)}$.
\State Construct the estimator $\tilde{\theta}_0$ as the solution to
$$\frac{1}{2K}\sum_{k=1}^K\sum_{j=1}^2 \mathbb{E}_{n,k,j}\big[\psi\big(W;\tilde{\theta}_0, (\hat{\eta}^{(1)}_{0,k}, \tilde{\eta}^{(2)}_{0,j})\big)\big]=0,$$
where $\mathbb{E}_{n,k,j}[\psi(W)]=n^{-1}\sum_{i\in I_k\cap \tilde{I}_j} \psi(W_i)$ is the empirical expectation over $(W_i)_{i\in I_k \cap \tilde{I}_j}$.
\end{algorithmic}
}
\end{algorithm}
Another option is to simultaneously calibrate all cross-fitted predictions $\hat{\eta}_0^{(2)}$ as described in Algorithm \ref{alg:dml2:calibration_all}. The main difference between Algorithm \ref{alg:dml2:calibration} and \ref{alg:dml2:calibration_all} is the dependency structure of the data used to calibrate the nuisance elements. In Algorithm \ref{alg:dml2:calibration} the recalibration is fitted on i.i.d. samples conditional on the corresponding ``training''-fold $I_k^c$, whereas in Algorithm \ref{alg:dml2:calibration_all} samples used for the calibration step have a complex dependency structure.
\begin{algorithm}[!htb]
\caption{(full-sample calibration) DML 2 Algorithm}
\label{alg:dml2:calibration_all}
{\footnotesize
\begin{algorithmic}[1]
\State \textbf{Input:} Data $(W_i)_{i=1}^N$. A $K$-fold random partition $(I_k)_{k=1}^K$ of $[N]=\{1,\dots,N\}$ such that each fold $I_k$ is of size $n=N/K$. For each $k\in [K]=\{1,\dots,K\}$ define $I_k^c:=\{1,\dots,N\}\setminus I_k$.
\State For each $k\in[K]$, fit a machine learning estimator
$$\hat{\eta}_{0,k} = \hat{\eta}_{0}\big((W_i)_{i\in I_k^c}\big)$$
of $\eta_0$, where $\hat{\eta}_{0,k}$ is a random element in $T$, where the randomness only depends on the $(W_i)_{i\in I_k^c}$.
\State Combine all estimated nuisance elements $\hat{\eta}_{0,k}^{(2)}$ to fit a re-estimation procedure
$$\tilde{\eta}^{(2)}_{0} =\tilde{\eta}^{(2)}_{0}\big((W_i)_{i=1}^N, (\hat{\eta}^{(2)}_{0,k})_{k\in [K]}\big)$$
of $\eta_0^{(2)}$, where $\tilde{\eta}^{(2)}_{0}$ is a random element in $T^{(2)}$.
\State Construct the estimator $\tilde{\theta}_0$ as the solution to
$$\frac{1}{K}\sum_{k=1}^K \mathbb{E}_{n,k}\big[\psi\big(W;\tilde{\theta}_0, (\hat{\eta}^{(1)}_{0,k}, \tilde{\eta}^{(2)}_{0})\big)\big]=0,$$
where $\mathbb{E}_{n,k}[\psi(W)]=n^{-1}\sum_{i\in I_k} \psi(W_i)$ is the empirical expectation over $(W_i)_{i\in I_k}$.
\end{algorithmic}
}
\end{algorithm}
\begin{assumption}[Calibration Complexity]\label{assumption::dml_calibration_complexity}
Let \( \{\tilde{r}_N\}_{N \geq 1} \) and \( \{\tilde{r}^{a}_N\}_{N \geq 1} \) be some sequences of positive constants converging to zero. The estimator $\tilde{\eta} = (\hat{\eta}^{(1)}, \tilde{\eta}^{(2)})$ is a random element in $T$ with nuisance realization set $\widetilde{\mathcal{T}}_N\subseteq \widetilde{\mathcal{T}}^{(1)}_N \times \widetilde{\mathcal{T}}^{(2)}_N$ such that
\begin{enumerate}
\item[(i)]
Let $\eta^{(1)}$ is a fixed element in $\widetilde{\mathcal{T}}_N^{(1)}$ and define
\begin{align*}
\mathcal{F}_2(\eta^{(1)}):&=\left\{\psi\left(\cdot; \theta_0, (\eta^{(1)}, \eta^{(2)})\right)- \psi\left(\cdot; \theta_0, (\eta^{(1)}, \eta^{(2)}_0)\right)\big|(\eta^{(1)}, \eta^{(2)}) \in \widetilde{\mathcal{T}}_N\right\}.
\end{align*}
Further let $F_2(\eta^{(1)})(\cdot)$ be a measurable envelope for $\mathcal{F}_2(\eta^{(1)})$ such that $$\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F_2(\eta^{(1)})(W_i)\|_{P,q} \equiv \sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F_2(\eta^{(1)})\|_{P,q} \le V_n$$
for $q>2$. For each $\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}$ define $\sigma^2_n$ as a sequence converging to zero such that
\begin{align*}
&\sup_{f\in\mathcal{F}_2(\eta^{(1)})} \mathbb{E}[f^2] \le \sigma^2_n\le \|F_2(\eta^{(1)})\|_{P,2}^2
\end{align*}
and $u_n$ such that
\begin{align*}
&\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}} J\left(\sigma_n /\|F_2(\eta^{(1)})\|_{P,2} , \mathcal{F}_2(\eta^{(1)}), F_2(\eta^{(1)})\right) \le u_n.
\end{align*}
Finally, assume the following growth condition is satisfied
\begin{align*}
&\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}} \Bigg(u_n \|F_2(\eta^{(1)})\|_{P,2} + \sigma_n \sqrt{\log(n)} \\
\quad &+ n^{1/q-1/2}V_n \left(u_n^2\frac{\|F_2(\eta^{(1)})\|_{P,2} ^2}{\sigma_n^{2}}\vee\log(n)\right)\Bigg) \lesssim \tilde{r}_N
\end{align*}
\item[(ii)] It holds
\begin{align*}
\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\left|\mathbb{E}\left[\psi^a\left(W;(\eta^{(1)},\eta^{(2)}_{0})\right) - \psi^a\left(W;(\eta^{(1)}_0,\eta^{(2)}_{0})\right)\right]\right| \le \tilde{r}^a_N
\end{align*}
with probability converging to one.
Further, the entropy conditions above also need to hold for $\psi^a$:\\
Let $\eta^{(1)}$ is a fixed element in $\widetilde{\mathcal{T}}_N^{(1)}$ and define
\begin{align*}
\mathcal{F}^a_2(\eta^{(1)}):&=\left\{\psi^a\left(\cdot; (\eta^{(1)}, \eta^{(2)})\right)- \psi^a\left(\cdot; (\eta^{(1)}, \eta^{(2)}_0)\right)\big|(\eta^{(1)}, \eta^{(2)}) \in \widetilde{\mathcal{T}}_N\right\}.
\end{align*}
Further let $F^a_2(\eta^{(1)})(\cdot)$ be a measurable envelope for $\mathcal{F}^a_2(\eta^{(1)})$ such that $$\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F^a_2(\eta^{(1)})(W_i)\|_{P,q} \equiv \sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F^a_2(\eta^{(1)})\|_{P,q} \le V_{n,a}$$
for $q>2$. For each $\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}$ define $\sigma^2_{n,a}$ as a sequence converging to zero such that
\begin{align*}
&0<\sup_{f\in\mathcal{F}^a_2(\eta^{(1)})} \mathbb{E}[f^2] \le \sigma^2_{n,a}\le \|F^a_2(\eta^{(1)})\|_{P,2}^2
\end{align*}
and $u_{n,a}$ such that
\begin{align*}
&\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}} J\left(\sigma_{n,a} /\|F^a_2(\eta^{(1)})\|_{P,2} , \mathcal{F}^a_2(\eta^{(1)}), F^a_2(\eta^{(1)})\right) \le u_{n,a}.
\end{align*}
Finally, assume the following growth condition is satisfied
\begin{align*}
&u_{n,a} \sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F^a_2(\eta^{(1)})\|_{P,2} + \sigma_{n,a} \sqrt{\log(n)} \\
\quad &+ n^{1/q-1/2}V_{n,a} \left(u_{n,a}^2\frac{\sup_{\eta^{(1)}\in \widetilde{\mathcal{T}}_N^{(1)}}\|F^a_2(\eta^{(1)})\|_{P,2} ^2}{\sigma_{n,a}^{2}}\vee\log(n)\right) \lesssim \tilde{r}^a_N.
\end{align*}
\end{enumerate}
\end{assumption}
Assumption \ref{assumption::dml_calibration_complexity} imposes high-level assumptions on the complexity of the calibration step. Assumption \ref{assumption::dml_calibration_complexity} (i) restricts the complexity of the class $\mathcal{F}_2(\eta^{(1)})$ via standard complexity measures.
If the function class is suitably measurable and the uniform entropy integral obeys
$ \log \sup_Q N\left(\epsilon \|F_2(\eta^{(1)})\|_{Q,2}, \mathcal{F}_2(\eta^{(1)}), \|\cdot\|_{Q,2}\right)\le C $,
it holds
\begin{align*}
J\left(\sigma_n /\|F_2(\eta^{(1)})\|_{P,2} , \mathcal{F}_2(\eta^{(1)}), F_2(\eta^{(1)})\right)
&\lesssim \sigma_n^{1/2} \|F_2(\eta^{(1)})\|_{P,2}^{-1/2},
\end{align*}
since
\begin{align*}
J\left(\delta, \mathcal{F}_2(\eta^{(1)}), F_2(\eta^{(1)})\right)
:= & \int_0^{\delta} \sup_Q \sqrt{1 + \log N(\epsilon\|\mathcal{F}_2(\eta^{(1)})\|_{Q,2}, \mathcal{F}_2(\eta^{(1)}), \|\cdot\|_{Q,2}})d\epsilon \\
\le & \int_0^{\delta} \sqrt{1 + C \epsilon^{-1}}d\epsilon
\le \delta + \sqrt{C}\int_0^{\delta} \epsilon^{-1/2}d\epsilon
\lesssim \sqrt{\delta}
\end{align*}
for any $\delta$ small enough and probability measure $Q$. The first part of Assumption \ref{assumption::dml_calibration_complexity} (ii) imposes a Lipschitz continuity condition on $\psi^{a}$ which is the first part of the linear score defined in Equation \eqref{assumption::linear_score}.
The second part of Assumption \ref{assumption::dml_calibration_complexity} (ii) provides similar complexity assumptions as for the function class $\mathcal{F}_2(\eta^{(1)})$ in Assumption \ref{assumption::dml_calibration_complexity} (i) and also the required growth rates. Remark that the conditions in Assumption \ref{assumption::dml_calibration_complexity} are quite similar to \citet{belloni2018uniformly}, but we build upon standard Donsker conditions. Again, it is worth noting that Assumption \ref{assumption::dml_calibration_complexity} (ii) depends only on the score $\psi^{a}$. In the case of a nonparametric causal model, also known as \textit{interactive regression model} (IRM), considered in Section \ref{ch:irm}, the first part of the linear score is given by $\psi^{a}=-1$, see Equation \eqref{score_IRM}, and therefore $\sigma^2_{n,a}=0$. Hence, in the interactive regression model, Theorem \ref{theorem::dml_asympt_dist_calibration} below holds with $\tilde{r}^a_N=0$ and only Assumption \ref{assumption::dml_calibration_complexity} (i) is required.
\begin{theorem}\label{theorem::dml_asympt_dist_calibration}
Let $\tilde{\theta}_0$ be the estimator according to Algorithm \ref{alg:dml2:calibration} to \ref{alg:dml2:calibration_all}. Assume $\delta_N \ge N^{-1/2}$ for all $N\ge 1$.
Under Assumption \ref{assumption::dml} and \ref{assumption::dml_calibration_complexity}, equation \eqref{eq:dml_normal_dist} in Remark \ref{theorem::dml_asympt_dist} holds with updated remainder, that is
\begin{align*}
\sqrt{N}\sigma^{-1} (\hat{\theta}_0 - \theta_0) = \frac{1}{N}\sum_{i=1}^N \bar{\psi}(W_i) + O_P(\tilde{\rho}_N) \rightsquigarrow \mathcal{N}(0,1)
\end{align*}
uniformly over $P\in\mathcal{P}_N$, where the size of the remainder term obeys
\begin{align*}
\tilde{\rho}_N:= \rho_N + \tilde{r}_N + \tilde{r}^a_N = N^{-1/2} + r_N + r'_N + \tilde{r}_N + \tilde{r}^a_N + N^{1/2}\lambda_N + N^{1/2}\lambda'_N \lesssim \delta_N.
\end{align*}
\end{theorem}
The crucial difference between Algorithm \ref{alg:dml2:calibration} and \ref{alg:dml2:calibration_all} lies in the assumptions on the nuisance realization set $\mathcal{\tilde{T}}_N$, which require convergence rates of the recalibration procedure. In Algorithm \ref{alg:dml2:calibration}, estimation properties are well known, e.g. for isotonic regression see Section \ref{ch:isotonic_theory}. These proofs heavily rely on the i.i.d. assumption of the samples used for recalibration, which is violated for Algorithm \ref{alg:dml2:calibration_all}. Nevertheless, cross-fitting might result in only weak dependencies between different samples, such that the convergence rates might still be sufficient for the calibration with Algorithm \ref{alg:dml2:calibration_all}. This algorithm closely mirrors the IC-IPW approach described by \citet{vanderlaan2024stabilized}, where calibration is applied to cross-fitted propensity scores on the full sample, retaining the theoretical guarantees established in \citet{vanderlaan2024stabilized} (Theorems 1–2).
\subsection{Calibration for Double Machine Learning Models}\label{ch:calibration_dml_models}
The results of Section \ref{sec:dml} can be applied directly to different regression models and causal parameters of interest.
In standard settings, a convergence rate of $\|\hat{m}_0 - m_0\|_{P,2}= o_P(N^{-1/4})$ is assumed. If isotonic regression is used for calibration, Lemma \ref{lemma::convergence_isotonic_regression} directly implies that the convergence rate of the calibrated propensity score will still satisfy the rate condition $\|\tilde{m}_0 - m_0\|_{P,2}= o_P(N^{-1/4})$. Furthermore, we state explicit assumptions to restrict the complexity of the nuisance calibration to avoid additional sample splitting.
\subsubsection*{Calibration in interactive regression models}\label{ch:irm}
Consider the fully heterogeneous or interactive regression model as in \citet{Chernozhukov2018dml}. This nonparametric regression model is often considered when augmented inverse probability weighting (AIPW) is used for estimation. Let $D\in\{0,1\}$ be a binary treatment variable and $W=(Y,D,X)$, where
\begin{align}
Y = g_0(D,X) + U,\; \mathbb{E}[U|D,X] &= 0, &
D = m_0(X) + V,\; \mathbb{E}[V|X] &= 0. \label{IRM}
\end{align}
A common parameter of interest is the average treatment effect $\theta_0 := \mathbb{E}[g_0(1,X) - g_0(0,X)].$
Let the score function for the augmented inverse probability weighted estimator be
\begin{equation}\label{score_IRM}
\psi(W;\theta,\eta) := (g(1,X) - g(0,X)) + \frac{D}{m(X)}(Y - g(1,X)) - \frac{1-D}{1-m(X)}(Y - g(0,X)) -\theta
\end{equation}
where $\eta=(g,m)$ denotes the nuisance functions for the outcome regression $g_0(D,X) = E[Y|D,X]$ and propensity score $m_0(X)=E[D|X]$.
\begin{assumption}[cf. Assumption 5.1 in \citet{Chernozhukov2018dml}]\label{assumption::irm}
Let $\{\delta_N\},\{\Delta_N\}\searrow 0$; $c,\epsilon,C>0$, $q>4$, $K\geq 2$ fixed; $N/K\in\mathbb{N}$. For $\eta=(\eta_1,\dots,\eta_\ell)$, define $\|\eta\|_{p,q}:=\max_{1\leq j\leq\ell}\|\eta_j\|_{P,q}$. For all $P\in\mathcal{P}$, the following hold: (a) Equations \eqref{IRM} are satisfied, (b) $\|Y\|_{P,q}\leq C$, (c) $P(\epsilon\leq m_0(X)\leq1-\epsilon)=1$, (d) $\|U\|_{P,2}\geq c$, (e) $\|\mathbb{E}_P[U^2|X]\|_{P,\infty}\leq C$ and (f) for a random subset $I\subset[N]$ of size $n=N/K$, the nuisance parameter estimator $\hat{\eta}_0=\hat{\eta}_0((W_i)_{i\in I^c})$ satisfies, with $P$-probability $\geq1-\Delta_N$:
\begin{itemize}[leftmargin=2cm,noitemsep]
\item[\textup{(\romannumeral 1)}] $\|\hat{\eta}_0-\eta_0\|_{P,2}\leq\delta_N$, $\|\hat{\eta}_0-\eta_0\|_{P,q}\leq C$
\item[\textup{(\romannumeral 2)}] $\|\hat{m}_0-m_0\|_{P,2}\times\|\hat{g}_0-g_0\|_{P,2}\leq\delta_N N^{-1/2}$ with $\|\hat{g}_0-g_0\|_{P,\infty}\leq C$, $\|\hat{m}_0-1/2\|_{P,\infty}\leq1/2-\epsilon$
\end{itemize}
\end{assumption}
Under Assumption \ref{assumption::irm}, Remark \ref{theorem::dml_asympt_dist} holds (\citealt[Theorem 5.1]{Chernozhukov2018dml}). Our assumption strengthens the original by requiring $q>4$ rather than $q>2$. For Theorem \ref{theorem::dml_asympt_dist_calibration}, we add:
\begin{assumption}[Calibration rate/complexity]\label{assumption::irm_calibration_rate}
The following assumptions hold: \textup{(\romannumeral 1)} With $P$-probability $\geq 1-\Delta_N$, we have: $\|\tilde{m}(X)-m_0(X)\|_{P,2}\lesssim\varepsilon_N\leq\log^{-1/2}N$, $\varepsilon_N(\cdot)\|\hat{g}_0-g_0\|_{P,2}\leq\delta_N N^{-1/2}$. Further, the predictions are well separated from zero and one, $\|\tilde{m}(X)-1/2\|_{P,\infty}\leq1/2-\epsilon$. \textup{(\romannumeral 2)} Let $\tilde{m}(\cdot)\in\mathcal{M}$, such that the covering numbers obey $\sup_Q N(\epsilon,\mathcal{M},L_2(Q))\leq C\epsilon^{-1}$.
\end{assumption}
Assumption \ref{assumption::irm_calibration_rate} imposes mild conditions on the calibration procedure. Assumption \ref{assumption::irm_calibration_rate}(i) ensures that the convergence rate of the calibrated propensity score $\tilde{m}(\cdot)$ is still sufficiently fast, while Assumption \ref{assumption::irm_calibration_rate}(ii) restricts the complexity of calibration, so that additional cross-fitting can be avoided.
\begin{theorem}\label{theorem::irm}
Under Assumptions \ref{assumption::irm} and \ref{assumption::irm_calibration_rate}(i) Remark \ref{theorem::dml_asympt_dist_naive} is valid. If additionally Assumption \ref{assumption::irm_calibration_rate}(ii) is satisfied, Theorem \ref{theorem::dml_asympt_dist_calibration} holds.
\end{theorem}
The proof of Theorem \ref{theorem::irm} is given in the Supplementary Material.
As mentioned in the beginning of Section \ref{ch:calibration_dml_models}, often convergence rates of $\|\hat{m}_0 - m_0\|_{P,2}= o_P(N^{-1/4})$ and $\|\hat{g}_0 - g_0\|_{P,2}= o_P(N^{-1/4})$ are assumed implying the conditions of Assumption \ref{assumption::irm}(f)(ii). Considering Lemma \ref{lemma::convergence_isotonic_regression} this immediately implies Assumption \ref{assumption::irm_calibration_rate} if isotonic regression is used for calibration as the convergence rate of the calibrated propensity score is given by $\|\tilde{m}_0 - m_0\|_{P,2}= o_P(N^{-1/4})$ and the complexity of monotone functions satisfies Assumption \ref{assumption::irm_calibration_rate}(ii).
\section{Simulation Study}\label{sec:sim}
In this section, we investigate the introduced calibrated propensity score models from Section \ref{ch:calibration_dml_models} through an extensive simulation study\footnote{The code for the simulation study is available at the following link: \url{https://github.com/JanRabenseifner/Causal-Propensity-Calibration.git}. The simulation is executed on an HPC cluster in parallel, using different seeds for the DGPs.}. We evaluate the impact of calibration methods (Venn-ABERS, Platt scaling, isotonic regression) on the performance of causal estimators (IPW, DML), supplemented by analyses of weight normalization and a comparison to covariate-balancing reweighting estimators (e.g., entropy balancing). Performance is assessed using calibration diagnostics (e.g., calibration plots, expected calibration error) and causal estimation metrics (RMSE, MAE, and variance) to unravel the interplay between robustness, forecast accuracy, and covariate balance.
To assess the contribution of potentially miss-calibrated propensity scores, we briefly introduce the causal estimators considered. The inverse probability weighting (IPW) estimator uses estimates of the propensity scores $\hat{m}(D=1|X)$ directly. Here, an estimate $\hat{\theta}$ of the ATE is computed as $\frac{1}{n}\sum_{i=1}^n \bigg( \frac{D^{(i)}Y^{(i)}}{\hat{m}(D=1|X^{(i)})} - \frac{(1-D^{(i)})Y^{(i)}}{1-\hat{m}(D=1|X^{(i)})}\bigg)$.
Especially treated units with low propensity scores and non-treated units with high propensity scores have extreme contributions. This can be critical if the underlying propensity score model is misspecified or overconfident.
The interactive regression model (IRM, Section~\ref{ch:irm}) allows for heterogeneous treatment effects without strong form assumptions. Contrary, the partially linear regression model (PLR, Supplement \ref{S-ch:plm}) imposes an additive structure\footnote{Both, the IRM and PLR model, are implemented via the \texttt{DoubleML} package \citep{DoubleML2022Python, doubleml2024R}.}. One obstacle in evaluating causal ATE models is that the true value of the causal parameter $\theta_0$ is not observed in observational studies. For a fair evaluation, we have selected four external sources for the data generating processes (DGPs), each proposing different challenges for the models. An overview and detailed descriptions of the DGPs are provided in the Supplement \ref{S-appendix:extradgp}. The DGPs are characterized by varying levels of noise (DGP 1, \citet{belloni2017program}), different dimensionality of observed covariates (DGP 1), and underlying nonlinearities (DGP 2, \citet{deshpande2023calibrated}; DGP 3, \citet{van2023causal}), overlap violations (DGP 2), or unbalancedness (DGP 4, \citet{ballinari2024calibrating, nie2020unbalanced}). The DGPs satisfy the unconfoundedness assumption, $Y(d)\perp D \mid X$, with $Y(d)$ indicating the potential outcome under treatment $D=d$. Hence, these settings allow for identification of the average treatment effect, $\theta_0 = \mathbb{E}[Y(d=1)-Y(d=0)]$.
\subsection{Learners and Calibration Methods}
We test different learners for the outcome regression and the propensity score estimation. For the outcome regression, we consider a simple linear regression along with the tree-based Machine Learning algorithms LightGBM (LGBM) \citep{ke2017lightgbm} and random forest. Both are flexible machine learning algorithms that perform well across a wide variety of datasets.
For the propensity score estimation, we consider logistic regression, LGBM classifier, and random forest classifier. All models are employed within their default settings. Employing different fine-tuning schemes could benefit either approach and distort the comparison.
For propensity calibration, we consider the three approaches introduced in \ref{sec:calib_methods}. First, we utilize \texttt{IsotonicRegression} from the \texttt{scikit-learn} package \citep{scikit-learn}. In addition, we employ the Inductive Venn–ABERS predictor (VAP) introduced by \citet{vovk2015vennabers} and available at \citet{Petej2024vap}. VAP builds on the groupings in the outcome space made by isotonic regression. It utilizes potential labels to fit separate isotonic regressions. Thus, simple isotonic regression receives a coarser partitioning of the outcome space. Lastly, we incorporate Platt scaling, implemented via the \texttt{CalibratedClassifierCV} in the \texttt{scikit-learn} package.
\subsection{Calibration Metrics}
For binary classification, the $\ell_p$ Expected Calibration Error (ECE) \citep{naeini2015ece, sun2024l2ece}, for $p \ge 1$, is defined as:
$\label{ece} \mathrm{ECE}_{p} := \mathbb{E} \left[ \mathbb{E} \left[\|D - m(X)\|^p \mid m(X) \right] \right]^{\frac{1}{p}}.$ To approximate the expected calibration error (ECE), the estimated propensity scores $\hat{m}(X)$ are divided across the probabilistic output range $[0, 1]$ into equally spaced intervals \citep{naeini2015ece} $\{I_0, I_1, ..., I_M\}$ or quantiles \citep{nguyen2015posterior} of $\hat{m}(x)$. This allows us to generate buckets $\{B_i\}_{i=1}^M$, where $B_i = \{(X, D) \mid m(D=1 \mid X) \in I_i\}$. Each predicted probability is assigned to the appropriate bin. The calibration error is then defined as the difference between the fraction of correct predictions (accuracy) and the mean predicted probability (confidence) within each bin: $\operatorname{ECE}_p = \sum_{i=1}^{M} \frac{n_i}{N} \left\| \operatorname{acc}_i(B_i) - \operatorname{conf}_i(B_i) \right\|_p,$ where $\operatorname{acc}_i(B_i) = \frac{1}{|B_i|} \sum_{j=1}^{|B_i|} D_j$ and
$\operatorname{conf}_i(B_i) = \frac{1}{|B_i|} \sum_{j=1}^{|B_i|} m(D=1 \mid X_j)$.
\vspace{-0.5em}
\begin{wrapfigure}{!r}{0.45\textwidth}
\centering
\includegraphics[width=0.43\textwidth]{paper/figures/figures_irm/paper_calib_errors/n_obs_00000_dim_x_020_learner_m_LGBM_learner_g_LGBM_clip_0.00_R2d_0.50_overlap_0.00_share_treated_0.00.pdf}
\fbox{
\scriptsize
\begin{tabular}{@{}l@{\hspace{2.0em}}l@{\hspace{4.0em}}l@{\hspace{2.0em}}l@{}}
\begin{tikzpicture}[baseline=0.3ex]
\useasboundingbox (0,0.3) rectangle (2em,0.9ex);
\draw[color=cbgreen2, line width=1.5pt] (0,0.9ex) -- (2em,0.9ex);
\node[color=cbgreen2] at (1em,0.9ex) {\Large$\times$};
\end{tikzpicture}
& \raggedright\arraybackslash Alg-1-uncalib
& \begin{tikzpicture}[baseline=0.3ex]
\useasboundingbox (0,0.3) rectangle (2em,0.9ex);
\draw[color=cbblack, line width=1.5pt] (0,0.9ex) -- (2em,0.9ex);
\node[draw=cbblack, circle, line width=0.8pt, fill=cbblack, minimum size=3pt] at (1em,0.9ex) {};
\end{tikzpicture}
& \raggedright\arraybackslash Alg-2-nested-cf \\
\begin{tikzpicture}[baseline=0.3ex]
\useasboundingbox (0,0.3) rectangle (2em,0.9ex);
\draw[color=cbdarkblue, line width=1.5pt] (0,0.9ex) -- (2em,0.9ex);
\node[draw=cbdarkblue, diamond, line width=0.8pt, fill=cbdarkblue, minimum size=3pt] at (1em,0.9ex) {};
\end{tikzpicture}
& \raggedright\arraybackslash Alg-3-cf
& \begin{tikzpicture}[baseline=0.3ex]
\useasboundingbox (0,0.3) rectangle (2em,0.9ex);
\draw[color=cbbrown, line width=1.5pt] (0,0.9ex) -- (2em,0.9ex);
\node[draw=cbbrown, regular polygon, regular polygon sides=3, fill=cbbrown, line width=0.8pt, minimum size=3pt] at (1em,0.9ex) {};
\end{tikzpicture}
& \raggedright\arraybackslash Alg-4-single-split \\
\begin{tikzpicture}[baseline=0.3ex]
\useasboundingbox (0,0.3) rectangle (2em,0.9ex);
\draw[color=cbmagenta, line width=1.5pt] (0,0.9ex) -- (2em,0.9ex);
\node[draw=cbmagenta, regular polygon, regular polygon sides=3, shape border rotate=180, line width=0.8pt, fill=cbmagenta, minimum size=3pt] at (1em,0.9ex) {};
\end{tikzpicture}
& \multicolumn{3}{@{}l@{}}{\raggedright\arraybackslash Alg-5-full-sample} \\
\end{tabular}
}
\vspace{-0.2em}
\caption{Quantile ECE, DGP 1, \\m = LGBM, n =2000, p = 20}
\label{fig:ece_n_obs}
\vspace{-1em}
\end{wrapfigure}
In Figure \ref{fig:ece_n_obs}, we can observe that the uncalibrated Algorithm \ref{alg:dml2:standard}, as well as the nested k-fold cross-fit Algorithm \ref{alg:dml2:naive} are poorly calibrated for small sample sizes. Additionally, a version of Algorithm \ref{alg:dml2:standard} clipped at the one percent level is included. This helps neglect some of the miscalibration, but still performs worse than Algorithms \ref{alg:dml2:calibration_split} and \ref{alg:dml2:calibration_all} for all sample sizes. In propensity weighting, severe deviations in the middle of the propensity score distribution are not particularly critical. However, deviations at the boundaries are crucial because they can lead to exploding weights. Therefore, it is generally advisable to include visualizations to assess both the overlap in propensity scores and their calibration properties. The overlap ratio plot, based on the reliability diagram, splits the probability space into equal parts. For each propensity bin, the plot displays the actual proportion of treated and untreated units separately.
\begin{figure}[!h]
\centering
\includegraphics[width=0.9\textwidth]{paper/figures/figures_irm/paper_ps_treatment/n_obs_02000_dim_x_020_learner_m_000_learner_g_000_clip_0.00_R2d_0.50_overlap_0.00_share_treated_0.00.pdf}
\vspace{-0.5em}
\fbox{
\scriptsize
\begin{tabular}{@{}l@{\hspace{1em}}llrr}
\textcolor{cbblue3}{\rule[+0.4ex]{2em}{3.5pt}} & Treated &
\textcolor{cborange2}{\rule[+0.4ex]{2em}{3.5pt}} & Untreated
\end{tabular}
}
\vspace{-0.2em}
\caption{Overlap Ratios, DGP 1, n = 2000, p = 20, m = LGBM}
\label{fig:overlap_ratio}
\end{figure}
\vspace{-1em}
No clear violation of the overlap assumption can be seen for the true underlying propensities displayed in the left panel of Figure \ref{fig:overlap_ratio}. The black dotted lines represent perfect calibration. The ratios illustrate the deviations for both the treated and untreated groups in the uncalibrated Algorithm \ref{alg:dml2:standard} in the middle panel. Notably critical are the substantial proportions of treated observations with estimated propensity scores near zero and untreated observations with propensities close to one. In contrast, the perfect calibration property of Algorithm \ref{alg:dml2:calibration_all} is displayed in the right panel.
\subsection{General Findings}
As expected, the nested cross-fitting Algorithm \ref{alg:dml2:naive} faces stability challenges in small sample size settings. Venn-Abers calibration relies on isotonic regression combined with sample splitting. Consequently, Algorithm \ref{alg:dml2:calibration}, when used with Venn-Abers, encountered similar instability. Additionally, the two-fold calibration Algorithm \ref{alg:dml2:calibration_split}, combined with isotonic regression, required clipping at the 1-percent level to maintain stability. This instability arises from the known limitation of isotonic regression, which is prone to overfitting, especially with small calibration sets \citep{vanderlaan2024selfcalibrating}. To ensure a fair comparison, the uncalibrated Algorithm \ref{alg:dml2:standard} is presented both unclipped and with a restriction at the 1-percent threshold (Alg-1-Clipped). Algorithms \ref{alg:dml2:calibration} and \ref{alg:dml2:calibration_all} were only clipped at a threshold of $10^{-12}$. Given the sample sizes, such violations seem unlikely. This serves more as a general recommendation, as the added clipping bias is negligible. Table \ref{tab:results_overview} provides a summary of the results across all treatment models and the algorithms discussed, in combination with isotonic regression.
\begin{table}[H]
\begin{center}
\begin{footnotesize}
\begin{singlespace}
\captionsetup{justification=centering}
\caption{Results Overview}
\label{tab:results_overview}
\scalebox{0.85}{
\centering
\begin{tabular}{lllccccccccc}
\toprule
DGP & Model & Method & \multicolumn{3}{c}{m = Logit} & \multicolumn{3}{c}{m = Random Forest} & \multicolumn{3}{c}{m = LGBM} \\
& & & MAE & RMSE & Std. dev. & MAE & RMSE & Std. dev. & MAE & RMSE & Std. dev. \\
\midrule
\multirow{12}{*}{1} & IRM & Alg-1-Clipped & 0.07 & 0.10 & 0.10 & 0.06 & 0.08 & 0.07 & 0.22 & 0.27 & 0.24 \\
& IRM & Alg-1-Uncalib & 0.08 & 0.13 & 0.13 & 1.85e+06 & 1.85e+07 & 1.84e+07 & 0.48 & 0.60 & 0.54 \\
& IRM & Alg-2-nested-cf-Iso & 0.12 & 0.15 & 0.14 & 0.10 & 0.12 & 0.12 & 0.10 & 0.12 & 0.12 \\
& IRM & Alg-3-cf-Iso & 0.06 & 0.08 & 0.07 & 0.06 & 0.08 & 0.08 & 0.06 & 0.08 & 0.07 \\
& IRM & Alg-4-single-split-Iso & 0.10 & 0.13 & 0.12 & 0.07 & 0.08 & 0.06 & 0.07 & 0.09 & 0.06 \\
& IRM & Alg-5-full-sample-Iso & 0.06 & 0.08 & 0.08 & 0.06 & 0.08 & 0.07 & 0.07 & 0.08 & 0.08 \\
\cline{2-12}
& IPW & Alg-1-Clipped & 0.08 & 0.11 & 0.11 & 0.18 & 0.19 & 0.06 & 0.94 & 0.98 & 0.31 \\
& IPW & Alg-1-Uncalib & 0.10 & 0.16 & 0.16 & 1.09e+06 & 1.09e+07 & 1.09e+07 & 1.57 & 1.77 & 0.83 \\
& IPW & Alg-5-full-sample & 0.07 & 0.09 & 0.08 & 0.12 & 0.14 & 0.06 & 0.14 & 0.16 & 0.07 \\
\cline{2-12}
& PLR & Alg-1-Clipped & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 \\
& PLR & Alg-1-Uncalib & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 \\
& PLR & Alg-5-full-sample & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.05 & 0.07 & 0.07 \\
\midrule
\multirow{12}{*}{2} & IRM & Alg-1-Clipped & 0.09 & 0.11 & 0.11 & 0.31 & 0.39 & 0.39 & 0.20 & 0.26 & 0.26 \\
& IRM & Alg-1-Uncalib & 0.09 & 0.11 & 0.11 & 2.22e+09 & 2.89e+09 & 2.86e+09 & 0.24 & 0.32 & 0.32 \\
& IRM & Alg-2-nested-cf-Iso & 0.18 & 0.23 & 0.22 & 0.16 & 0.22 & 0.22 & 0.19 & 0.24 & 0.24 \\
& IRM & Alg-3-cf-Iso & 0.09 & 0.11 & 0.11 & 0.09 & 0.12 & 0.12 & 0.10 & 0.12 & 0.12 \\
& IRM & Alg-4-single-split-Iso & 0.15 & 0.20 & 0.20 & 0.09 & 0.11 & 0.11 & 0.09 & 0.11 & 0.11 \\
& IRM & Alg-5-full-sample-Iso & 0.09 & 0.11 & 0.11 & 0.09 & 0.11 & 0.11 & 0.09 & 0.12 & 0.12 \\
\cline{2-12}
& IPW & Alg-1-Clipped & 0.09 & 0.12 & 0.11 & 4.45 & 4.59 & 1.11 & 2.67 & 2.74 & 0.61\\
& IPW & Alg-1-Uncalib & 0.09 & 0.12 & 0.11 & 1.17e+10 & 1.37e+10 & 7.23e+09 & 2.86 & 2.98 & 0.81 \\
& IPW & Alg-5-full-sample-Iso & 0.09 & 0.11 & 0.11 & 0.20 & 0.23 & 0.10 & 0.18 & 0.20 & 0.10\\
\cline{2-12}
& PLR & Alg-1-Clipped & 0.09 & 0.11 & 0.10 & 0.12 & 0.14 & 0.10 & 0.10 & 0.11 & 0.10 \\
& PLR & Alg-1-Uncalib & 0.09 & 0.11 & 0.10 & 0.12 & 0.14 & 0.10 & 0.10 & 0.11 & 0.10 \\
& PLR & Alg-5-full-sample-Iso & 0.09 & 0.11 & 0.10 & 0.09 & 0.11 & 0.10 & 0.09 & 0.10 & 0.10 \\
\hline
\multirow{12}{*}{3} & IRM & Alg-1-Clipped & 0.05 & 0.07 & 0.07 & 0.08 & 0.10 & 0.10 & 0.11 & 0.13 & 0.13 \\
& IRM & Alg-1-Uncalib & 0.05 & 0.07 & 0.07 & 9.97e+07 & 2.74e+08 & 2.73e+08 & 0.11 & 0.14 & 0.14 \\
& IRM & Alg-2-nested-cf-Iso & 0.10 & 0.13 & 0.12 & 0.08 & 0.10 & 0.10 & 0.10 & 0.13 & 0.13 \\
& IRM & Alg-3-cf-Iso & 0.05 & 0.07 & 0.07 & 0.06 & 0.07 & 0.07 & 0.05 & 0.07 & 0.06 \\
& IRM & Alg-4-single-split-Iso & 0.10 & 0.12 & 0.11 & 0.05 & 0.07 & 0.06 & 0.05 & 0.07 & 0.06 \\
& IRM & Alg-5-full-sample-Iso & 0.05 & 0.07 & 0.07 & 0.06 & 0.07 & 0.07 & 0.05 & 0.07 & 0.07 \\
\cline{2-12}
& IPW & Alg-1-Clipped & 0.06 & 0.08 & 0.08 & 0.54 & 0.58 & 0.21 & 2.00 & 2.02 & 0.31\\
& IPW & Alg-1-Uncalib & 0.06 & 0.08 & 0.08 & 3.07e+08 & 8.49e+08 & 7.91e+08 & 2.05 & 2.07 & 0.34 \\
& IPW & Alg-5-full-sample-Iso & 0.06 & 0.08 & 0.08 & 0.37 & 0.37 & 0.07 & 0.42 & 0.43 & 0.07\\
\cline{2-12}
& PLR & Alg-1-Clipped & 0.06 & 0.08 & 0.08 & 0.06 & 0.07 & 0.07 & 0.08 & 0.10 & 0.07 \\
& PLR & Alg-1-Uncalib & 0.06 & 0.08 & 0.08 & 0.06 & 0.07 & 0.07 & 0.08 & 0.10 & 0.07 \\
& PLR & Alg-5-full-sample-Iso & 0.07 & 0.09 & 0.08 & 0.06 & 0.07 & 0.07 & 0.06 & 0.07 & 0.07 \\
\hline
\multirow{12}{*}{4} & IRM & Alg-1-Clipped & 0.04 & 0.06 & 0.06 & 0.07 & 0.09 & 0.09 & 0.17 & 0.21 & 0.20 \\
& IRM & Alg-1-Uncalib & 0.04 & 0.06 & 0.06 & 1.74e+08 & 2.70e+08 & 2.68e+08 & 0.36 & 0.46 & 0.43 \\
& IRM & Alg-2-nested-cf-Iso & 0.05 & 0.07 & 0.07 & 0.05 & 0.07 & 0.07 & 0.06 & 0.07 & 0.07 \\
& IRM & Alg-3-cf-Iso & 0.04 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 \\
& IRM & Alg-4-single-split-Iso & 0.05 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.04 & 0.06 & 0.06 \\
& IRM & Alg-5-full-sample-Iso & 0.04 & 0.06 & 0.06 & 0.05 & 0.06 & 0.06 & 0.04 & 0.06 & 0.06 \\
\cline{2-12}
& IPW & Alg-1-Clipped & 0.13 & 0.15 & 0.06 & 0.54 & 0.55 & 0.13 & 6.14 & 6.17 & 0.52 \\
& IPW & Alg-1-Uncalib & 0.13 & 0.15 & 0.06 & 4.64e+08 & 6.89e+08 & 5.11e+08 & 8.37 & 8.45 & 1.21 \\
& IPW & Alg-5-full-sample-Iso & 0.07 & 0.08 & 0.05 & 0.11 & 0.12 & 0.05 & 0.10 & 0.11 & 0.06 \\
\cline{2-12}
& PLR & Alg-1-Clipped & 0.08 & 0.09 & 0.06 & 0.05 & 0.06 & 0.05 & 0.05 & 0.06 & 0.05 \\
& PLR & Alg-1-Uncalib & 0.08 & 0.09 & 0.06 & 0.05 & 0.06 & 0.05 & 0.05 & 0.06 & 0.05 \\
& PLR & Alg-5-full-sample-Iso & 0.08 & 0.09 & 0.06 & 0.07 & 0.09 & 0.05 & 0.08 & 0.09 & 0.05 \\
\bottomrule
\end{tabular}
}
\caption*{\scriptsize For all DGPs: g = LGBM, and for Algorithms 2 - 5: Calibration = Isotonic Regression; DGP 1: n = 2000, p = 20, R2\_d = 0.5; DGP 2: n = 2000, p = 3, overlap = 0.5; DGP 3: n = 2000, p = 4; DGP 4: n = 4000, p = 20, share treated = 0.1}
\end{singlespace}
\end{footnotesize}
\end{center}
\end{table}
\vspace{-1.5cm}
Across all DGPs and settings, we can observe that calibration improves the inverse propensity-based IPW and IRM especially in combination with the tree-based propensity learners. The PLR model produces stable results across all settings, with minimal improvement from clipping or calibration for most DGPs. In general, the doubly-robust calibrated IRM and the PLR outperform the IPW. Algorithm \ref{alg:dml2:calibration_split} is biased in the PLR model for the tree-based methods random forest and LGBM in combination with VAP or isotonic regression\footnote{For more details on the influence of sample size on the proposed Algorithms, we refer to the Supplement \ref{S-appendix:extrasim}, Figures \ref{S-fig:ate_errors_n_irm}, \ref{S-fig:ate_errors_n_drug}, \ref{S-fig:ate_errors_n_nl}, \ref{S-fig:ate_errors_n_unb}.}. This bias appears to persist regardless of the sample size, as demonstrated on the right-hand side of Figure \ref{fig:ate_errors_n}.
\vspace{-0.5em}
\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{paper/figures/figures_irm/paper_ate_errors_IPW_IRM_PLR/n_obs_00000_dim_x_020_learner_m_LGBM_learner_g_LGBM_clip_0.01_R2d_0.50_overlap_0.00_share_treated_0.00.PDF}
\vspace{-0.5em}
\fbox{
\scriptsize
\begin{tabular}{@{}l@{\hspace{0.8em}}llllllllll}
\raisebox{0.5ex}{\begin{tikzpicture}[baseline=-0.5ex]
\draw[color=cbgreen2, line width=1.5pt] (0,0) -- (2em,0);
\node[color=cbgreen2] at (1em,0) {\Large$\times$};
\end{tikzpicture}} & Alg-1-uncalib &
\raisebox{0.5ex}{\begin{tikzpicture}[baseline=-0.5ex]
\draw[color=cbblack, line width=1.5pt] (0,0) -- (2em,0);
\node[draw=cbblack, circle, line width=0.8pt, fill=cbblack, minimum size=3pt] at (1em,0) {};
\end{tikzpicture}} & Alg-2-nested-cf &
\raisebox{0.5ex}{\begin{tikzpicture}[baseline=-0.5ex]
\draw[color=cbdarkblue, line width=1.5pt] (0,0) -- (2em,0);
\node[draw=cbdarkblue, diamond, line width=0.8pt, fill=cbdarkblue, minimum size=3pt] at (1em,0) {};
\end{tikzpicture}} & Alg-3-cf &
\raisebox{0.5ex}{\begin{tikzpicture}[baseline=-0.5ex]
\draw[color=cbbrown, line width=1.5pt] (0,0) -- (2em,0);
\node[draw=cbbrown, regular polygon, regular polygon sides=3, fill=cbbrown, line width=0.8pt, minimum size=3pt] at (1em,0) {};
\end{tikzpicture}} & Alg-4-single-split &
\raisebox{0.5ex}{\begin{tikzpicture}[baseline=-0.5ex]
\draw[color=cbmagenta, line width=1.5pt] (0,0) -- (2em,0);
\node[draw=cbmagenta, regular polygon, regular polygon sides=3, shape border rotate=180, line width=0.8pt,fill=cbmagenta, minimum size=3pt] at (1em,0) {};
\end{tikzpicture}} & Alg-5-full-sample \\
\end{tabular}
}
\vspace{-0.2em}
\caption{DGP 1, n =2000, p = 20, R2D = 0.5, m = LGBM, g = LGBM}
\label{fig:ate_errors_n}
\end{figure}
\vspace{-1.5em}
The impact of calibration is strongly dependent on the underlying propensity score learner. In contrast to \citet{blasiok2023doesoptimizingproperloss}, we generally observe good calibration properties for logistic regression and the least improvements through calibration.
\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{paper/figures/figures_irm/paper_ps_calib/n_obs_02000_dim_x_020_covariate_x_1_clip_0.00_R2d_0.50_overlap_0.00_share_treated_0.00.pdf}
\vspace{-0.5em}
\fbox{
\scriptsize
\begin{tabular}{@{}l@{\hspace{1em}}ccccc}
\textcolor{cbblue2}{\tikz[baseline=-0.7ex]{\fill (0,0) circle (0.4em);}} & Observed treatment &
\textcolor{cbgreen2}{\rule[+0.5ex]{2em}{1.5pt}} & True $m_0$ &
\textcolor{cbbrown}{\rule[+0.5ex]{2em}{1.5pt}} & Predicted $\hat{m}$
\end{tabular}
}
\vspace{-0.2em}
\caption{DGP 1, n =2000, p = 20, R2D = 0.5}
\label{fig:ps_calib_iso_2_irm}
\end{figure}
\vspace{-1.5em} Random forest can be both under-confident, as seen in Figure \ref{fig:ps_calib_iso_2_irm} for DGP 1, and over-confident for DGPs 2 and 3\footnote{The corresponding figures are located in the Supplementary Material, Figures \ref{S-fig:ps_calib_iso_2_drug} and \ref{S-fig:ps_calib_iso_2_nl}.}. As shown by \citet{johansson2023wellcalibrated}, random forest tends to be under-confident for the minority class in unbalanced settings (DGP 4, Figure \ref{S-fig:ps_calib_iso_2_drug} in the Supplementary Material). In general, the combination of random forest with calibration performs well across various settings. On the other hand, boosting-based methods, such as LGBM, tend to be over-confident.
Figure \ref{fig:learner_m} displays the distribution of ATE estimates for different propensity score learners under the IRM across 100 repetitions with different seeds in DGP 1. Calibration can correct for both under-confident and over-confident learners. However, the impact appears strongest for over-confident learners. For more details, we refer to the Supplementary Material, where the robustness of our algorithms is tested with respect to the propensity and outcome learners, different clipping thresholds, and various levels of signal-to-noise ratio (DGP 1), overlap (DGP 2), and share of treated units (DGP 4).
\vspace{-0.5em}
\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{paper/figures/figures_irm/paper_learner_m/n_obs_02000_dim_x_020_learner_m_003d_learner_g_LGBM_clip_0.01_R2d_0.50_overlap_0.00_share_treated_0.00.pdf}
\vspace{-0.7em}
\fbox{
\scriptsize
\begin{tabular}{@{}l@{\hspace{1em}}llllll}
\textcolor{cbblue}{\rule[+0.5ex]{2em}{1.5pt}} & Alg-1-uncalib &
\textcolor{cborange}{\rule[+0.5ex]{2em}{1.5pt}} & Alg-2-nested-cf &
\textcolor{cbgreen}{\rule[+0.5ex]{2em}{1.5pt}} & Alg-3-cf \\
\textcolor{cbred}{\rule[+0.5ex]{2em}{1.5pt}} & Alg-4-single-split &
\textcolor{cbpurple}{\rule[+0.5ex]{2em}{1.5pt}} & Alg-5-full-sample &
\textcolor{cbyellow}{\rule[+0.5ex]{2em}{1.5pt}} & Oracle \\
\end{tabular}
}
\caption{DGP 1, n =2000, p = 20, R2D = 0.5, m = LGBM, g = LGBM}
\label{fig:learner_m}
\end{figure}
\vspace{-1.5em}
\section{Discussion}
We extended and adapted some of the simulation settings to gain further insights from the studies implemented by \citet{deshpande2023calibrated}, \citet{Gutman2024platt}, \citet{ballinari2024calibrating}, and \citet{vanderlaan2024stabilized}.
In our simulation study, our objective was to determine whether calibration works and to explore the best methods to achieve effective calibration.
\citet{deshpande2023calibrated} employ a non-cross-fitted version of Algorithm \ref{alg:dml2:calibration_split}, where the data is split for calibration only. They achieve significant improvements, even in logistic regression, through calibration. This result seems counterintuitive. However, given their deterministic propensity scores in the drug effectiveness DGP, it is unsurprising that logistic regression is miscalibrated. In our adapted version, DGP 2, with non-deterministic propensity scores, such large improvements are not observed. In particular, Algorithm \ref{alg:dml2:calibration_split} without clipping, or in combination with PLR, is not recommended based on our findings.
\citet{ballinari2024calibrating} implements the nested cross-fitting Algorithm \ref{alg:dml2:naive}. As demonstrated, nested cross-fitting combined with isotonic regression is effective only when clipping is applied. Generally, Algorithm \ref{alg:dml2:naive} shows poor calibration properties and is not recommended for small sample sizes, where calibration has the most significant impact. The unbalanced and nonlinear settings labeled "difficult" and "extreme" for DGPs 4 and 5 by \citet{ballinari2024calibrating} were tested under DGP 4 in this study. We show that either using a larger share of observations for calibration or adding 1-percent clipping allows isotonic regression to lower RMSEs for boosting-based methods and remain stable for other learners. The instability observed in the results of \citet{ballinari2024calibrating} was also addressed by \citet{vanderlaan2024stabilized}. The authors concluded that, in a discretized version, it is important to calibrate the treated and untreated observations separately. However, as we demonstrate, the instability of isotonic regression is more likely due to small sample size issues. Our calibration algorithms, \ref{alg:dml2:calibration} and \ref{alg:dml2:calibration_all}, are both stable without clipping or separate calibration for treated and untreated units. The latter follows \citet{vanderlaan2024stabilized}'s recommendation to calibrate on the full sample using cross-fitted propensity scores.
In summary, our findings emphasize the critical role of method selection and sample size in calibration procedures. While (nested) cross-fitting for the calibration step is not always necessary, it may require supplemental clipping in settings prone to overfitting. Crucially, propensity score calibration enhances the robustness of inverse propensity-weighted ATE estimates.
\section{Acknowledgment}
This work was partly funded by the Bavarian Joint Research Program (BayVFP) – Digitization
(Funding reference: DIK0294/01). Economic AI kindly thanks the the VDI/VDE-IT Munich for the organization and the Free State of Bavaria for the financial support.
\bibliographystyle{abbrvnat}
\bibliography{paper/paper-ref}
\clearpage
\ifarXiv
\foreach \x in {1,...,\the\pdflastximagepages}
{
\includepdf[pages={\x}, fitpaper=true]{supplement.pdf}
}
\fi