Density Ratio-Free Doubly Robust Proxy Causal Learning
Authors: Bariscan Bozkurt
Houssam Zenati
Dimitri Meunier
Liyuan Xu
Arthur Gretton
Venue: NeurIPS 2025 poster ()
Keywords: Proxy causal learning
doubly robust estimation
kernel mean embeddings
Primary area: probabilistic_methods

Abstract: We study the problem of causal function estimation in the Proxy Causal Learning (PCL) framework, where  confounders are not observed but proxies for the confounders are available. Two main approaches have been proposed: outcome bridge-based and treatment bridge-based methods. In this work, we propose two kernel-based doubly robust estimators that combine the strengths of both approaches, and naturally handle continuous and high-dimensional variables. Our identification strategy builds on a recent density ratio-free method for treatment bridge-based PCL; furthermore, in contrast to previous approaches, it does not require indicator functions or kernel smoothing over the treatment variable. These properties make it especially well-suited for continuous or high-dimensional treatments. By using kernel mean embeddings, we propose the first density-ratio free doubly robust estimators for proxy causal learning, which have closed form solutions and strong uniform consistency guarantees. Our estimators outperform existing methods on PCL benchmarks, including a prior doubly robust method that requires both kernel smoothing and density ratio estimation.

Density Ratio-Free Doubly Robust Proxy Causal Learning
Bariscan Bozkurt!
Houssam Zenati
Dimitri Meunier!
Liyuan Xu²
Arthur Gretton
'Gatsby Computational Neuroscience Unit, University College London, 2Secondmind, 'DeepMind.
{bariscan. bozkurt. 23, h. h.zenati, dimitri dimitri.meunier mewnier. 21l-211...e.. ac.uk 1iyuan99880gmaill com arthur arthur. grathrrgggetttn@gaall.... com
Abstract
We study the problem of causal function estimation in the Proxy Causal Learn- ing (PCL.) framework, where confounders are not observed but proxies for the confounders are available. Two main approaches have been proposed: outcome bridge-based and treatment bridge-based methods. In this work, we propose kernel-based doubly robust estimators that combine the strengths of both ap- proaches, and naturally handle continuous and nigh-dimensional variables. Our identification strategy builds on a recent density ratio-free method for treatment bridge- based PCL; furthermore, in contrast to previous approaches, it does not require indicator functions or kernel smoothing over the treatment variable. These properties make it especially well-suited for continuous or high-dimensional treat- ments. By using kernel mean embeddings, we propose the first density-ratio free doubly robust estimators for proxy causal learning, which have closed form so- lutions and strong uniform consistency guarantees. Our estimators outperform existing methods on PCL benchmarks, including prior doubly robust method that requires both kernel smoothing and density ratio estimation.
Introduction and related works
Estimating the effects of interventions-alss referred to treatments-on outcomes is central goal of causal learning. This task is particularly challenging in observational settings, where randomized experiments
feasible. major difficulty stems from confounding variables that influence both the treatment and the outcome. In practice, the commonly made unconfoundedness (or ignorability) assumption that there are no unobserved confounders- -otendderr-ottnnnn fails hold, as it is unrealistic to expect that all relevant confounders can be accounted for. As a result, one may instead assume that the observed covariates serve as proxies for latent unmeasured confounders.
One classical line work for addressing unobserved confounding is instrumental variable (IV) regression, which assumes access instruments that affect the treatment but independent of the unobserved confounders 2
A more recent and promising direction is the Proxy Causal Learning (PCL), which leverages two auxiliary variables: (i) treatment proxy, denoted by Z. which is causally related to the treatment, and (ii) an outcome proxy, denoted by W, which is causally related to the outcome. The associated causal graph illustrated in Figure ( In this setup, bidirectional arrow indicates that either causal direction between the two variables is plausible, or that they may share an unobserved common cause. Within this framework, Miao et al. 5 showed that the causal effect can be identified via an outcome bridge function, without the need
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

to explicitly recover the latent confounders. This contrasts with approaches such as Kuroki and Lee et al. 88, which attempt estimate the latent confounders directly.
Pearl 66, Louizos al.
Specifically, Miao et al. 5 demonstrate that outcome bridge function-a function of the outcome proxy W and the treatment-whose conditional expectation equals the regression function, can be integrated over the distribution of W recover the average causal effect. Building on this idea, several methods have been proposed estimate the outcome bridge function and the causal effect, utilizing sieve expansions [10]. reproducing kernel Hilbert spaces (RKHS) [1]
2
W
neural networks 133 14 and minimax learning [S] 16
Y
In particular, Mastouri et [T] proposed two kernel- based methods for efficiently estimating the outcome bridge function: a two-stage regression approach called kernel proxy variable (KPV) and a one-step estimator based on maximum moment restrictions [17] called proximal maximum moment restriction (PMMR). Although outcome bridge function-barn based methods have received widespread attention, alternative line of work focuses on identifying causal effects using ideas inspired by inverse propensity score (IPS) models 19 20 The seminal work by Cui et al. [21) introduced a complementary identification strategy based on a treatment bridge function-a function of the treatment proxy Z and the treatment. This approach, however, is limited to binary treatments and relies on an indicator function over the treatment of interest in the identification formula. Moreover, it requires density
Figure I: An illustrative causal graph for proxy causal learning (PCL), consistent with Assumption (2.2) (5] Observed vari- ables are shown as yellow nodes: de- notes the treatment, Y denotes the out- come, denotes the treatment proxy, and W denotes the outcome proxy. The un- observed confounder U depicted as a white node. Dotted bi-directional arrows indicate potential bidirectional causality (ambiguous directionality) or the existence of a shared latent ancestor between vari- ables. For other causal models that satisfy this assumption, see Table A. 1 in 9
18
ratio estimation to recover the treatment bridge function. Extensions to continuous treatments have
been proposed by Wu et al. 22 and Deaner [10] In particular, Wu et al. [22 generalized the binary treatment framework of Cui al. 21 by replacing the indicator function with kernel smoothing. However, density ratio estimation remains a significant challenge, particularly in high-dimensional settings. To overcome the limitations posed by indicator functions, kernel smoothing, and explicit density ratio estimation, Bozkurt et al. [23) introduced density ratio-free identification strategy based on a novel formulation of the treatment bridge function called kernel alternative proxy (KAP). Their method simplifies the least-squares objective to bypass the explicit density ratio estimation, following ideas similar to those in Kanamori et al. (24), and also avoids both kernel smoothing and indicator function. As result, it is more scalable for settings with continuous and high- dinen-iimeniional treatments.
In practice, the treatment- and outcome- bridge approaches each have different and complementary strengths. Bozkurt et al. [23] illustrate in experiments that the relative performance of each method depends on the informativeness of the proxies [see Section (13) in 233. Specifically, when the outcome proxy Wi is more informative about the unobserved confounder U than the treatment proxy Z, their treatment bridge-based method outperforms outcome bridge-based alternatives. Conversely, when Z is more informative about U than W, outcome bridge-based methods such as those of Mastouri et al. I and Singh [12 yield better performance. However, in practice, it is often difficult to know in advance which scenario applics. This motivates the development of a doubly robust (DR) approach for the PCL setting that combines the strengths of both classes of methods in a way that the resulting approach will remain robust, even if one of the methods fails to identify the causal effect.
DR methods are well established in fully observed causal settings due to their appealing property of consistency if either the outcome model or the propensity score model is correctly specified [2] 26
Such estimators have been generalized to a wide range of functional estimation tasks [27 28 and
are often grounded in the theory of efficient influence functions 22 30 31. These techniques allow
for bias reduction and valid inference even in high-dimensional or nonparametric settings
In the PCL setting, the first doubly robust (DR) approach was introduced by Cui et al. 21 who combined their novel treatment bridge function-based method with the earlier outcome bridge

---

function-bass based method of Miao et al. [5]. Their method is built on the theory of efficient influence functions and achieves the semiparametric efficiency bound. However, as noted earlier, is limited to the binary treatment case. An extension to continuous treatments was proposed by Wu et al. [22]. who replaced the indicator function over the treatment of interest with a kernel smoothing technique, resulting in a nonparametric DR estimation procedure. Nonetheless, learning the treatment bridge component of their DR estimator still requires explicit density ratio estimation. Additionally, the use of kernel smoothing can introduce further difficulties when the treatment is high-dimensional.
In this work, we develop two novel doubly robust estimators for PCL. Specifically: (i) our first DR estimator integrates the two-stage regression method KPV of Mastouri et al. [T]] with the treatment bridge function formulation KAP of [23). and (ii) our second DR estimator leverages the maximum moment restriction strategy KAP. Both estimators retain consistency as long as either the outcome bridge function or the treatment bridge function is correctly specified, thus achieving the doubly robust property.
algorithm PMMR in Mastouri et al. [I]] with the same treatment bridge
Our key contributions are: (i) We propose two novel doubly robust algorithms for causal effeet estimation in the PCL setting that circumvent the need of explicit density ratio estimation; (ii) We leverage conditional mean embeddings (CMEs) derive simple, closed-form estimators using matrix and vector operations-SSalabll to continuous and high-dimensional treatments; (iii) We establish uniform consistency of our proposed estimators, which is stronger than the pointwise convergence typical in doubly robust causal learning; and (iv) We empirically demonstrate that our methods outperform existing PCL baselines in challenging scenarios.
The remainder of the paper is organized as follows:
Section (2) introduces the problem setup and our doubly robust identification result. Section (3) presents nonparametric estimation algorithms. Sec- tion (4) outlines the consistency results for our proposed methods. Numerical experiments are reported in Section 5 and we conclude in Section (6) Our implementation code is available on GitHut"]
2 Problem setting and identification
We consider the problem of estimating causal effects, defined as counterfactual outcomes under hypothetical interventions. Let A E A denote the treatment and Y E e C R the observed outcome. An unobserved confounder U E L affects both A and Y. Our goal is estimate the dose-response curve, as formalized in the following definition.
The dose-response is defined as OΛTE(a) = EY A a, U]), representing the counterfactual mean outcome if the entire population received treatment a. The subscript ATE signifies that its semiparametric counterpart the average treatment effect (ATE).
Definition 2.1.
The key difficulty in estimating dose-response stems from the fact that the confounding variable U is unobserved in most of the applications. To address this issue, proxy causal learning setup assumes the availability of two proxy variables: Z, which serves as a proxy for the treatment, and W, which acts proxy for the outcome. The underlying causal graph is depicted in Figure ( This graphical model admits the following conditional independence assumptions which we leverage throughout:
Assumption 2.2. (Conditional Independencies). The following conditional independence statements are implied by the causal graph illustrated in Figure (7): i-) 1 ZUU, A (Conditional Independence for Y), ii- W 1 ZIU, A and W 1 LU (Conditional Independence for W).
In the following subsections, we review two main approaches to estimate the dose-response curve in PCL setting: (i) outcome bridge function methods (5) [ 12 [3 34], and (ii) treatment bridge
function methods TO, 21 22 2 35 These methods rely on completeness Assumptions (2.3) and
(2.5) that formalize the requirement that the proxy variable be sufficiently informative about the unobserved confounders, highlighting the importance of collecting rich sets of proxy variables in observational studies to mitigate unobserved confounding. While completeness is generally not testable (even when all relevant variables are observed), it is known hold for a wide range of semiparametric and nonparametric models [5) For an extended discussion on completeness within the PCI. framework, we refer the reader to Miao et al. [ Section (S.2)).
3

---

2.1 Outcome bridge canction-basede identification
To identify the dose-response, Miao et al. [S) Wang Miao and Tchetgen 34
show that the
causal effect is nonparametrically identifiable under the model in Figure I given the following completeness assumption [[3] Assumption 2]:
Assumption 2.3. (Completeness Assumption). For any square-integrable functtoo é: U + R. and for any a E A, E[{(U) Z, A = a| = almost surely if and only if E(U) = almost surely.
The completeness condition ensures that the proxy variable Z exhibits sufficient variation relative to the unobserved confounder U, enabling the identification of the treatment effect.
The theorem below identifies (a) via an outcome bridge function.
Theorem 2.4 (Causal Identification with Outcome Bridge Function 34).
Let Assumptions (.. and 2.3 hold. Furthermore, suppose that there exists outcome bridge fincciin by (ow a) satisfying
Y Z, A]
ho(w, A)p(w Z, A)dw.
(1)
Then, the dose-response can be identified by (ATE(a) = E[,,
The proof of Theorem 22.4 can be found in 34 Proposition (3.1)), 13 Corollary (I)).
2.2 Treatment bridge funetion-based identification
Cui et al. [2 and Deaner [0] proposed an alternative identification results using a treatment bridge function, but its learning requires density ratio estimation. Therefore, instead follow Bozkurt et al. [23]. whose method avoids this issue as we describe
Section (3) and forms the basis of doubly robust estimator. The structure of this treatment bridge function-based identification is analogous to inverse propensity weighting principle [19 20]. We adopt the following completeness assumptions identify the dose-response via the treatment bridge function
(23 Assumption 3.3]:
Assumption 2.5. (Completeness Assumption). For any square integrable function E 4 + for all E A, ([((U)) W,A a| 0 p(W) almost surely if and only if E(U) = 0 p(U) a.c.
The completeness condition ensures that the proxy variable W exhibits enough variation in relation to the unobserved confounder U, enabling the identification of treatment effect.
Theorem 2.6 (Causal identification with treatment bridge function (23). Let Assumptions (2.2 and
2.5g. Furthermore, suppose that there exists a treatment bridge function Polz, a) such that
E. (Z, W, A a] p(W)p(a)
Va A,
p(W. a)
then the dose-response can be identified by (TTE(() E (2(, A A= a).
Theorem (.6 is proved in 23 Theorem 3.4]. The existence of the treatment bridge function po relies on Assumption 2.3). along with mild regularity and integrability conditions. In particular, note that Assumption (2.3) underpins the identification of the causal effect via the outcome bridge function. On the other hand, Assumption (2.5) ensures the existence of hg and identification via the treatment bridge Pp- Thus, these completeness assumptions serve distinct yet complementary roles. See Supplementary Material (S.M.) DD) for further discussion on the existence of the bridge functions.
2.3 Doubly robust proximal identification
With two distinct identification approaches hand, now present doubly robust identification result that combines outcome bridge function and treatment bridge function-based identification. Theorem 2.7 (Doubly robust causal identification). The dose-response can be identified with
((RR (a; h, (p) ((h=ha El Z, -Y boo(W,) a)) | A = a) + E[ho a)),
(3)
where he and Po the outcome and treatment bridge functions that satisfy Equations | and 2. respectively. Furthermore, ((DR) (a; h up) admits double robustness such that (DDR) identifies the dose respones curve if either h solves Equation (1 ip solves Equation (2) but not necessarily both.

---

We prove Theorem 2.7 in S.M. in S.M. B
A and relate our estimator to semiparametric efficiency theory
Remark 2.8. In S.M. B we provide a background in semiparametric efficiency theory and efficient inftuence functions (EIF) 129 30 28), and we derive the EIF of the dose-response in the discrete trealment case. Under the mild Assumption (B3). Theorem (B44) shows that the EIF is WATTE = =po(Z, A)(Y - ho,,, ((),,))
p(A] 1(A A))
hhh Woo (TTE((), where
1 ag = agj otherwise.
O ZY, A), and 1[:] is the indicator function defined 1[a a,l
0
Since the expectation of the influence function is zero, we note that:
= )=Eaaar(OO;a] = Elpo a) ho( a)) A = a] B(h(W, a)] AATE(A).
which implies the identification 0 (aT = Eipo(Z, a) - hooW,, a)) A = + EW a)}-
While this ElF-based identification is derived for discrete treatments, Theorem 2.7 shows the same form holds for continuous treatments and preserves double robustness.
Remark 2.9. Our identification result in Theorem (2.7) differs from the doubly robust identification in Cui et al. (27]. They identify the dose-response as EE [q AY ho ho A) 1) al] E,ho(W,o, a)}, where qo is the treatment bridge function solving the Fredholm integral equation E[go(Z. (Z, a) W A = a] = = 1/p(AA = a W) Va E A, and ba solves Equation Q Unlike our formulation based on Equation (2) their identification uses the indicator function 1 a). Wu et al. 122) approximate this indicator with kernel smoothing in estimation. While both methods use joint expectations over (Y, Z, W. A). ours relies on conditional expectation under the distribution p(Y, Z, W. I A=a), This structural difference enables applicability of our methods to high-dimensional trealments, where the seereel-mmotthiggbbee -bsed DR approach fails.
3 Nonparametric kernel methods
We now present a nonparametric estimation algorithm for approximating the function Equation 3 using RKHS theory and conditional mean embeddings. The next subsection reviews the necessary RKHS background, followed by our proposed algorithms' derivations.
3.1 Reproducing kernel Hilbert spaces
In this section, we briefly review reproducing kernel Hilbert spaces (RKHSs), as concepts such kernels and mean embeddings are used throughout our algorithm derivations and consistency analyses. For each space F E {A, W, 2}, we denote the associated positive semi-definite kernel by k>(+-): F F - R. which induces the RKHS H. {e: F R). We denote the corresponding canonical feature map by o-(f) k- f) E Hl. The inner product and norm in the RKHS H,y are denoted by HHz and respectielll The tensor product space is denoted by & and for convenience, we use the shorthand notation H.g. This space is isometrically isomorphic to the Hilbert space of Willerr-SSchidt operators from Hgg to Hr [36), denoted by S2 Hr)) Analogously, we denote the tensor product feature map df( @ dg(g) H,r & Hgg by @Fo(f.g) Throughout the paper, we impose the following assumption on the domains and the corresponding kernels:
Assumption 3.1. We assume that (i) each F (A, W, 2} is a Polish space; (ii) kz(s,-), is continuous and bounded by K, i.e., supjeF [k I K, for almost every f e
For given distribution p(f) on F and kernel ky such that E[k F, F)] < 0o, the kernel mean embedding of F(F) is defined as pr = S (fr(:))((f conditional distribution p(Flg) for each g E 9. the conditional mean embedding (CME) of (Flg)
E Hr 337 38 Furthermore, for a given
is defined as the operator fr (g)) f kr(,, ffggg Hr 39
3.2 Doubly robust algorithm for dose-response curve estimation
We note that the doubly robust function 0 ((DR) consists of three components:
0 (DR) (a) E|y (Z,, | A a)
a) | A a) B(ho(W,)))
Let D= {yis Zi, W;, a; J= be i.i.d. samples from the distribution p(Y, Z, W, A) that are used to estimate 0 (DR) We develop two doubly robust algorithms for estimating the dose-response curve: (i) ATE
5

---

combining the two-stage regression method KPV from Mastouri et al. [1 with the density ratio- free treatment bridge approach KAP from Bozkurt et al. [23], and (ii) combining the maximum moment restriction method PMMR from Mastouri et al. [1]] with the same KAP strategy. Following Mastouri et al. [1 and Bozkurt al. 23 we assume hg E Hw @ HAA and yo E Hz 8 H.A.
Two RKHS-based methods for estimating ho are Kernel Proxy Variable (KPV) and Proxy Maximum Moment Restriction (PMMR), both proposed by Mastouri et al. [ KPV estimates E[h(Waa) in three steps: (i) Given samples {w;; (,,a)))m C D. it estimates Mw (2, a) [[) Z 2, A = a) via regularized least squares. (ii) With second-stage data {U:; 2,, 0=)== C D, it optimizes the sample-based counterpart of the loss CKPV (h) E (Y Eih(W, A) A1)
+ >n.2ll/bll using the representer theorem 44 and the first-stage estimate. (iii) The dose-response at the treatment a is then estimated via the sample mean 1 Σm fuww; a) where { c D can be considered as the third-stage samples, and h de- notes the bridge function estimate. Here ruy ma, and tn denotes the number of samples in each stage. One may either split the training samples across stages reuse them. In the original KPV implemen- tation, data is split for stages one and two, and the entire dataset is used in stage three, i.e., tu = [T]..
PMMR instead uses a one-step estimation via the empirical version of the loss LEMMR (h) g [ (WW, (Z) A)]2 + AMMR|//t|/a closed form using the representer theorem. The entire dataset t samples is used directly for training. Dose-response estimation mirrors the same procedure as in KPV. Additional details and pscudo-code are provided in S.M.
with g E "ZA, which can be solved in
(C.I) and Algorithms 2) and (3).
Bozkurt et al. 23 propose a two-stage regression algorithm in RKHSs, called the Kernel Alternative Proxy (KAP), to estimate the treatment bridge function Pu. followed by third-stage regression to approximate the dose-response curve via @[ a) A = a]. where & approximates Po- A key advantage of this approach is that it avoids explicit density ratio estimation, in contrast to the methods proposed by Wu et al. [22 and Cui et al. [21] Specifically, KAP simplifies the regularized least squares objective as:
(
27 (P())(()) Elp(Z, W. .1.11 | P(W,A)
(K)P(()) = =E
+ 12.011l1l
p(W)p(A) ()) W, W. + plpplll2 + const.
= E [[ A) W. + E
p(W, AA)
= B 4) W, A]]] + EA [E(s (2 A) W, A]] + A2 11 1 + const.
(4)
Here, Eal denotes the decoupled expectation under p(W) ((A)) "const." refers to terms in- dependent of i., and 22 is a regularization parameter. Notably. the objective in Equation 4 involves no density ratio terms. The KAP method proceeds in three stages: (i) Given i.i.d. sam- C D. the first-stage regression estimates the conditional mean embedding (w (2, a) = (W) Z = 2, A = a) via regularized least squares. (ii) Using this estimate
ples (w; 2, a, 1 ng
to approximate the conditional mean E[ A) | W. = (pp. uw aw a) @ da(a)) the second-stage regression minimizes the empirical counterpart of Equation (4) using second-stage samples (w,, āz Th 3 C D. A closed-form solution is derived via the representer theorem [44). (iii) With the estimate S the dose-response is estimated via kernel ridge regression to approximate E[Y a) A a] using third-stage samples (in,, å., d C D. Here, ne, myp and bep denote the sample sizes for the first, second, and third stages, respectively, which may be obtained by splitting or reusing the original training samples. In the published implementation, data is split for the first and second stages, while all t samples are used in the third stage, i.e., tpp = t. Further details are provided in S.M. (C.2). and the algorithm is outlined in Algorithm
As a final step to estimate ((DR) (a), we develop a procedure for estimating El (Z, a (W, A a). completing the doubly robust algorithm when combined with estimators for E|Y (o( a) A ==a aj and E[ho( W. a)). Using the properties of the tensor product, we notice that
E Elspo (2, a)oo A= al ~ ()(Z,))((W,)) | al (h, 3 plll)) HWOHA IA=d a
= E (0.0)(2) 8 (a(()) HHOHA

---

=E { o h (02(Z) 8 )a(a)) & (ow(W) & (Λ(a))) NI0H, |A= (5)
a
As derived in S.M. (C. the expectation in Equation (5) can be approximated using a combina- tion of kernel ridge regression and the properties of the inner products in RKHSs. In essence, the approximation procedure requires learning the conditional mean embedding (2,wIA(a) = [oz(Z) & dw(W)] A = a]. We approximate this term with vector-valued kernel ridge regression. In particular, under the regularity condition (Z W) A = E H. for all g E Hzw, there exists a Hilbert-Schmide operator Cz,wJA E S2(HA, Hzw) such that E [l (Z) & @w((W) A = a] = C2.WIA@u(a)) This operator is learned by minimizing the regularized least-squares function:
EDR( = Σ (z)( & @w(w,; - Cál(aa))iee + >DDRR *
where {2,, W/, amli D denote observations from training set, is the number of given training samples, and ADR is the regularization parameter for the vector-valued kernel ridge regression. The minimizer is given by (zz,ila()) = Cz.WIadA(a) = Σ, (i(a)oz(z)) @ Φw(w:) where (a) (K + (ADR!) Aa KAA is the kernel matrix over {a. Aaa is the kernel vector between training points {a:)i=1 and the target treatment Using this learned conditional mean embedding, we show that the approximation of Equation (5) results in the following closed-form:
Ele(Z, I A a] Σ ((a))(2z,
(6)
We summarize the full procedure in the pseudo-code presented in Algorithm TT. We name our methods Doubly Robust Kernel Proxy Variable (DRKPV) and Doubly Robust Proxy Maximum Moment Restriction (DRPMMR). Specifically, DRKPV uses the KPV algorithm (summarized in while DRPMMR employs the PMMR method (summarized in Algorithm B).
Algorithm 2).
Algorithm 1 DRKPY DRPMMR Algorithms
Input:
Training samples (yA, Wi, 2s, 4
Parameters: Regularization parameter ApR and the parameters of Algorithms 0 and 44p
Output: Doubly robust dose-response estimation g(DR) (a) for E A.
I: Collect outcome bridge and dosc-response curve estimates, h and (.)), with either Algorithm
KPV)
(3] PMMR).
2: Collect treatment bridge and the dose-response curve estimates, & and 02(-), with Algorithm 4
3: Let ds()) be given by (-(-) ) Ei= ((())(z,, (oo0),)))(((w) >)(w.,)) +), where (-) = [ [K + tDu I) K KAo)+
4: For each
A, return the doubly robust dose-response estimation 0 ((DR) (a) (1(a) (a( (.(a).
Consistency results
we present non-asymptotic uniform consistency guarantees for our proposed doubly robust dose-response curve algorithms, DRKPV and DRPMMR. Theorems (4.1) and (4.2) below are a consequence of our non-asymptotic uniform consistency and demonstrate that our estimators converge to the true causal function.
In the S.M. EE]
Theorem 4.1. Suppose Assumptions (3. E.
and E. hold. Then, for the DRKPV algorithm with given training samples (ve. Wi, 2g. a. we obtain that
A A (0DR)) (a) 0xre(a)] 0 with t + almost surely by reducing the regularizer ApR and the regularizers of KPV and KAP algorithms at appropriate rates.
Theorem 4.2. Suppose Assumptions B.7. E..
16). and (E.17) for DRPMMR algorithm with given training samples (yi, Wi, Ro, a4 ai we obtain that
E.6
E.8
hold. Then, for
sup,cA ( (a) OATe(a)) -> with > oo almost surely by reducing the regularizer ApR and the regularizers of PMMR and KAP algorithms at appropriate rates.
7

---

The precise high probability inite-sample bounds on the error in the supremum norm and the corresponding optimal regularization parameters are given in Theorems (E.27) and (E.28). Note that both KPV and KAP involve multiple regression stages, each potentially using different sample sizes (i.e., U4. mi: Ii Dug mig, hip). These samples are derived from the original training set {yi, wi 2i, a,) ) i=1- either through data splitting or reuse across stages. In particular, KPV uses Ra-- and my, samples from the training set for its first- and second-stages, respectively, and t, samples from training set to estimate the dose response. Similarly, KAP uses Rup, 1ig,, and by samples from the training set for its first-, second-, and third-stages, respectively. We detail the data-splitting procedure in S.M. F.2 Consequently, the convergence rate in Theorem (E.27) depends on the sizes of these stage-specific subsets as well as the total training size t. In contrast, PMMR uses the full dataset of size t, and the convergence rate for DRPMMR Theorem (E.28) depends on t and the sample sizes used in the KAP stages. A comprehensive summary of the consistency results for KPV, PMMR, and KAP appears in S.M. (E) and complete proofs for methods' convergence provided in S.M. (F.3)
Remark 4.3. In Theorems (4.1 and (4.2). we establish uniform consistency of DRKPV and DRPMMR, which ensures control of estimation error across the entire trealment domain. These results hold under smoothness and effective RKHS dimension assumptions.
We note that while kernel ridge regression provides optimal rates in Coboley-Matén classes 145), the rates of convergence are slower for larger input dimensions (45) 46). Our results therefore apply in high- digensiinnni settings under dümension-deeenden smoothness assumptions, but addressing increasing- dimnension regimes is left for future work. For further details, see Remark (E.29).
Furthermore, while our estimators are consistent and derived from the EIF, they are not classical one-step estimators 747 (8) and do not automatically achieve local efficiency or asymptotic normal- ity [28]. This remains an interesting topic future work and is discussed further in Remark (E.30).
5 Numerical experiments
In this section, we assess the performance of our proposed estimators for dose-response curve estimation using both synthetic and real-world datasets. We benchmark our methods against several recent state-of-the-art PCL algorithms, including Proximal Kernel Doubly Robust (PKDR) [22], Kernel Negative Control (KNC)
12 Kernel Proxy Variable (KPV) [11, Proximal Maximum
Moment Restriction (PMMR) [.. and Kernel Alternative Proxy (KAP)
23 Except for experiments
exp( exp(-113 (11f f, 11/ (2(2) for each F E (W, Z, A). where denotes the kernel bandwidth. The bandwidth is selected using the median heuristic based on pairwise distances. For PKDR, we follow the original implementation by Wu et al. [22] and use the Epanechnikov kernel. We determine the regularization parameter ADR by utilizing the closed-form expression for leave-one-out
involving PKDR, we use Gaussian kernel of the form ky (fi; ff)
validaaiii (LOOCV) in kernel ridge regression. Either LOOCV or held-out validation set are applied for the regularization terms in the treatment and outcome bridge methods, in line with prior approaches in [11] 12 23). We
provide additional experimental details-including ablation studies on hyperparameter selection, and scalability analysis with Nyström approximation- in the S.M. (F).
Synthetic Low Dimensional: We adopt the synthetic data generation process from Wu et al. 22 which simulates a confounded, nonlinear, and noisy treatment-outcome relationship:
U1 ~ 21 1, 2), U2 ~ (1 10 < U, S 1), W = [U2 + 4111,11), + +N(0, 1)]
Z = UU2 + (N +2[-1, I]], A U, + 10,
Y = cos(2 + 0.3U1 + 0.2) + 1.5A) + ()).
Here, 2(aa b) denotes the uniform distribution on the interval [a, b). and WN,, o²) is the Gaussian distribution with mean j and variance a². We run experiments using training sets of sizes 500, 1000, and 2000. Figure (2a) presents the mean squared error (MSE) of different PCL benchmark methods, averaged over 30 independent runs. Our proposed methods, DRKPV and DRPMMR, consistently outperform competing algorithms and demonstrate improved accuracy with increasing data.
dSprite: We use the Disentanglement testing Sprite dataset (dSprite) dataset 49 a collection of 64 x 64 grayscale images characterized by latent variables: scale, rotation, posX, and posY. Originally designed for disentanglement studies [50), it was recently adapted by Xu et al. [3 as benchmark for proxy causal learning. In this setup, the treatment is high-dimensional vector obtained by flattening each image and adding Gaussian noise. The target causal function is defined as 0 AATE(A) ((vec(B) -32-3000)// /500, where A E R 1096 and B _61x61 with entries of B given
8

---

dSprite Experiment
Synthetic Low Dimensional Data Experiment
45
Algortenn
Risbut kPY
0.35
Doub Risbut poom
Mean Squared Error
0.25
0.30
Kerne Negathe Contri
Alterative
Mean Squared Error
40
KPV
Kernel Aitenaa
Kernel Negattee Cormel
KPV
0.20
PUNE
30
0.15
25
0.10
0.05
0.00
500
1000 Size
2000
500
1000
2000
Synthetic low-dimensional setting
(h) dSprite datasct
Legalized Abortion Experiment
Grade Retention Experiments
0.010
0.07
Doubly
Mean Squared Error
0.06
Negatise
Pray
Mean Squared Error
0.008
PKD
Kemal Abteraatvv
0.05
Remel Regachve Contral
0.006
RPV
0.04
PN
0.03
0.004
8
0.02
0.002
0.01
0.00
0.000
1500 Data Size
IM
Dataset
(c) Legalized abortion and crime dataset
(d) Grade retention cognitive outcome datasets
Figure Dose-response curve estimation across various datasets and algorithms: DRKPV and DRPMMR (Ours). PKDR (22). KAP (23). KNC [12), KPV [IT) and PMMR (T]] (a) Synthetic low-dimensional setting, (b) dSprite dataset, (c) legalized abortion and crime dataset, and (d) grade retention and cognitive outcome datasets.
by By 132 j//32. The outcome is generated via Y = 12(posY 0.5) + where (O 0.5). The treatment proxy Z E F comprises the latent variables scale, rotation, and posX. The outcome proxy W is another dSprite image sharing the same posY' value as the treatment image, while the other latent factors are fixed to scale 0.8, rotation = 0. and posX = 0.5. We run our evaluations using training set sizes of 500, 1000, and 2000. Figure (2b) reports the MSE results averaged across 30 independent runs. Since the implementation of Wu et al. [22] does not support high-dimensional treatments, we omit PKDR. Our methods outperform others on this high-dimensione benchmark.
Legalized Abortion and Crime:
We evaluate our methods on the Legalized Abortion and Crime
dataset [5T), following the preprocessing and setup from 11 22 23
(22) We use the version of In the causal graph, the the murder rate. The treatment proxy
the dataset available from the GitHub repository Mastouri et al.
treatment the effective abortion rate, and the outcome Z is the generosity of aid to families with dependent children, while the outcome proxy W include beer consumption per capita, the logarithm of the prisoner population per capita, and the presence of a concealed weapons law. Figure (2c) reports the MSE results averaged over 30 runs, where each of the 10 data files² is evaluated using three different random runs (leading to different data-splits for different regression stages in KPV and KAP). Our DRKPV and DRPMMR outperform their non-doubly robust counterparts (KPV and PMMR) and all other baselines.
Grade Retention: We evaluate the effect of grade retention on long-term cognitive development
using data from the ECI.S-K panel study 110 53 following the setup of Mastouri et al. [11 The

---

treatment variable A indicates whether a student was retained in grade, and the outcome variable Y corresponds to cognitive test scores in math and reading measured at age 11. The treatment proxy Z consists of the average scores from 1st/2nd and 3rd/4th grade assessments, while the outcome proxy W includes cognitive and behavioral test scores recorded during kindergarten. Figure (2d) reports the MSE results averaged across 30 realizations (3 independent runs per each of the 10 dataset), and compares our method with other proxy-based approaches. On the IM (math retention) dataset, DRKPV and DRPMMR outperform their non-doubly robust counterparts (KPV and PMMR) and all baselines. On the IR (reading retention) dataset, DRKPV outperforms KPV and all other methods except PMMR, while DRPMMR performs on par with PMMR, both outperforming the rest.
Doubly Robust Estimation in Misspecified Setting: We further evaluate the robustness of our methods when the bridge functions are misspecified. By representer theorem [4]. the bridge functions in KPV, PMMR, and KAP can be expressed as linear combinations of feature maps in their respective
RKHSs: h = Exm1 _jm1 Osz (w;) & A (ā,) for KPV, h = Ei a.dw( @ (a(a.) for PMMR, and e --= DES "Mz W.A (Wi, āy) 6 ( for KAP (see S.M. C.1 and (C.2) To simulate misspecification, we first train the DRKPV (or DRPMMR) in the synthetic low-dimensional setup and then perturb either the outcome (KPV/PMMR) or treatment (KAP) bridge coefficients by
adding Gaussian noise. Figures (3ab and (3b) show DRKPV results averaged over five independent runs with standard deviation bands. In Figure (3a. we illustrate the result when the outcome bridge coefficients are perturbed via ai < Osy + Ezj- where Eij ~ (O0 0.2). Similarly, Figure (bb presents results when the treatment bridge function is misspecified by jittering 7ij "ay Esp with Esj ~ (M000000 0.2). In these plots, the term slack prediction denotes the empirical estimate of
Blø(Z, (fflZ)))ho) a)) a) A = a] in Equation 66. Despite misspecification, DRKPV recovers the true causal function, as the slack term offsets the error. Figures (3c) and (3d) show analogous experiments for DRPMMR. Once again, we perturb one of the learned bridge functions by injecting Gaussian noise 0.2) into its coefficients. As the DRKPV case, DRPMMR continues to accurately recover the true causal effect, demonstrating its robustness to misspecification of bridge functions. For enhanced legibility, a larger version of this figure is provided in S.M. (F.4).
5d in S.M. (F.4) illustrate the same experiments under higher perturbation, where the coefficients are jittered with Eij (O, 0.5). Furthermore, additional robustness evaluations are provided in S.M. (F.3) where we adopt the semi-synthetic setups from (23] Section 13] based on
Figures 55
the JobCorps dataset [54 (55). These experiments vary the informativeness of the proxy variables to
challenge the outcome and treatment bridge completeness assumptions (Assumptions (2.3 and (..
DEPOMT Trestrent Erides Moupeciled
(a)
(c)
(d)
Figure Experimental results in bridge function misspecifications with the synthetic low-dimensional
data: (a, b) DRKPV estimates under outcome and treatment bridge misspecifications, respectively; (c. d) DRPMMR estimates under outcome and treatment bridge misspecifications, respectively.
Discussion and conclusion
We introduce two doubly robust estimators-DDKKD and DRPMMR-in the proxy causal learning framework that address unmeasured confounding without requiring explicit density ratio estimation. Built on kernel mean embeddings, both estimators have closed-form expressions and effectively combine outcome and treatment bridge functions recover the dose-response curve. They are provably consistent and remain robust under misspecification of either bridge. Empirically, our methods outperform recent baselines on challenging benchmarks and scale well to high-dimensional treatments. However, the curse of dimensionality poses a subtle theoretical limitation when the input dimension grows with the sample size; addressing this is an interesting area for future work. A key limitation of our work is the computational cost of kernel methods; while we present an initial step toward scalability using Nyström approximation, exploring stronger scalable methods remains a primary direction for future research.
10

---

Acknowledgments and Disclosure of Funding
Bariscan Bozkurt, Houssam Zenati, Dimitri Meunier, and Arthur Gretton are supported by Gatsby Charitable Foundation. We are grateful for the constructive feedback and insightful discussion provided by the anonymous reviewers NeurlPS 2025.
References
[1] J. Pearl and Robins. Probabilistic evaluation of sequential plans from causal models with hid- den variables. In Uncertainty in Artificial Intelligence: Proceedings of the Eleventh Conference on Artificial Intelligence, pages 444 453, 1995.
[2] Olav Reiersøl. Conftuence analysis by means of instrumental sets of variables. PhD thesis, Almqvist & Wiksell, 1945.
[3] James Stock and Francesco Trebbi. Retrospectives: Who invented instrumental variable regression? Jourmal of Economic Perspectives - J ECON PERSPECT, 17:177-194, 09 2003. doi: 10122570895330037662444441
[4] Whitney K. Newey and James L. Powell. Instrumental variable estimation of nonparametric models. Econometrica, 71(5):1565-1578, 2003. ISSN 00129682, 14680262. URL http:
//www.jstop org/ httg//saale/// 15555112
[5] Wang Miao, Zhi Geng, and Eric Tchetgen Tchetgen. Identifying causal effects with proxy variables of unmeasured confounder. Biometrika, 105(4):987-99)) 2018.
[6] Manabu Kuroki and Judea Pearl. Measurement bias and effect restoration in causal inference. Biometrika, 101(2):423-437, 2014.
[7] Christos Louizos, Uri Shalit, Joris M. Mooij, David Sontag, Richard S. Zemel, and Max Welling, Causal effect inference with deep vatiil-v models. In Advances in Neural Information Processing Systems, volume 31, 2017.
[8] Changhee Lee, Nicholas Mastronarde, and Mihaela van der Schaar. Estimation of individual treatment effect in latent confounder models via adversarial learning. In Advances in Neural Information Processing Systems, volume 32, 2018.
[9] Eric J. Tchetgen Tchetgen, Andrew Ying, Yifan Cui, Xu Shi, and Wang Miao. An Introduction to Proximal Causal Inference. Statistical Science, 39(3):375 390, 2024. doi: 10.1214/23-STS111 https:://///t org/10../110 1214/23-SSS1911
[10] Ben Deaner. Proxy controls and panel data, 2023. URL 00283
[11] Afsaneh Mastouri, Yuchen Zhu, Limor Gultchin, Anna Korba, Ricardo Silva, Matt J. Kusner, Arthur Gretton, and Krikamol Muandet. Proximal causal learning with kernels: Two-stage estimation and moment restriction. In International Conference on Machine Learning, 2021.
[12] Rahul Singh. Kernel methods for unobserved confoundiig:: Negative controls, proxies, and instruments, 2023. BBL.itttpss//rrxvvorrrr/aa//
10315
[13] Liyuan Xu, Heishiro Kanagawa, and Arthur Gretton. Deep proxy causal learning and its appli- cation to confounded bandit policy evaluation. In Advances in Neural Information Processing Systems, 2021. URL http::s///pperreeiirr. htt.:///eeerooeoooooooooooooaoaa
[14] Benjamin Kompa, David Bellamy, Tom Kolokotrones, Andrew Beam, et al. Deep learning meth- ods for proximal inference via maximum moment restriction. Advances Neural Information Processing Systems, 2022.
[15] Nathan Kallus, Xiaojie Mao, and Masatoshi Uehara. Causal inference under unmeasured confounding with negative controls: A minimax learning approach, 2021.
11

---

[16] Amiremad Ghassami, Andrew Ying. Ilya Shpitser, and Eric Tchetgen Tcheigen. Minimax kernel machine learning for a class of doubly robust functionals with application proximal causal inference. In Gustau Camps-Valls, Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, Proceedings of 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 7210-7233 PMLR, 28-30 Mar 2022. URL
mlr ppwps/eesseepepp hht.:////ggaaaai2222 htmm
[17] Krikamol Muandet, Wittawat Jitkrittum, and Jonas M. Kübler. Kernel conditional moment test via maximum moment restriction. In Conference on Uncertainty in Artificial Intelligence, 2020.
htttp::///pppp. www.entiiscchol.eeoo
Swagg/CorpusI:::00 211252264
[18) Daniel G Horvitz and Donovan J Thhompson. A generalization of sampling without replacement from a finite universe. Journal of the American statistical Association, 47(260))6633666 -68666666666666- 1952.
[19] P. Rosenbaum and D. Rubin. The central role of the propensity score in observational studies for causal effects. Biometrika, 70:41-55, 1983.
[20] J.M. Robins, M.A. Hernan, and B Brumback. Marginal structural models and causal inference in epidemiology. Epidemiology, 1105)550-550, 2000.
[21] Yifan Cui, Hongming Pu, Xu Shi, Wang Miao, and Eric Tchetgen Tchetgen. Semiparametric proximal causal inference. Journal of the American Statistical Association, 119(546):1348 1359. 2024. doi: 10. 1w..1000016214992222222222 01621459 2023. 2191817
URL https://doi. org/
[22] Yong Wu, Yanwei Fu, Shouyan Wang, and Xinwei Sun. Doubly robust proximal causal learning for continuous treatments. In International Conference on Learning Representations, 2024. URL wtt.httpss//ppprrvviieee htt.://wweeee/oormmmmaddddcCcc..
[23] Bariscan Bozkurt, Ben Deaner, Dimitri Meunier, Liyuan Xu, and Arthur Gretton. Density ratio-based proxy causal learning without density ratios. In The 28th International Conference on Artificial Intelligence and Statistics, 2025.
[24] Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance estimation. J. Mach. Learn. Res., 10:1391-1445, December 2009. ISSN 1532-4435.
[25] James M. Robins and Andrea Rotnitzky. Comments on "inference for semiparametric models: some questions and an answer Statistica Sinica, 1 1(4):920-936, 2001.
[26] H. Bang and J. M. Robins. Doubly robust estimation in missing data and causal inference models. Biometrics, 61(:))962 97099629933900 2005.
[27] Oliver Hines, Oliver Dukes, Karla Diaz- and Stijn Vansteelandt. Demystifying statistical learning based on efficient influence functions. The American Statistician, 76(3):292-304, 2022.
[28] Edward H Kennedy. Semiparametric doubly robust targeted double machine learning: review. Handbook of Statistical Methods for Precision Medicine, pages 207-236, 2024.
[29] Peter Bickel, Chris AJ Klaassen, Peter Bickel, Ya' acov Ritov, J Klaassen, Jon A Wellner, and YA' Acov Ritov. Efficient and adaptive estimation for semiparametric models, volume 4. Johns Hopkins University Press Baltimore, 1993.
[30] Anastasios A Tsiatis. Semiparametric theory and missing data, volume 4. Springer, 2006.
[31] Aaron Fisher and Edward H Kennedy. Visually communicating and teaching intuition for influence functions. The American Statistician, 12- 2021.
[32] Mark J. van der Laan and Daniel Rubin. Targeted maximum likelihood learning. The Interna- tional Journal of Biostatistics, 2(1), 2006.
[33] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), 2018.
12

---

[34] Yilin Li Wang Miao, Xu Shi and Eric J. Tchetgen Tchetgen. A confounding bridge approach for double negative control inference on causal effects. Statistical Theory and Related Fields, 8(4): 262-273, 2024. doi: 10. N0.100002775522992222229220 24754269 2024. 2390748
URL. http:::tp/// org 10. 1080/
[35] Bingxi Zhang, Tao Shen, and Yifan Cui. Neural estimation of treatment bridge functions for proximal causal inference. Statistical Analysis and Data Mining: An ASA Data Science Journal, 18(5):e70045, 2025. doi: https://doi.org// http://oo...rrr0..000///....mmmm
aws/aaass.100 10. 1002/000000 sam.0mmaa...... 70045
[36] Jean-Pierre Aubin. Applied functional analysis. John Wiley & Sons, 2011.
[37] Alex Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf. A hilbert space embedding for distributions. In Marcus Hutter, Rocco A. Servedio, and Eiji Takimoto, editors, Algorithmic Learning Theory, pages 13-31, Berlin, Heidelberg. 2007. Springer Berlin Heidelberg. ISBN 978 978-3-540-75225-7.
[38] Arthur Gretton. Introduction to RKHS, and some simple kernel algorithms. Advanced Topics Machine Learning lecture, University College London, 2013.
[39] Le Song, Jonathan Huang, Alex Smola, and Kenji Fukumizu. Hilbert space embeddings of conditional distributions with applications to dynamical systems. In International Conference on Machine Learning. 2009.
[40] Steffen Grünewälder, Guy Lever, Luca Baldassarre, Sam Patterson, Arthur Gretton, and Mas- similano Pontil. Conditional mean embeddings as regressors. In International Conference on Machine Learning, 2012.
[41] Junhyung Park and Krikamol Muandet. A measure- measure-theorreic approach kernel conditional mean embeddings. Advances in Neural Information Processing Systems, 2020.
[42] Ilja Klebanov, Ingmar Schuster, and Timothy John Sullivan. rigorous theory of conditional mean embeddings. SIAM Journal on Mathematics of Data Science, 2(3):583-606 2020.
[43] Zhu Li, Dimitri Meunier, Mattes Mollenhauer, and Arthur Gretton. Optimal rates for regularized conditional mean embedding learning. Advances in Neural Information Processing Systems, 2022.
[44] Bernhard Schölkopf, Ralf Herbrich, and Alex J. Smola. A generalized representer theorem. In David Helmbold and Bob Williamson, editors, Computational Learning Theory, pages 416-426, Berlin, Heidelberg, 2001. Springer Berlin Heidelberg. ISBN 978-3-540-44581-4
[45] Zhu Li, Dimitri Meunier, Mattes Mollenhauer, and Arthur Gretton. Towards optimal sobolev norm rates for the vector-valued regularized least-squares algorithm. Journal of Machine Learning Research, 25(181):1-511 2024.
[46] Konstantin Donhauser, Mingqi Wu, and Fanny Yang. How rotational invariance of common kernels prevents generalization in high dimensions. In Marina Meila and Tong Zhang edi- tors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 2804-2814 14. PMLR, 18-24 Jul 2021. URL htttps://proeeedingg mlr press/v139/dee dww.//seraaar1....mm .ttml
[47] Whitney K Newey and Daniel McFadden. Large sample estimation and hypothesis testing. Handbook of econometrics, 4:2111-22453 1994.
[48] Johann Pfanzagl and Wolfgang Wefelmeyer. Contributions to a general asymptotic statistical theory. Statistics & Risk Modeling, 3(3-4):37--388, 1985.
[49] Loic Matthey, Irina Higgins, Demis Hassabis, and Alexander Lerchner. dsprites: Disentangle- 2017.
ment testing sprites dataset.
[50] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Mathew Botvinick, Shakir Mohamed, and Alexander Lerchner. beta- VAF: Learning basic visual concepts with constrained variational framework, In International Conference on Learning Representations, wtt.://wwww...ttrooooorrrsssyaaa
2017. URL
13

---

[51] III Donohue, John J. and Steven D. Levitt. The Impact of Legalized Abortion on Crime* The Quarterly Journal of Economics, 116(2):379 420, 05 2001. ISSN 0033-5533, doi: 10. 0..166/003355301511 1140001111144000 144050. URL
[52] Spencer Woody, Carlos M. Carvalho, P. Richard Hahn, and Jared S. Murray. Estimating heterogeneous effects of continuous exposures using bayesian tree ensembles: revisiting the impact of abortion rates on crime, 2020. URL lattts://aar/i org/ abs/2007.... 2007. 09845
[53] Jane Cooley Fruehwirth, Salvador Navarro, and Yuya Takahashi. How the Timing of Grade Retention Affects Outcomes: Identification and Estimation of Time- Mire-ngyying Treatment Effects. Journal of Labor Economics, 34(4):979-102) 1021, 2016. doi: 10. 10.1086/68626.. URL https
//ideas rwpec. org/ http////aac//e//o/iiaaaaaaaaaa/a
(2166668622.225 html
[54] Peter Schochet, John Burghardt, and Sheena McConnell. Does job corps work? impact findings from the national job corps study. American Economic Review, 98(5):184486, December 2008. doi: 10. 18..1255e....5.11555... 100.0.125//aeee 98 5. 1864
URL/ttts:://wwt aeaweb org/ articles?
[55] Carlos A. Flores, Alfonso MoreseLaggunes, Arturo Gonzalez, and Todd C. Neumann. Estimating the effects of length of exposure to instruction in training program: The case job corps. The Review of Economics and Statistics, 94(1):153-171, 02 2012. ISSN 0034-6535, doi: LALL1///RETT__a0000
htttss////oo org 11w.:.1111////1_11111 __017777
[56] Mises. On the asymptotic distribution of differentiable statistical functions. The annals of mathematical statistics, 18(3):309-348, 1947.
[57] Hidehiko Ichimura and Whitney K Newey. The influence function of semiparametric estimators. Quantitative Economics, 13(1):29-61) 2022.
[58] Charles Stein et al. Efficient nonparametric testing and estimation. In Proceedings of the third Berkeley symposium on mathematical statistics and probability, volume 1, pages 187-195, 1956.
[59] Edward H. Kennedy, Zongming Ma, Matthew D. McHugh, and Dylan S. Small. Non- parametric Methods for Doubly Robust Estimation of Continuous Treatment Effects. Journal of the Royal Statistical Society Series B: Statistical Methodology, 79(4):1229-1245, 09 2017.
[60] Kyle Colangelo and Ying-Ying Lee. Double debiased machine learning nonparametric inference with continuous treatments. Journal of Business Economic Statisticx, pages 1-26, 2025.
[61] Houssam Zenati, Judith Abécassis, Julie Josse, and Bertrand Thirion. Double debiased machine learning for mediation analysis with continuous treatments. In International Conference on Artificial Intelligence and Statistics, volume 258, pages 4150-4158, 2025.
[62] Liyuan Xu and Arthur Gretton. Kernel single proxy control for deterministic confounding, 2024. URL httt.t:s///:asst org/ abs 2308. 04585
[63] R. Kress. Linear Integral Equations. Applied Mathematical Sciences. Springer New York, 2013, ISBN 9759462 URL htt.::tttttooootoo http:////oook..oggl..... co. -oookss7dd=---BBAAAA
[64] Dimitri Meunier, Zhu Li, Tim Christensen, and Arthur Gretton. Nonparametric instrumental regression via kernel methods is minimax optimal. arXiv preprint arXiv:2411.19653 2024.
[65] Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer Publishing Company, Incorporated, 1st edition, 2008. ISBN 0387772413.
[66] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331-368, 2007.
[67] Simon Fischer and Ingo Steinwart. Soboley norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 21 21 2020.
[68] Bharath Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert Lanckriet. Hilbert space embeddings and metrics on probability measures. Journal of Machine Learning Research, vl1 1517-1561 (2010), 11, 11 2009. doi: 10.48550%arxiv.0007..555
14

---

[69] R Singh, L Xu, and A Gretton. Kemel methods for causal functions: dose, heterogeneous and incremental response curves. Biometrika, 111(2)497-516, 07 2023. ISSN 1464-3510. doi: 10. 10. 1093//boomet/ssa004.. URI. http:::pp////.
//doioor////o/////////////oo htt.httss://di..oo www./../............. 10.
[70] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005. ISBN 026218253X.
[71] Milton Abramowitz. Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables,. Dover Publications, Inc., USA, 1974. ISBN 0486612724.
[72] Bharath K. Sriperumbudur, Kenji Fukumizu, and Gert R.G. Lanckriet. Universality, characteris- tic kernels and rkhs embedding of measures. Journal of Machine Learning Research, 12(70):
-410, 2011. URL http://jmlr. wttp:/rrrgggpappps//2222 http:ppppprrs/////pprrprruuuuuur
sriperumbuuurrrr 1 1a. html
[73] Christopher Williams and Matthias Seeger. Using the nyström method to speed kernel machines. T. Leen, T. Dietterich, and Tresp, editors, Advances Neural Information Pro- cessing Systems, volume 13. MIT Press, 2000. URL/tttpss//prrcceedngg.. neurips. cc/ haper_files/pppeeere /0tpapper2000///Fee httpp///w....teeeeeeocoooeeaeeee
Paper .df
[74] Sanjiv Kumar, Mehryar Mohri, and Ameet Talwalkar. Sampling methods for the nystrom method. Journal of Machine Learning Research, 13(34):911-1006, 2012. URthttp::/mmlr
org/papers/viae htt////ar///////// kumar html
[75] Giacomo Meanti, Luigi Carratino, Ernesto De Vito, and Lorenzo Rosasco. Efficient hyper- parameter tuning for large scale kernel ridge regression. Gustau Camps-Vall Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 6554 6572. PMLR, 28-30 Mar 2022. URL https::/ www.innnns........mmmmmm....mm.m press/
v151/meanti22a.htn
15

---

Supplementary Material for Density Ratio-Free Doubly Robust Proxy Causal Learning
Section (A) reviews our doubly robust identification result and proves Theorem (.7). Section (B]
outlines semiparametric theory and derives the efficient influence function for discrete treatments. In
Section C we present derivations for our proposed algorithms: Doubly Robust Kernel Proxy Vari-
able (DRKPV) and Doubly Robust Proxy Maximum Moment Restriction (DRPMMR)) Section (D)
discusses conditions for the existence of outcome and treatment bridge functions. Section (E) proves
consistency of our methods, Section details the experimental setup, including hyperparameter tuning and additional results.
A Doubly robust identification
Here, we prove Theorem 2.7 For completeness, we restate the theorem below.
Theorem A.1 (Doubly robust causal identification; Replica of Theorem 2.7 The dose-response can be identified with
( (a; h, (P) ((m=hu 1la aY -oo a)) A a) + E(ho(W,a))
where hg and Po the outcome and treatment bridge functions that satisfy Equations (2) and 2 respectively. Furthermore, 0 0/DR) (a; h, up) admits double robustness such that ( identifies the dose-nesponse curve if either h solves Equation (7 ip solves Equation 2 but not necessarily both.
Proof. To understand why (a) exhibits double robustness, we consider the following two cases:
Case i: Suppose that h is correctly specified, i.e. h = hg. Then,
9 (DRR (a; ho, 4p) El (Y ho( a)) al + ((,
= ((2,a) (f(V- ho(w, a) ) w/ d)ddd p(=la)dz ([, (a)
=U due Equation
where the final equality due to Theorem 2.4
Case it: Now, suppose that yo is correctly specified, i.e., = 40. Then, we note that
ATE (a: h po) B[y so (o, a) A a] Elpoo(Z, a) A + B[)(Waa)]
First off, note that the first component EE[ (o(( a) A a] identifies the dose-response due to Theorem 2.6p. Next, consider the remaining two terms
E3 El Z (W, a) A= a] + E(M(W,a))
-- /
f P0(2, a)h(w, wla)dzdu + B(h(W,
/ / (P([, a)h(w, ()
(W,,)]
/ (/ ( a) a)p(2 a)dz h(w, ((w.a)p(wladdu BB/((
= / (E(pu(Z, w) a)) h(w, a)p(wla)dw ((ho
p(w) hww, a))(wla)dw + E[A(W,,)) p(wla)
EEM(W, (BW(W, a)] 0,
where the first equality in the last line above is due to Equation
The following lemma the direct implication of Theorem 2.7 which shows that the dose-response curve also be identified with the conditional expectation of the multiplication of bridge functions.
16

---

Lemma A.2. Let he and Po be outcome and treatment bridge functions satisfying Equation ! and (2) respectively. Then, the dose-response can be identified by the following slack function OATE Rlla (W A a).
Proof. As we observe in the second case of the proof of Theorem 2.7
Elso(2, A A=d = /
p(w)
ho a)p(wla)dw E[[o(WW,)) = (ATE(a).
p(wla)
B Derivation of efficient influence function in discrete treatment setting
B.1 Background efficient influence functions
Let P denote nonparametric statistical model, and let v/2 + be target functional of interest such the dose-response curve (a). In semiparametric efficiency theory, statistical parameter such as the average treatment effect is viewed as a functional mapping the underlying denote by IP the true distribution with Pf = E[ Pa the empirical distribution with IP, (E) := E ((Z)). The C2(P) norm is written ||4|| := f e dIP. We are interested in how the functional & changes as the data tata-geenrating distribution p E P varies. When 0 is sufficiently smooth, it admits Gâteaux derivative, defined
data-generating law p E P to a real number. In this section,
d 0(p + cH) de
0' (p: H)
for H in the tangent space P. The functional is said to Gâteaux differentiable P if 0' H) exists and is linear and continuous in H. The Gâteaux derivative formalizes the local sensitivity of distribution. extends the notion of a directional
functional to infinitesimal perturbations derivative to the space of probability laws and provides linear approximation of how 0(p) varies when p displaced along an admissible direction H within the tangent space.
This gives rise to the von Mises expansion [56] 28
(4)) = (p)) + (p: - p) + ()
where R2(p, q) second-order remainder depending on products or squares of differences between q and The von Mises expansion plays analogous role to first-order Taylor expansion in finite dimensions, expressing the change in the functional as a linear term term- canturd by the influence function-ple a higher-order remainder. This expansion enables the definition of the influence function.
Definition B.1 (Influence function influence function of f at
29
g
A measurable function ((:p) C2(p) is the
(ppp - p) = / (1:pp) d(a-dd((p))) - (z), for all q E P
with f (2::))p((2) = and Ep (Z: < oc.
The influence function 1b is the Riesz representer of the derivative 0' ('(p:)). as La(p) is Hilbert space. It is also referred to as the pathwise derivative gradient [57] or Neyman orthogonal score Intuitively, the influence function measures how sensitive the parameter of interest is to infinitesimal changes in the dlta-generating distribution. It can be seen as the "directional derivative" of the functional with respect perturbations of the underlying distribution like gradient in function space. This makes it a key tool for constructing estimators that correct first-order bias and achieve minimal asymptotic variance.
30 33
To compute influence functions, we analyze parametric submodels {pt] C P such that Pe=o = p and + P. smooth. A canonical choice is the tilted model:
dp dp
(z) = 1 Es(2), where ||s|loo S M, < 1/M
17

---

The associated score function is sz(2) = s(2) = de log pc) 58 31
Applying the Mises expansion to Pc. obtain:
((p.) - + f ((2:))s(=)dpp())) + R2(p,
Since (p, P.) is second-order, it satisfies de R2(p 0, condition known Neyman orthogonality [33 Thus,
d
- =
((pc)
de c-=0
This defines pathwise differentiabilityy and the function & satisfying this identity the efficient influence function.
Definition B.2 (Pathwise differentiability 129 28)). A functional 0 is pathwise differentiable at p if for every smooth parametric submodel {p+} P with score se sz(z), we have:
d
= / (x)
8(pe)
de
c=0
Once di is known, it can be used to construct an efficient estimator. Suppose P estimates p. and let Px. denote the empirical measure. The plug-in estimator 0p := 0()) admits the expansion:
p - (p) = EE 1 (Z: p)] + (p.
This motivates the one-step estimator, which corrects the plug-in estimator with an estimation the first-order bias (1/1//22 (Z;
8 PF.
When analyzing one-step estimators, is common write the following decomposition 28
-(p) (P, (p) + (P, P) (Z:z))
(2
empirical process
second-ords remainder
The first term is sample average and converges by the central limit theorem. The second term, known the empirical process, vanishes at rate oz (1/2) if the estimated influence function & is fa-consistent (i.e (Z) v (Z)]l2 >p 0, where | = f (dp..) which is the case it lies in
Donsker class or is estimated via cross-fitting (see 28
The third term, R2 (f, p), is the second-order remainder, which captures nonlinear errors. This term often exhibits a rate which depends on product of nuisance error rates. Therefore, under standard product-rate conditions, such as nuisance errors converging at rate n 1/4, we obtain R₂ or(n 1/3).
Together, these results imply that 0 asymptotically linear and efficient. Specifically, in nonpara- metric models, the best achievable precision for estimating a target parameter 0 is governed by local asymptotic minimax bound which is attained by this estimator.
This bound extends the classical Cramér-Rao inequality using pathwise differentiability and para- metric submodels. For any smooth submodel {p-} } with score s(2), the Cramér Rao lower bound is:
E 2 (Z: S [[((Z;))
[[s(2)2]
where &b is the influence function. Equality is achieved when s = v. yielding the nonparametric efficiency bound:
var(s(Z:p)}
Since the one-step estimator satisfies
Vn 0) ~ N(0, var w/( (i(Z:p)) ())))),
it achieves this bound and is thus locally efficient.
18

---

Eventually, note that the efficient influence function (EIF) is defined relative to the model's 's tangent space, which formalizes the allowable directions of perturbation around the true distribution. The tangent space is the set of valid perturbations around the true law and thus determines which directions contribute efficiency.
In nouparametric models, the tangent space is the entire Hilbert space of square-integrable mean-zero functions:
T = {e E (F( [=[[()]] = 0}.
This ensures that any regular, pathwise differentiable functional admits unique influence function e that also serves as the EIF.
In contrast, a semiparametric model imposes structural constraints on the distribution (e-g.. conditional independence, smoothness, dimension reduction). These constraints reduce the tangent space to strict subspace (F)). In this case, there may exist multiple functions e satisfying the von Mises expansion, but only those that lie within the tangent space correspond to valid perturbations. Among them, the EIF is defined as the unique d that also lies in the tangent space and minimizes variance. Therefore, careful attention must be paid to ensure that the candidate vé lies in the tangent space, as only then does it correspond to the efficiency bound [28].
B.2 Derivation of the efficient influence function
In this section, we consider discrete treatment setting, i.e., A E 10..1....,44 We define the conditional expectation operator T (2 C2(PZ.A) by T(l) E(e(W,,) Z, A], and its adjoint 7" (P(,,))) - (PWw.)) by (c) (2, W. Al. Pz.4 and Pw.4 denote the joint distributions of the random variables (Z, A) and (W, A), respectively. Here, L2(F. denotes the space of square- integrablee functions on domain F with respect to the measure p. We impose the
following regularity condition on the conditional mean operators following the setup in Cui al. [21::
Assumption B.3. T and T are surjective.
Under this setting, we derive the efficient influence function (EIF) of the dose-response curve, following the framework of [21] Theorem 3.1].
Theorem B.4. (Efficient influence function) The efficient influence function of ATE((a) where Assumption B.3 holds, and Equation ( holds at the true data generating law, is given by
1 VATE a) = =o( A)(Y hoW,, A)) p(A) 1[A= - a) ho a) - OATE(U)
where O is the collection of the variables AY W, 2), and 1[:] the indicator function. Therefore, the corresponding semiparametric local efficiency bound of (TTE(a) equals to E [MATE(O; a))]
Remark B.S. We note that our identification result for the dose-rew rose-response curve in the continuous treatment setting can be viewed as generalization of the identification result implied by the efficient influence function in Theorem B.4). Specifically, observe that, since the expectation of the influence functioo is zero, we have:
p(a) 0 = ((.p((,a))( - ho ho(W, a)) 1[A a) ho (W, a) ATT((a]
=
ho(W, a)) A = a] + B[h(W, a)] - (ATE(a)
Therefore, in the discrete treatment case, the dose-response curve can be identified by
(ATE(a) TElyo (Y - ho,W))) A a) + R(h(W,a)) a)]
In Theorem (2.7 we show that the same representation identifies the dose-response curve in the continuous treatment case, while also exhibiting the double robustness property.
Remark B.6. We note that the formulation of our efficient influence function (EIF) differs from that of Cui et al. (2I) particularly the treatment bridge function term. In their framework, the treatment bridge function Z x A + R is defined as the solution to
P(W)
A= E (z (W
Va A,
(8)
P(A
W)
p(W,a)
19

---

and the outcome bridge function solves 7 Under the Assumption [B.3) and assuming that Equation 77 holds under the true data generating law, Cui et al. (2I) shows that the EIF of the dose-response given by
ATE a0)
ho (W, a] + ho a) - (Ar
which leads to the identification formula
( TE [go A) - ho (W,
a]] E(h(W, a)..
Instead, our identification result in Equation 7 differs from Equation 2 in two key ways: (i) We use a treatment bridge function so that solves Equation ( which includes an additional factor of p(a) in the numerator compared to Equation (8); (ūi As result, the first term in our identification formula involves the conditional expectation E. El Z. aY hW,, a))
a) over the distribution p(Y. Z. W A == 0.), whereas in Equation 99 the expectation is taken over the joint distribution p(Y, Z, W, A) with an indicator function enjorcing the treatment level. Because of this difference, our framework extends naturally to continuous treatments without the kernel smoothing step of Wu et (22), whose technique has been extensively used in the doubly-robust literature to handle continuous treatments [59 60 6ly. Furthermore, this structural difference extends the applicability of our methods to high-dimensional treatments, where kerne- semee-smoothin based DR approaches are ineffective as demonstrated with dSprite experiment in Section
parametric sub-model for nonpara- metric P indexed by which includes the true data generating law at 0 for some mean-zero function É 0 + R,
Proof of Theorem (B.4). We consider the following 28 31 57
(O p(0) cl(O))
(10)
where O = (Y, Z, W, A) is the collection of variables of interest, |1e| S M < co, and < 1/M
so that p. (O) 2 O. We have the score function s=(O) de
# log (O) -= To find the efficient
influence function for 0 (TEE()), we need to find the curve ATE that satisfies
0 _ 9ATE (a:p.)
(1I)
De
First of all, we recall that
E|Y A] / h(w, A)dp(w A)
and define the bridge function h2 so that:
Ep, (Y | [. = / f he(w,, A)pa(w Z, A)dw.
Therefore, we have
I (u neww,, Pe),, a) le-o dwdy O,
which implies that
I @h.(w,o a) de
p(w|z, a)dw= : / (u --h(w, a)) alogpe (w, uız, a) Pe(w (w, 812, a)l
dwdy (12)
de
Now, we recall that
@ATE(a;PE) = / he (w,, ())x
Hence,
a0 AT
(PATE(a:pp
(w,(wa)) a)
p(w)dw + f N(w, a) @p((w) c==
dw
de
c=0
de
de c=0
/
ah, a)
p p())dw |blu, 2 blog (w
pl pww)dw
(13)
Oe
l==0
20

---

We notice that
f h(w, a)
① logg Pe (w) de
log p. (Z, Y, W)
p( = E (h(W, a) - (ATT((a)
(14)
de
One can show this as follows:
Elog Pe (Z, Y, Y, W.
(h(W.a) a)- (ATE(a))
dc
f h(w, a) h(w, a)
Op.(z.B.d',u
/ Plogpe
dp(2.y. y. d', w) - Oaro(a)
dp(z, y. d. w)
de i (l)
de
dp(z,y, a.y.dd.w) w)+ h(w,a) 2 log (.
dp(z 2 y, y, w)
DE (logP((w) de
6=0
de
c==0
/
I PlogPe(Z,Be (w)
-ATTE(a)
dp(z U. d' w) - -xTE(a))
dp(= y, w)
de
h(w, a)
i (w)
dp(z, dp(z,y,d',)) / (((wwa)) (llwwa)-0AAaa) OATx(a))
(lo
dp(z, y, ,
de a Ologp.(w)
==0
de
l==0
f (w,aa) a)
dp(z,v.d', V. ww)
dc
+ / (M(w,, a) - 0rTE(a))
(
log Pe (z, Pe a'
dp(z, v,a',lw) dp(w)
de
= =f =/M h(, a)
d Ologp. dp(w).
de
e=0
Next, we consider the first component in Equation 13 (w.a)
(h,(w,a)
p(w)
p(w)dw=
- o p(w, a)dw
-f @w, a')
1/a a'] p(w)
plw, a' d')dwdda
de (h, a') de
6=0
p(w, a²)
1/a a') p(a') BlPu(Z,a') W w. A= a']p(w,,'') a')dwda'
cm0
(he a')
1/a a') p(a") pol2, a' (l(2, w, a' ddwda'
dE
k=0
/ 1/a a'}
(fo (h.(w,a)
p(aa) polzza')
p(w|z, al )dw plz, a' 'dzdd'
Oe
/ 1/a a') p(a') p polz, al) (ʃ (1w-h(w..))) kw,,a)
Ölogp. 10.y12,00 a) p(w, uız, a' ddwdy) p(z, af ddddd
de
(15)
/ 1/a = a'] p(a') po ('y- h(w, a))
Mogg. (w, v12,0l)
p(w, y, z, a' ddwdydzdd
de
+ / 1/a a') 1 p0 (v Mw,a)) Alogp.(z,
plw, y, 2, ydddydzda
p(a')
Oe
2 log Pe (W,Y,
E 1/A p(A) Pol Po(Z, h(W, A))
(16)
de
where equality in Equation
[5 is due to Equation Now, combining the results in Equations
14h
and 16 obtain
00 AT 10BATTE(O (c;
de
==0
21

---

E ( [|A (1AA=a) -
1
Olog P. (W, Y. Z, A.
p(A) Po A) h(W, A)) + h(W.a) (
de
which shows that the efficient influence function for 0 (a) is given by
VAT=(O;;) - ()(2, h(W 1
p( (A) [ - a) + ho OATE(O)
(17)
Next, show that (a) belongs tangent space Ti + Tz using Assumption
(B.3) where
T = (s(2, A) E (2 x A,PZ,) B[s (2, A) 0}
T2 W A) E (2 AZ, E [(Y hW, A), Z, Z, e cd((T)))).
Here, (2 x Az,z, is the orthogonal complement of C2 x A, A,,z,a), R(T) is the range space of the conditional mean operator T, and cl(.) denotes the closure. We decompose Equation as follows:
B[h(W, a) (TT () Z,A]
+ h(W, a) (A) - (l, a) 0AI=(a) A]
+
1[A = -a
h(W,
p(A)
Notice that ([((W,,) a) OTTE (a)
[, E Ti due to Equation
For the remaining two terms,
notice that
E ST( a) AIL (a) [W, a) (x)=(a) A. Z,A]=
(18)
E
1[4 - ao Z A)(Y h(W, A)) A = 0. P(A)
(19)
Hence, using Equations 18 and 19 and Assumption B.3 we observe that
E [{ ho(W, A)) (h(W,a) a) (a) E[h(W, a) OATE(a) Z, Z, A] c(R(T)),
1
E
P(A) 1[A= all==aoo((2 A)(Y h(W, A))² Z, cdR(T))
C
Algorithm derivations
In this section, we derive the doubly robust PCL algorithms, called doubly robust kernel proxy variable (DRKPV) and doubly robust proxy maximum moment restriction (DRPMMR). We begin by reviewing the results from Mastouri et al. [1] and Bozkurt et al. 233 which establish estimation procedures for the outcome and treatment bridge functions: namely, the Kernel Proxy Variable (KPV), Proxy Maximum Moment Restriction (PMMR), and Kernel Alternative Proxy (KAP) methods. Building these, we introduce fully kernelized doubly robust algorithms for estimating the dose- response curve from observational data, without requiring an explicit density ratio estimation step. This property makes our methods particularly suitable for continuous and high-dimensional treatment settings.
C.1 Outcome bridge function-based approaches
Two kernel-based methods KPV and PMMR have been proposed in ||| to estimate the outcome bridge function h W x A + R, which satisfies
/ ho(w, A)p(w Z, A)dw.
FI 2 =
Once hy satisfies this equation, the dose-response curve can be identified by E Ei (W, a)) 5
as stated in Theorem (2.4). In both approaches proposed by Mastouri et al. I1 ho is assumed to lie the RKHS Hw 8 H.A. Therefore, we note that
0 (ATE(a) EE/ho a)] = E (ho, (w() @ H
22

---

= (ho, Blow (W) @ HWSHA ~ (ho, @ @ HwSHA
where jow- E. Σラ= ow((;;) the sample-based estimate of the mean embedding [W)] As result, by replacing ho with its approximation h in the above expression, the dose-response can be estimated using the inner product:
0 OxTσ(A) ss (hr (ir @ @a(a))
Hw@HA
In the following, we review the estimation procedures of both KPV and PMMR proposed by Mastouri et al. П. Both methods yield closed-form solutions for the estimation of the bridge function and the causal function, which we present in the subsequent subsections.
C.1.1 Kernel proxy variable algorithm
To solve the outcome bridge function equation in Equation [ the KPV method aims to minimize the following regularized population loss function, assuming that h Hwaa:
Ckrv(h) E (Y E[h(W, A) Z, A)) + An.all/i//
where 26, is the regularization constant for the second stage regression of KPV. Since this loss function involves the conditional expectation ([,(W, A) Z, A), cannot be directly minimized. To address this, the KPV algorithm decomposes the problem into two-stage regressions. Note that, under Assumption
3.I). the conditional expectation can be rewritten as
Bl)(W,A) Z, A]
E [(,, ow (W) (A(A))) | Z,
(h, ((ow(W((W) A @ (A(A))
Thus, the first step in KPV is to approximate the conditional mean embedding (w]z,a(z,a) [low(W)) Z 2, A al using vector-valued kernel ridge regression. In particular, under the regularity condition that [Z) = A = lies in Hz, for all g E Hww there exists an operator CWiZ, E S2(Hzz,, Hw) such that μw Z. a) Cwiz, A, (0z @ @a(a)).
Let {w; 2,, ) be i.idd samples from the distribution p(W, Z, A), representing the first-stage regression samples, where nn is the number of first-stage samples for KPV algorithm. The conditional expectation operator is estimated by minimizing the following regularized least-squares objective:
1 Thk
CxPV(( (C)
(ow()))) - C (bz(z,) @ (a:)) || + An. ||C|| (,, (HW.HHa))
where An,1 is the regularization constant for the first stage regression of KPV. The minimizer ČwIZ.A has closed form solution and is given by [1]
Wh
Cwizz z. (
3 &d)))) = (w1z,((2,0) = > (z).)))(
where
B(=,a)= -K C KĀā + (K₃z 2 Kūa)
By substituting this approximation into the sample-based version of the loss function CVPV(h) (h), the KPV algorithm estimates the bridge function ho. Specifically, let (ū:, ž,, ai be i.i.d. samples from the distribution p(Y, Z, A), denoting the second-stage data, where my is the number of second-stage samples for the KPV algorithm. The sample-bass based KPV objective then given by
mi
EKPV (1)
mi Σ
(Us - (h, (Z (a.)))) An.llh P
(20)
By the representer theorem
h admits the representation
Σ
anow(Wu,)) ( (a))),
(21)
1=1
23

---

for some coefficients a:j. Substituting this representation into the sample loss in Equation (20) KPV finds a closed-fome solution for the bridge function by minimizing for the coefficients atj- For further details, we refer the reader to
Below, summarize the pseudo-code for the KPV algorithm. Unlike the original implementation in [T we use more numerically stable version that optimizes fewer parameters, as presented in [62 (see Appendix F in (62).. This version of the KPV algorithm optimizes my parameters rather than Bi X mA parameters.
Algorithm 2 Kernel Proxy Variable Algorithm [11
Input: First-stage samples {1x &, 0s Ji Second-stage samples {ūr, & âs 1m
Parameters: Regularization parameters As,1 and 2, and kernel functions kz(:)) for each F (A, 2, W). Output: Bridge function estimation h. I: For variables F' E {W, Z, A} with domain F, define the (required) first and second stage kernel matri- ces/vectors as follows:
K F = |k (fr, Filys ER"AX", = = (F. fillss E E xma,
K FF = [kr(f. @Rmaxmn, KFy = (Fr f)], E RMT,
Ki = [k>(fs ERT
2: Define the following matrices and coefficient vector:
B AK 2 K2z + 1) (KAA & K K2z),
M K KAΛ O (B KWWBB
a (M + madn,1)) Y.
3: The bridge function estimation is given by h(w, a) a (Kio (B Kww)). The dose-response estimation given by
(TT = Σ h a) = Σ a (Kio @ o (B'k..)))
where we take {rit (/m
C.1.2 Proximal maximum moment restriction algorithm
The PMMR algorithm relies on the following conditional moment restriction result:
Lemma C.1 (Lemma (1) in [1). measurable function h on the domain W x A solves
E Y Y Z = A a- f h(w, a)p(w 2, a a)dw
if and only it satisfies the conditional moment restriction (CMR): E[ N(W, A)| A| =
Lemma [C.1) implies that E((Y h)W.))))(( Z, A)) = 0 for all measurable defined on the domain Z A. Thus, solve the CMR. the PMMR algorithm solves the CMR by minimizing the regularized maximum moment restriction (MMR) objective to find the function h E HwA:
LMMR (h) sup E(Y h(W, ((W.)))(ZZ.)) A))2 9EH2A 1
AMMR [|l||Haa-
(22)
Under the integrability condition E [XY - h(W, (oz, A), ÓZA (Z) < co, Mas- touri al. [ Lemma (2)] shows that the regularized MMR objective admits the equivalent form:
LMMR (h) = E [[Y h(W, AY) h(W', A3) (z (Z)) + AMMR || [
(23)
where V denotes an independent copy of E A, 7, W. and Φ is the canonical feature map of H.Z.A. In this paper, adopt the tensor product structure for the feature map and RKHS by defining
24

---

= & @ and Hz.= = Hz @ H.. This choice allows us to construct H...1 from separate kernels over Z and A. However, also possible define a kernel directly on the joint domain 2 A for PMMR algorithm, in which case the feature map and RKHS do not need to take the form of tensor product.
As result, given training samples {vi, Wi, 24, a. ai)i the empirical objective for PMMR is defined
LMMR (1)
(y K5) kig + AMMR NAMR//////wwar
(24)
1=1
where h; h(w,, a:) and kij Hw @ H.s By the representer theorem 44). the solution Equation (24 admits the following representation:
(0za &i, a:), (ZA a,)=,,, Recall that we assume
= Σ a.w((ws) @ da(a;).
Substituting this representation into the sample loss in Equation (2)) the empirical objective becomes
1
á = arg min 2 Y Lox T W(Y La) + AMMR La,,
where Y = |y 1/2 Vel a = [ai 02 as} [| = kww Wj ( ay), and
[W]; kz 23 12 (a,, a,). The PMMR algorithm solves for the coefficient vector We
summarize the procedure in Algorithm (3), following the variant presented by Xu and Gretton 662
Algorithm 3 Proximal Maximum Moment Restriction Algorithm 62
Input: Training samples [w;, Fi, di Ji J==1,
Parameters: Regularization parameter AMMR. and kernel functions (r(,)) :) for each F E {A, W. Output: Bridge function estimation h. I: For variables F {W, Z, A} with domain F, define the (required) kernel matrices/vectors as follows:
KFF [k>(St Rxt Kr, (Kr (f, R
2: Construct the following kernel matrices and the coefficient vector:
Kzz, a a=va c=va(VGLVC L VGY. where /G is the square of the matrix G.
L = KAa @ KwW, G
KAa G=
3: The bridge function estimation is given by h(w, a) a (KAo 2 Kww).
4: The dose-response estimation given by dATE(a) R El= h(w,, a) = Ei a (Kaa @ Kww)).
C.2
Treatment bridge funstion-based algorithm
An alternative identification result in the proxy causal learning setting, based on density ratios, is presented in [23] Unlike previous works 22], this approach bypasses the explicit density ratio
estimation step by leveraging a novel bridge function definition. Specifically, 23 proposes Kernel Alternttve Proxy (KAP) algorithm, two-stage regression method to estimate the treatment bridge function PP satisfying
E (Z, W. A = a) p(W)p(()
p(W,a)
This new approach aims to minimize the following regularized population loss function, assuming that & E HZa
(RAP((> = E (r(W, A) [b((ZA)) W. A))" + >xp.2lvlllzz
25

---

E (o W.
W.A]") 21Ew E (2( W, 222211lllzzz const.
(25)
where 2 is the regularization parameter for the second stage regression, the notation Ew BAll denotes the decoupled expectation taken over W and A, i.e., for any function E, Bu A ( A)] J f(w, ()p(w))p() (p))dddd, and 'const'' represents the terms independent of ş. As in the outcome bridge function estimation, this loss function cannot be directly minimized due to the presence of the conditional expectation term IE (ip(Z, A) A]. Notably, appealing to Assumption 3.1. this expectation can be rewritten as
E ( A)
W. A]
E [(p, 2z (Z) @ QA (A)) W W;A] A] (p,, El (Z) W. (a(AA)) Hz
Thus, the first step in KAP method is estimate the conditional mean embedding MZ (w,, a) (Z)z(Z) W = w, regularity condition that ([g(2) [) A,A=
using vector valued kernel ridge regression. Specifically, under the lies in Mwa for all E Hz. there exists an operator Czww,A (2HW..: 2) such that Mz|w,A (w, a) CzJw,A (#w(w) @ @a(a)) Let [i, 2,, 0g ==1 be i..d samples from the distribution p(W, Z, A), denoting the first-stage samples, where ne is the number of first-stage samples for KAP algorithm. The conditional expectation operator is learned by minimizing the regularized least- squr-s function:
(CPP)
Σ (0)())) - C (dw(m)) & (a()))) + A Ap.: n
1
12
where Xx,1 is the regularization parameter for first-stage regression of KAP method. The minimizer C2/w, admits closed form solution and is given by 23
Az a) = (22)) = Σ B,w a))z(z))
where
3(w, a) (Kww @ Kā + 1 (Kw 2 K Kūa)
By substituting this approximation into the sample-based version of the loss funetion ( Bozkurt et al. (23) derives an algorithm to estimate the bridge function Po- Let { {w, ā. }= be i.i.d. samples from the distribution p( (W,A), A), denoting the second-stage samples, where my is the number of samples for the second-stage regression of KAP algorithm. The sample-based version of the loss function in Equation
25 can be written as
((P))AP-
1
Σ
(p, (ws, ā.) o (a(()))
mip
==1
2 mm 1)
p. (. a.) da(i)) (
+ +1p21lPl1 (26)
Hzd
The representer theorem
44
ensures that yo admits the representation
mp m
6
muzz ( (1,, a.) (a()))
7m mm 1)
Pzw. ( (
for some coefficients mh By substituting this representation into the sample loss in Equation (26). the minimizer & admits a closed-form solution, as given in Bozkurt et al. [23 Furthermore, the dose-response curve can be estimated via (ATE(a) x E.Y A a). Noting that
which outline in Algorithm
@ E[Y = a) = E Y (Z) @ A=
= EYz(() = a) 8 (oll)))
(27)
26

---

evaluating the above inner product requires estimating the conditional mean embedding Hyzja B| (2 A == a). which can be obtained via vector-valued kernel ridge regression, similar to the first-stage regression. In particular, let {úr, ž;, du )= denote i.i.d. samples from the distribution p(Y. Z, A), denoting the third-stage samples for KAP algorithm which we take to be {ūe- &2, de ) (īi, 2i, Jiji=1 = (Mi: =,, j- Then, the conditional mean embedding μy 2)A(a) CYZlAΦA())) is learned by minimizing the regularized least squares objective:
(C) - - +Apall01Bs KAP
Chulle
where Xx.3 denotes the regularization parameter for the third-stage regression of KAP method. The minimizer given by
(a) -YZIAadA(a) ( dd + mndpal] -KK Aa
where $z (0z(&1)
(z(2,)) Y M
9x. KÅA E RR is the kernel matrix
over [4 and K in 's the kernel vector between the training points (á, (e and the target treatment a. As result, by substituting the estimation (1) in Equation (27 the dose-response curve can be estimated with a closed-form solution that is efficiently implemented via matrix operations 233, as summarized below Algorithm 4 For further details of this algorithm, we refer
the reader to 23).
27

---

Algorithm 4 Alternative Kernel Proxy Variable Algorithm 123
Input:
First-stage, second-stage, and third-stage: {ii, Re, ag Rie. {ūr, Rs, âe Ji=i (ūr, 2, d Parameters: Regularization parameters Ap,1, Ap,2, and X and kernel functions kz (,,), for each J (A, Z.W).
Output: Bridge function and dose response estimations: .), (xIL-))
I: For variables J' E {W: Z, A} with domain F, define the (required) first- and second- and third-stage stade kernel matrices/vectors as follows:
K F = (r Bilis E R
FFF = (N fillss E ER
K = [k> A] ERTY, FF = [ filis E CR
K [k=(B; flls E R Kip [k)(f, Riueey
Kpp - -Ix(fr. files E Ro
Kpp - [k f(] E 1R'v
2: Define the following matrices:
B (Kw Kad ppd (Kww @ Kaa).
m
B e R where B.,
(Kww Kāā + medppa)) (K - @ K An,)
B Kz B KA
T
L
T B K 22 @ K AA 17
E R
(m (
map
2 K A
M
(m)
B" K 22 @ KA
mag
N
M] E ( (m
L +
M E
7
m&p
3: The bridge
function estimation is given by ((2,a) Tim (B* K2.) o Aal
"LB"" Kzz) Kinl
7m my
4: The dose response estimation is given by
@ATR(A) 7imm (K2))diagg()) (xdp,a1) 'Ka)) Kin)
-K
7m+1 B (K2gdiag(Y) AA + top 1] K Aa)
2 K Aa
1
mi
C.3 Doubly robust kernel proxy variable algorithm
Armed with the estimation procedures for the outcome and treatment bridge functions, we now present the doubly robust estimation algorithm for the dose- respones curve. As discussed in Section (3.2), need to derive an algorithm to estimate the term B| 2, hoo a) A = a). Let {zi, wi, aj C D
denote i.i.d. observations from the distribution PP Z, W, A). We observe the following
Ei Z, a hou (W,
A=
a] = E p Z, (W,
A=a a)
E
(1)0((Z) 2 (o())
(h) @ ea(a)) Hw8HA
Hz@HA
E
(p&
h (0z (Z) & (a)) (dw(W) & da()))
Hz HzSH.
As detailed in Section (S.M. C.I.I the outcome bridge function in the KPV algorithm has the following representation (see Appendix B.3 in 11
h
3 8 (0,)))
28

---

where the set of variables i2 and âj denote the samples from first and second stage regressions of KPV algorithm, respectively. Furthermore, as detailed in Appendix (C.2). the treatment bridge function admits the representation
ma
mip
7mo+1 my (m
HaZZIaA ā.) 2 ((())))) +
(2 (5) 00 (x(āg),
which we can write without loss of generality as
W. āi) @ da(ā))
by appropriately choosing the coefficients Tu from the set { U (0). Next, notice the following
( a))(,, A=
(
ΣΣ (u ai) 8 sea(ai) R
ΣΣ (1)) (
=E
(0z @ (a())) @ (dw(W) @d((a))
HEAWHWA
ΣΣ
any mB (2ii:a āt) (á(i)) & ow(w,) 2 (at))
1,
2z @ (a) @ dw @ )dl())) HZa8HWA A A=q
ΣΣ
ans an
ai), Hz
(da(at), (Qllat)_(xao)) Ra(a)) H4
(ow(()),
A=
=Σ
any (iziwaa āt) @w(w:), (z)00)))2((22 @ dw((w)) How
13 Le
(da(m) & ddlla)) dA(a) & da(a))
A= d
Han
=ΣΣ
asym (Ziw:a āt) @ (w, (Z) @ ów(W) | A d=al) Hew
13
1.
(a)( 2 (a (a) @ da(a))
(28)
HAA
In order evaluate the of inner products above, need the conditional mean embedding (z) = E [0 @ (W) A == We approximate this term with vector-valued kernel ridge regression. In particular, under the regularity condition EZg( W) = E H. for all e Hzw, there exists operator W S2(HL,, Hzw) such that E (Z) & (w A (2)) This operator can be learned by minimizing the regularized least-squares function:
(b)(C)= - Σ (zz(zz) & dw, Co +ADR|CIIl SH
The minimizer is given by
(z,wia(a) = =Cz.)1uu(a)= = Σ (zz((z)) & @w(w;;) = -2w (KAA - tADRI) -! KAas (29)
where & =
(KAA+ + tADR!) KAa and 2zw dz((z) @ @w(w:)).
[0z(21) @ ddw(w1)
By plugging this approximation to conditional
29

---

mean embedding in Equation
28
we observe that
ΣΣ 1. 1.t
nn (Z)W: āi) @ ow(())), Σ (022(2e) @ @ow(ws)) (
New
(pa(a)) (alda). (Qa())) Har
= =ΣΣΣ no AZIWwa ār) @ (o(w()) @ @w(ws)) Hew
W 1,
(((()) (x(0,) ds(a) (aaaaa))llll) )l(a)) Han
Σs E.
Mukz ai) (ll()i))22(2)) & ella) Hzz
1,t
Σ anjbw((x)) & (w(w)) @ edla)
HWA
Σ E.p(zs, a)) a).
In conclusion, we estimate the term El (2, (2)oo(2 a] with the following expression:
Blpo(Z,a)/o A a] ~ Σ & ()))))),..))((....
(30)
where (a(a) (KAA + (ADRI) Aa The same procedure can be applied with the PMMR
representation the outcome-bridge function that will lead to the same form of estimation in Equation (30). For completeness of this section, we repeat the pseudo-code of doubly robust dose- response estimation algorithms that we give in Section (3.2) below. Note that DRKPV uses the KPV algorithm (summarized in Algorithm whereas DRPMMR employs the PMMR method (summarized in Algorithm
Algorithm Doubly Robust Kernel Proxy Variable Algorithm (Replication of Algorithm in Section 3.2)
Input: Training samples (yi, wi, 2s, 0: ki=1 C D,
Parameters: Regularization parameter ÄpR and the parameters of Algorithms and or Output: Doubly robust dose-response estimation (i0k) (a) for all a E A.
I: Compute the estimated outcome bridge function h and dose-response curve 0i(-) F(,(W, -)] with either
Algorithm 2 (3)
2: Compute the estimated treatment bridge function & and the dose-response curve dz(-) = ELY ((Z-)
= with Algorithm (4)
3: Let (=()) be given by Os(-) = E, (()))((;;; ()( :), where {:(a) [(KAA LAna KAo], for all A.
each a A. return doubly robust dose-response estimation as
(DR) ATE (a) = 0, (a) + (a) (a(.).
D Discussion on the existence of bridge functions
In this section, we present the conditions under which the outcome and treatment bridge functions
exist their respective RKHSs, following the formulations in Miao al. (5), Deaner [10], Xu et al.
13.. and Bozkurt et al. 23 To this end, we consider the conditional expectation operators defined by
E, :L2(W, Pw PW (A=a) (2,2,22144
F, (2, + L2(W, PW PW)))=))
30

---

such that
Bal(=) B(((W) 2 A a}. Fol(w) E ( - =w, ad aj,
2 2
w E W.
Here, Awm and Az]Amo denote the conditional distributions of W and Z given = a, respectively,
To ensure the existence of the bridge functions, we impose the following assumptions on the condi- tional expectation operators.
Assumption D.1 (Assumption (4) in [3]. For each E A. the operator Ea is compact operator with singular value decomposition {ME,, Dr FE )
Assumption D.2 (Assumption (10.3) in [2]). For each a E A, the operator Fa is compact operator with singular value decomposition (1F... Fa,i, Fe
To apply Picard's Theorem
63 Theorem 15.8]-as used in Lemma (2) of Xu et al. [3 and Theorem
(10.6) of Bozkurt et 23
to establish the existence of bridge functions, we make the following
additional assumptions.
Assumption D.3 (Assumption (5)
13
For each a E A, the conditional expectation lyja =
E| Z , a) satisfies
Σ
WEL
oc,
Baa
where the singular system {ME,,, &E,ns te is given in Assumption D.7.
Assumption D.4 (Assumption (10.4) in 23 For each A. the density ratio ra (W) satisfies
p(W)p(o) p(W.a)
(ras c......aaa 00,
where the singular system {MF.... FF th F. 1 is given in Assumption D.2
We now ready state the results the existence of bridge functions. The following lemma establishes the existence of the outcome bridge function.
Lemma D.5 (Lemma (2) in [3].. Suppose that Assumptions
2.2
22.3
D.
and D.3 hold. Then,
for each a E A, there exists solution to the functional equation
E[Y Z. = hhow,, A)p(w Z, dww.
The next lemma establishes the existence of the treatment bridge function. Lemma D.6 (Theorem (10.6) in (23). Suppose that Assumptions (2.2). 22.3
D.2 and D.4 hold.
Then, for each A, there exists solution the functional equation
())o(Z A p(W) p(a)
A.A.
p(W,a)
Note that the existence of the outcome bridge function hinges on Assumption (2.5 which also ensures the identifiability of the dose-response curve when using the treatment bridge function (see Theorem (266). Conversely, the existence the treatment bridge function hinges on As- sumption (2.3). which guarantees identifiability of the dose-response curve through the outcome bridge function (see Theorem (2.4)). Therefore, these completeness assumptions play distinct yet complementary roles in establishing both existence and identifiability for the respective bridge functions.
Consistency results
In this section, we provide the convergence results of the algorithms that we proposed. First, the following two sections respectively review the consistency results for outcome bridge function- and treatment bridge function-based methods from I|I 12121 23
31

---

E.1 Consistency results for treatment bridge function-based method
Here, we overview the consistency result of the treatment bridge function based on the results in For expositional clarity, we first review the results from [3]] as we leverage similar developments in presenting the consistency results for the outcome bridge function in the next section. We begin by
introducing the assumptions regarding the noise between the treatment and the outcome.
Assumption E.1 (Assumption (12.1-iv) in (23) We assume that there exists R. such that for all 2 2, PA allmsst surely,
E B|Y | A 29902
(31)
Next, we assume that the problem is well-specified, as stated in the following assumption:
Assumption E.2 (Assumption (5.1) in (23]). (1) There exists Czww,1 E S22HWA Hzz) such that
A A) = Cziw.A A (W, A); (2) There exists a solution yoo E Hz. of Equation (2) (3)
There exists Cy 2IA S2 Hz) such that E oz (Z)1A) = (
Following proof techniques from [64). Bozkurt et al. 23 show that convergence to the minimum
RKHS norm solution of the treatment bridge equation-characterizee in the following definition-is sufficient for consistency in estimating the dose-response curve using the Kernel Alternative Proxy
(KAP) method. We follow the same setting here.
Definition E.3 (Treatment bridge function solution with minimum RKHS norm, Definition (12.2) in 23). Under Assumption EE. (2)), the minimum RKHS solution to Equation ( is defined as
Po- = arg min 1|ęl|
s.t.
A),
PEHEA
with (W,A) ((A)) A(WA)
Next, we define the (uncentered) covariance operators corresponding to the three-stages of the KAP method:
Definition E.4 (Definition (12.4) in 23 The covariance operators are defined as
(Stage 1) Ex,1 E [ówa (W, A) (w, Ówa(W, (w(W))@ @ (a(())
(Stage 2) L. = E [[[(HZ]w.a A) @ (x())) @ (HZ),,aWW @ (
(Stage 3) E, E(oA(A) & (A)).
The operators Ep. Ex, and Σg.3 are self-adjoint and positive semi-definite. With Assumption (3.1 they belong to the trace class, which ensures their compactness and implies that they have countable spectrum 65
Furthermore, state the source condition (SRC) and eigenvalue decay (EVD) assumptions [66 67
Assumption E.5 (Assumption (12.6) in (23). We assume that the following conditions hold:
There exists constant B.p,1 such that for given 8p, (1,3).
Ba1-1
Isa(Hwa $ B There exists a constant B,,2 <, 0o such that for given B,, E (1,3), Bx
||Zz,,
E
GolluZA I Bp...
There exists a constant B3, <, such that for given 8p, E (1,3).
Pe.a-
[|Cy 4.3
Isa(HA. S B
Assumption E.6 (Assumption (12.9) in (3)).
We assume the following conditions hold
Let (,g,1,))22 be the eigenvalues of Exp, For some constant Cp,1 > and parameter Pe, E (0, and all 1.
Ap,1.1 S cp,12 1/pe. 1/Poa
32

---

Let (A)22,)) be the eigenvalues Ep.2. For some constant Cp.2 > and parameter Pp,2 (0, fnr all 1,
Xp,2, I Gp.21 1/p
Let (A)23,))121 be the eigenvalues p. For some constant 3. > and parameter Pp,3 E (0, fnr all 1,
1 S Cp,32" 1/P2.33
Now, we are ready to state the convergence result of the treatment bridge function. We note that Assumptions (E.2 3), (E.5/3) and (E.6.3) are not required for the bridge function consistency. However, the consistency of our proposed double robust estimators rely on these assumptions as we show in Section (S.M. E.3
Theorem E.7 (Theorem 12.21 in 23 Suppose Assumptions 3.1
E.1
E.2
1
& 2), E.S &
Pp. to 1p.1-1
2) and E.6-1 1 & 2) hold and set Xp,1
nup
Pp.1t and n&p map
where up 0.
Then,
- Up,a11 Pp.24I
bp S Bo.2-Po.2
3p,t1 then 11 &Polltzz
=0p
map
with Xx, = -
Mip
Pwi2-1 Ap.2tPoz
i. bp 2
Pela then || PollHza = Op
Seatl
mpp
with A..2 O
Map
hp controls the ratio of stage stage 2 samples. In Theorem E.7 (ii), when byp
Pe.2+Po.2 i.e.
P.2+1
when we have enough stage samples relatively stage 2 samples, the convergence rate optimal in map under the assumptions of the theorem (see Section 5 (23]
E.2 Consistency results for outcome bridge function- methods
In the following subsections, review the consistency results for the outcome bridge functions of both the KPV and PMMR methods.
E.2.1 Consistency result for kernel proxy variable
We first review the consistency result of the outcome bridge function with the Kernel Proxy Variable (KPV) algorithm that have been obtained concurrently by [2]. We assume that the problem is well-specified stated in the following assumption:
Assumption E.8 E.. (Assumption (2.1) and (8.1) in (I2)). (I) There exists CwiZ,A E S2(HZZ, Hw)
such that PwIz,a (Z, A) Cw (Z, A); (2) There exists solution by E HwA of Equation (7).
Following a similar convention to
23 we define the minimum RKHS norm solution to the outcome
bridge function equation:
Definition E.9 (Outcome bridge function solution with minimum RKHS norm). Under Assump-
tion E.8 (2)), the minimum RKHS solution to Equation Z is defined as
hg = arg min [[h]]/w
s.f. E[h(W, () 2, yp(y 2, s,a)dy,
hEHWA
Similar to 23 establish the convergence the minimum RKHS solution for the outcome
bridge function. Now, we define the covariance operators corresponding to different stages of regression of KPV.
Definition E. (Covariance Operators). The covariance operators are defined as
(Stage 1) Σn,1 E (dzı & (z4(Z, A)], Oz = Oz @d) (A() @ (wWIZZ A) 8 (A))];
(Stage 2) Σp.2 = E. [(uw]z.A
Similarly to the previous section Σ,,1 and Σn.2 are self-adjoint, positive semi-definite and trace class under Assumption (3.T). We state the source condition (SRC) and eigenvalue decay (EVD)
assumptions 67]
33

---

Assumption E.11 (SRC for KPV). We assume that the following conditions hold:
There exists constant BN,1 such that for given Bh,, (1,3).
Ph.. [C0]2,,2hh (1s.(H2A S Br.1
There exists a constant B3,2 < do such that for given Bn.2 E (1,3).
Ph. hollHA S Bb,.-
|
Assumption E.12 (EVD for KPV), We assume the following conditions hold
Let 1.),1 be the eigenvalues of Eh,1 For some constant Ch,1 and parameter Ph,1 (0, and all i 122 1,
N S 1/m.a,
Let (2n,2,)i21 be the eigenvalues of Eh. For some constant Ca,2 >> and parameter Ph,2 E (0, and for all 1211
2,,2, S Ch.2i-1/ph.
Theorem E.13 (Theorem (7) [1]). Suppose Assumptions 4
3.1
E.
E.8
E.
and
E.
hold
Ph.1 1
and set An,1 = 11 db,1 =pn.1 and nn = m2
where in >> O. Then,
P6.2- -
ih S
2h then Ii hollzwaa = Op 3
m
74,2-T
= with 2n,2 m
Ph.2+7
-h,2-1 Fh.2+Ph.2
i. if 2 Ph. then I hollzwwa Op
mp
with An.2 A m
01,2+Ph.2
Similarly to Theorem E. controls the ratio of stage to stage 2 samples. In Theorem E. E.13 (ii)..
when we have enough stage samples relatively to stage 2 samples, convergence rate is optimal in ma under the assumptions of the theorem. A similar result is obtained by Singh [12] Theorem (3))
with slightly worse requirement on stage samples in order achieve the optimal rate my..
Technically, Mastouri et al. [1] require Y to be almost surely bounded. However, we can resort the weaker assumption of sub-exponential noise, Assumption (E.I). by using the concentration
inequality of [67 Theorem (26)] as done in 64 23
E.2.2 Consistency result for proximal maximum moment restriction
Here, state the convergence result of the outcome bridge function for the PMMR algorithm, originally established by Mastouri et al. (TT]. This result plays key role in our consistency analysis of the DRPMMR algorithm. Similarly to KPV, assume that the problem is well-specified but only require Assumption (E.8 (2)).
Definition E, 14 (Outcome bridge function solution with minimum RKHS norm PMMR). Under Assumption (2),, we define
hg = arg min ||4]] ||]]Hwa bERWA
s.t.
A)] = B ((2,
Note that under Assumption (2)) there always a solution to Ej&(W, A)Φ( (,((Z,A) A)] E[Y (2, A)]. Next, we introduce the cross-covariance operator associated to PMMR.
Definition E.15 (Cross-covariance Operator). T ElozA A) @ (W,
To ensure the consistency of the PMMR outcome bridge function estimator, we require the following assumption.
Assumption E.16 (Assumption 6 in surely and E < Cy.
We assume that there exist Cy such that [Y] < Cy almost
34

---

Assumption E.17 (SRC for PMMR such that for given (1,3),
Assumption 16 in [TI). There exists a constant B < 0o
I(TTT) -5
"hol: hollowaa S B3.2.
Theorem E.18 (Theorem in IT). and set APMMR TT then
Suppose Assumptions (3.1)
16 E.8 and E.17 hold
[1 hollHaw = Op +ナ
Using a refined analysis as done in 64 it should be possible improve their result with rate depending on the eigenvalue decay of the operator T as well as relaxing Assumption (E.16 to Assumption EE.I.
Remark E.19. Notice that the convergence rate given in Mastouri et al. [11] Theorem 3] is l/h ho Haw = =pp (t min( ( However, since ? E (1,3), we have 2-1 S 2. Therefore, the
min(-) operation this consistency result is not necessary.
E.3 Consistency results for doubly robust kernel methods
In this section, we provide the consistency results of our proposed methods DRKPV and DRPMMR. Our results will rely on the consistency results of the outcome bridge function- and treatment bridge function-basee based methods that we discussed in Section (S.M. (S.M.E.2) and (S.M. E. Furthermore, we have seen that our doubly robust estimation procedures require learning the conditional mean embedding (A Cz.wIA
E [ & (w( A a). Recall that we approximate this CME via kernel ridge regression. We assume the following conditions:
Assumption E.20. There exists Cz,wIA E A2 Hzw) such that pz,wia(a) CZ.WIAỘA(o).
Assumption E.21. There exists a constant BDR < oo such that for a given 3DR (1,3).
EDr
(C2.wiA 23
(Hs, S BDR-
First, recall that we showed in Theorem 2.7 that the dose-response can be identified with
0 (a)=(a) = E [oo( aY ho a)) A a] + E(ho(W, a)). We define an intermediate quantity ( that is given by
(a) = E ((Z,a) (y i(W,a)) AA + ((((W,) a)). Noticing, by triangle inequality, that
(a) (are(a)) S (aΛTε(a) A (a)r((a) (a))(a))
will subsequently derive bounds for the two terms above in order to achieve the consistency result of our proposed methods.
Lemma E.22. Let ho and po be outcome and treatment bridge functions that satisfy Equations and 2 respectively. Then,
ATE((A) (aTe(a) = [B((p0(Z,)) (190(2,))($(Z.)))) (ho(W,a) (((W,a)--(W.a)). I
Proof.
0ATE(A)- - OATE(O) (( (W),)) = E,h((W,a))
35

---

E (+(2,a) a) (x - h(W A a) B(((Wa)]
-Blso(2,a)) - -1 a=
E -ttp:h/h..KkAA-------666220001.0 A= a|
( A A= == a +Blho(W, h(W,a)
= =11600/0000020000144444444000000 ((,,))))6((22))0000)))1) a)) a) | A = a] h(W,a)]
a) - - -E - AA + =
+ E [ - a)) AA= - E
pa a) (ho - i(W. a)) A =
= E (oo0z,,))-(Z,aa)))) - ho(W))) A a]
+ ((0(Z (lZ. - (((,a)) a)) A = a + Eloo -WWWa)
A= a]
= E E
[r ho. ( - ()
+ -
-(W,a)) | A = a] + (W, - -))
E (W A A=q a)
E (po(2,a) - (2(,,)) (bo(W,, a) - -h(w.a) a)) |A a a)
oo( (ho(Waa) -((W.) a)) |A AA=al a a) + E[, - i(W, a)
E
E ()) ()(,,a)) (Ro(W, a) - (No(W,a)-h((WW))) i((,a) A aa a)
(/ (00( ())(=//, a)dz (ho(w,,)) a) h(w h(w,,)))) plucla)dw
p(we)/((w(a)
(W, a) k(W,)) a)
E [100(Z, (2,a) a(),a)) h((,,)))(((Wa.) A a) a(( plw)dw + BBha(W, - h(w, a)]
/ (No(w,, a) h(w.a))
= E ((po(Z,)) a)) (No(W,a) (bo(W,a)h((Waa)) h().a)) A = a)
Lemma E.23. Suppose Assumption 3.1 hold and let ho: so be as in Lemma E.22 with ho E HWA and po Hz, Then,
(ATE(a)= @sTE(a)) S < @111//o -oll.
Proof. By Lemma E.22
(ATE(O) OATE(A))=
E (v)( ((2,a)) - ((W,))))((Waa)) A -
S [ a) - 2(z (W, - =q]
(ho h ohw (W) 30 ( | a s ||vo - *llloz (2) 2 (a - Alllow 3 @aaa)l]|A
[Keo (z) 0 (a(()))
x I| @l11/o hll.
36

---

Next, we need to derive a bound for (QATE(a) - (ATE(a))-
Lemma E.24. Suppose Assumptions 3.1 and E.20 hold and let ho; yo be as in Lemma E.23 Then,
[0ATE(O) - @ATe(a))
In² < - -Poll + IIsoll)) I|CyzJA - Arzj + (11- holl + l1bool) Ilāw mwll
+ n² (I.P0 (1:o0l1/Aol) |IAool + 11 - 40 ho 11 ho holllspall + 11:0- wollll - - |1C2wia Czwiall-
Proof.
We start by noting that we can write
@AT = E ( i(W,))) | A = a] + E(i(W. a)]
= (1tp::/0....000---------.//////- A A=d a]
(2, Cr ( ()) + (h. )o((a)) (p h c h callaala)))
where we define
clu) ( = E (Z) @ da(a)) @ (dw(W) 8 (a()))) | A = aq]
(32)
Similarly, we observe that we can write our doubly robust estimate
(TTE(a)= (1), Crzueλ(a) @ ()(()) (h. @ ox(l)) -
and define C(a) C to be the augmented version of the conditional mean operator estimate in Equation (29.. that is given by
(2w)aaa(a)- Σ ( & da(() @ @ @ @s((a
where E:(a) = ((KA++ODDR) _Kaa 'Ka]; Therefore, we have that
BATE(a)= - (ATE(O)= = (p. Crz da(a)) + (h pw @ pa(a))
cl
(0, (y) 2 oa) (h @ da(a)) + (p h Clo) 041
(A, (CYzIA- CYZIA) 9u 3 x)(1))
+ h (uw - kww 3 (a ( h (CheiA
(aoiaa )olla)
= p - pu: (CrzjA - CyziA) (l(() 2 (l(a))
+ po. (CrzIA- -
(4 2 @d(a))
+ ho. (uW faw) @ (oalo) + (ho, (uW faw) & @l(a)
( s h, (ClOriA Ela) c2011) @a(c))
As result,
( - @ATE(a))
37

---

~ soll + ([ ||y2 Cyzjall (lh- (1f-hin holl + |12oll) l/āw awil
(11ooll1/boll | polll|holl + hhlllllllop + 11$5 ppllll - holl [
C 22m1)) é사(a)
Note that HzAwA and HzwAA are isometric Hilbert spaces. Let be the canonical isomorphism.
V Hz @ Hw @ H. H. - Hz H. @ Hw @ HA
such that it acts on the elementary tensors as (u @ a2) = u X an 2 a2. V is a unitary isomorphism, hence,
Cha c (2wnn8a(a) Cla) (a)IRu(aa)
HZAWA HZAWA
= || ((CZwIAA))())) & (A(A) @ (al(a)) P ((zzmarka(a)) & @a(a) s )alla) (Czwladalaa) @ da(a)) @ dala)) HEWaN
Hzawa
-[(C2wIabA(U)) )d4(4) 8 da(a)
S x² [22wnsd((a) CZwlaRΛ(a)
S <n² C2lcciaa- SI(HA.H2W)
Haw
Using this inequality, we obtain the desired bound as
(T) 0ATE(A) xx² (ll@ -ool + ll |Cyz|д Crzuul + (11h- holl [lhol) 1lāw uwll
x (1s00ll/bbol |1@ holl + || hollll 11 -oolllllll ( ||zzwia Czwiall
Combining the bounds from Lemma (E.33) and E.24 we observe that
(a)T(()) (a) S x *1 *11l//0 Alll
(33)
+ x² (11+- soll + (1 ||Cyz| Crzzall + (1 holl + (lholl) llāw uwll
+ x3 (Iseolll/bol ||@ :ooll|/oll + 114 hollllsoll + 11 -olllh holl) |Czwia C2wlall
Hence, the convergence of our methods will depend on the factors 1 holl- 11 poll. -( |[CyzJa Cyziall, and ||C2wJa O2wialll We have already established the consistency results for |h holl and @ spoll, and [|ūw /w|l converges with Op M. Hence, we need to establish the converge results for |CyZJA - Õy 2)zz all and [[Z2w14 Czwjall.
Theorem E.25 (Theorem (12.22) 23 Theorem (3) 45
Suppose Assumptions
E.
E.2 3), (3.1
E.5 3) and 3E hold and take X + tp Re.3-Po.3
There is a constant J3 independent
t 2 1 and of (0, 1) such that
8ea-1
CYzIAA Crzia (.HAA)2)) ( S J/ log(5/8)
1
(6, Bp.3. Pp.s),
Vtip
is satisfied for sufficiently large tip with probability at least
5.
Theorem E.26 (Theorem (3) [55). Suppose Assumptions
E.
20
E.21 and E.6 E.63), hold and
take ADR = θ ( Tpatppa
There is constant JDR > independent of 2 1 and 5 E (0,1) 1) such
that
Spg EDRRPPS.
1|Z2II- CZwia / (sstHa)))))) < JoR log(5/8) is satisfied sufficiently large t > 1 with probability at least - 8.
( (,,,,,,,,nnn Po.3),
38

---

Recall that KPV and KAP use the number of samples (ny, ma, ) and {n., mep tip in their first, second, and third stage of regression, respectively. Following the original implementations from Mastouri et al. [1]] and Bozkurt et al. [3 we will reuse the data in their third-stage, i.e., ln l1p I,
as pointed out Algorithms (2) and (4)
Theorem E.27. Suppose the assumptions in Theorems E.
E.25 and (.2 hold. For a given training dataset DD = {Ui Wi; 2i; 0z) 1. let {nn, mh, Ln) and Amee mugs Sissshpp hup) denote the number of samples used in first-, second-, and third-stage regressions for KPV and KAP, respectively,
with 1p = top =t. Then, for DRKPV algorithm with high probability,
1f Un
Asat1
and lip
9p.2+1 Pp.2+
Set Xn.2
O
m,
Phat
and X
Ph,2+Ph.2
θ mip
Ap,27
Then,
(0TT(()) OATE(A)] = Op
AopatRop + m
80,3-1
-# 0h.2-1
Phi211
mp
te 8p.2-1
Po2t1
ii.
If 4n I
B1.+1 and lup B
B2241 De.2+Do.2"
Set Xn.2 O m,
Pa,a+1
and Xp,2
θ map "oattiin ) Then,
(0TE((a) DATE(A) = =Op (t -
Aea-1
- 26.21
80.2-1
1pppttoo + m,
9h.211 mip
Set Xn.2 1 10(1))) m,
ii.
If Up 2
B1 and up S
Bootl Bu.2t
@m/,
Fn.2 phi
and X
Ph.2
up
+ m&p Pp.2
) Then.
- -1
1/2 Ph.
tor 80.2-1
10ATE(a) DATE(a) = =0p
f + mp
Ph..2 my
iv. If
f1.t1 and up 2 Bo2tPoo2 P1.2+Pa.2
Bott1
Set Xn.2 A (mhh
Fn, 1b,
and Xp.2
+ m&p न,. ). Then,
$ Ph.thh myp
--
- 89.2-1
0Anz(a) Bare(a)) = Op Pnppttees + mp
Po.2
Proof. We combine the bounds in Lemmas E.23 and E.24 by using the convergence to the minimum norm RKHS norm bridge function solutions and discarding the faster terms, obtain
10ATE (a) OATE(() Sllço *1 hll + ||Oy ČYzial
|láw uw/l ||C2w14 ĆZwiall-
Note that the term l(åw mw/l converges with 1-2 rate hence be discarded as well. The 42,3-1
terms |CyzAA Cyziall and ||CZwJa Čzw All converges with Op t
ap.attppp with the
given regularizer parameters Àp.3 and ÅDR Theorems (E.5 and E.26 Hence, convergence will be governed by
8 (ATE(()) -ATT((a) S T0 dell Po.3-P.33
39

---

Appealing to the conditions in Theorems (E.7) and (E.3), we will have four conditions depending on the first and second stage data splitting conditions of KPV and KAP algorithms:
Suppose i. S
2 and uip
BeatPei Then, condition [i] Theorem E. and E.7 apply. Combining these bound gives
Beat1
Da 2 Ph.2
condition [i] in Theorem
(Q)T((a) ((a))-AATTEE (ATR(A)) = Op t
Bea-1 Po.3 + +m
Pb.2-2 2 Ph,2
4
30.2- Pp.2
mi
with A1.2 = 1 (m Ph.oti and A222 =
to Pp.t7
ma
ii. Suppose S 94.4 and hp
B Then, condition [i] in Theorem 13 and
condition [ii] in Theorem E.7
apply. Combining these bound gives
(0)T((a) (xTE(a)) = =0p (1
Pp.a1
-
Pb.2 1- Ph,2 mip
Ap.2-
t -mp
with An.2 1 (m Th.27T and Ap.2 = Ap.2= mpp
Jv.2 1
ii.. Suppose t 2 Pt.
91 and up S
Bupt Then, condition [ii] in Theorem EE.3 and
Buatl
condition [i] in Theorem E.7 apply. Combining these bound gives
Re.a-1 (0) (xTE(a)) = =0p 12gatoop + -m
(
Fh.2 8p.2-1 Pa.2 1h,2 my
-4
with An.2 = m
and Xp,2 6
mip
Pv.2
iv. Suppose th 2 Ph.2
9h.t1
and up BapttPe Then, condition [ii] in Theorem E.13 and
Beatl
condition [ii] in Theorem E.7 apply. Combining these bound gives
Be.-1 (6- (a)TE(a) (ATT(()) Op
2 Bu.2-1 0h 1hh2 mip
1p.atppa + +m
[m Ah.2tpbi2 and X O mm
with An.2 =
Theorem E.28. Suppose the assumptions in Theorems E.7
E.17
E.25 and (.266 hold. For
given training dataset = {ys, W4, 24, 2,, ai) let {nip: Mips too denote the number of samples used in first-, second-, and third-stage regressions of KAP algorithm with up = t.. Then, for DRPMMR algorithm with high probability
hf S Bo.atPo.a' set Ao.2 = Ap.2= -(pp) (m -
Boat1
Then,
( + 2 - mi
Pea-1
42.2-1
(0 (aTz(a))= =0)
40

---

i.. 1f Pu
Bp set 21,2
map
fo.atpoia Then,
10ATE(a) OATE(A)]= (1 Sua
2p.2-1 2 Dv.2 1x.33
SpiatRoa + 2 2-1
mp
Proof. Combining the bounds from Lemmas (F.23) and (..24), and using the convergence of the estimators to the minimum-norm RKHS bridge function solutions, we obtain the following by discarding higher-order terms:
Similar to the proof of Theorem E.27 the bound will be governed by
Re.3- Pp.3.
(ATE (a) (x) 50 =
b||l L
Now, using the conditions in Theorem EE.7 and the rate in Theorem E.17 we will have two conditions:
Suppose up I
Bp2t1 then condition [i] Theorem E. 13 applies. Hence, Bo.ttPo2"
(ATE(a)) - -AAI((a)== =0, 1 ( + 1 7-1 mip
(1
Bp.3-1
- Y Pei Pe.2tpo.2
with Av, O m&
ii Suppose up 2
1 Poott
then condition [ii] in Theorem (E.13 applies. Hence,
6,.3-1 t pp.tt + t T4T mpp
- P
00xm((a) [ATE(a)] Op
with A,, O map
Pei Paiz
Remark E.29 (Curse of Dimensionalityy. Our proposed algorithms (DRKPV and DRPMMR), along with KPV and KAP methods, rely multi-stage kernel ridge regressions. Kernel Ridge Regression known to achieve minimax-optimal rates moderate dimension, but very high-dimensiooal regimes, particularly when the input dimension d grows with the sample size t (i.e., d/t² - c some B E (0, 1)), the situation becomes more subtle 46 45
Specifically, the effectiveness of KRR suffers from two primary issues in high dimensions:
When functions lie in Sobolev classes, the achievable rate explicitly depends both smoothness and dimension. Intuifively, the regression function must be smooth enough relative to the ambient dimension d for KRR to remain consistent. More specifically, Fischer and Steimwart [67] Corollary 5 & 6] (for standard scalar KRR) and Li et al. (45) Corollary 1 & 2] (for vector-valued KRR) show that if the target function has smoothness parameter (in Sobolev sense), the optimal rate of =-2722, convergence in L2 norm is O(t==a/a)) This rate clearly demonstrates the curse of dimensionality, for fixed the bound becomes vacuous as d - cc.
Donhauser et al. [46] show that for rotationally invariant kernels (such as RBF Matérn), polynomial approximation barrier arises: the learned function effectively restricted to low-degree polynomials d grows, regardless of eigenvalue decay. This implies that consistency in high dimensions limited unless additional structural assumptions are imposed.
In our work, we do not address the curse of dimensionality: Instead, our contribution is to show that doubly robust PCL estimators can be constructed without density ratio estimation and kernel- smoothing, thereby extending practical applicability continuous and high dinenssimns treatments
41

---

where prior DR methods 121 221 fail. Our theorems remain valid in high-dimensional settings provided the smoothness and effective RKHS dimension assumptions hold. However, we acknowledge that when input dimension grows with sample size, methods based on standard RBF or Matérn kernels may indeed fail. Addressing this deeper theoretical limitation left as a future work.
Remark E.30 (Asymptotic Efficiency and Normality). While our identification formula in Theorem 2.7 is derived from the Efficient Influence Function (EIF), our resulting estimator is not classical one-step estimator [47 48) and does not automatically inherit local efficiency or asymptotic normality 128]. Unlike standard ElF-based estimators that leverage first-order orthogonality [30] between two coupled nuisance components (e.g., outcome regression and propensity score), our method involves three components- the outcome bridge, treatment bridge, and a correction term- each estimated separately. As a result, the bias terms accumulate additively, rather than canceling multiplicatively, in traditional doubly robust estimators. Establishing local efficiency would thus require further analysis, including proving asymptotic linearity, implementing cross-fitting, and verifying attainment of the semiparametric variance bound-steps we leave for future work. Similarly, asymptotic normality is not guaranteed; although the estimator is consistent and derived from an EIF, central limit theorem would require controlling higher-order remainder terms across all three nuisance
pareentttsprtrlaarrry challenging in our framework. We therefore refrain from claiming asymptotic normality the present version. Given the nonregularity typical of continuous treatment settings (59) 600 6ly, a slower-than vn convergence rate is generally expected for the parameter, point supported by both theory and our convergence results.
F
Supplementary on numerical experiments
Here, we provide additional details the numerical experiments presented in Section (5) including hyperparameter optimization procedures and supplementary experimental results.
F1
Kernel
We utilize the Gaussian (RBF) kernel for our experiments, defined as
-ll.ft 1
ks(fe, = exp
(34)
212
for fi, fy E RdF This kernel is widely used due to its boundedness, continuity, and characteristic property (68) The parameter 1 > O, known the length scale, controls the smoothness of the kernel. We set / using the commonly adopted median heuristic, which sets 12 to half the median of the pairwise squared Euclidean distances in the dataset { that is:
2 2 median 1,11-11111111 <i sisiissnn
This approach has been frequently used in causal inference applications, including Singh et al. Mastouri et al. [11), Singh [12), Xu and Gretton [62]. Bozkurt al. 23
69
In addition, we consider dimension- wise variant of the Gaussian kernel, defined as the product of one-dimensional Gaussian kernels applied to each coordinate:
d5 exp
kf
(35)
21
where f.(b) denotes the k-th coordinate of vector f. Each length scale ((k) can be set independently using the median heuristic applied to that specific dimension. We refer to the kernel in Equation (35) as the columnwise Gaussian kernel.
Following the setup in Bozkurt et al. [3], we apply the columnwise Gaussian kernel to the out- come proxy variable W within the KAP algorithm for the synthetic low-dimensional experiment in Section (5). For all other experiments and for all variables used both outcome and treatment bridge-based methods, use the standard Gaussian kernel defined in Equation 34
42

---

We also consider the Matérn kernel which provides a flexible alternative to the Gaussian kernel by introducing a smoothness parameter. The Matérn kernel between two points fi and fj is given by 70
21-V kF (F. = T(u) 11fe fzll2
222
22v K2 llfe 1551-f11a)
(36)
Here, the parameters // and are positive variables, r(-) is the gamma function, and (-)(-) is modified Bessel function of the second kind [71]. The Matérn family of kernels satisfies the
assumptions on kernels, including those in Assumption (3.1) 72. The parameter v directly controls the differentiability of the kernel function; as - oc the kernel converges to the Gaussian kernel in Equation
For practical implementation and to explore varying levels of smoothness, use cases where is half integer in the form / = p--122 which yield simplified, closed-form polynomial expressions [70)
(/ exp (- V2v
-2V (-7 fI - fsll
[(p + 1) [[2p + k=0
(p + k)! kl k)!
<VV II - fill
p-
(37)
In our ablation studies in Section (S.M. F.7 we test the performance of our method across Matérn kernels corresponding to different integer values of p.
F.2 Discussion the data splitting of KPV and KAP
We detail the data splitting strategy employed for the multi-stage KPV and KAP algorithms, which derive their samples from the full training observations (y/, Wo 2s, a )=
Algorithm Stage
Samples Used
Implementation Details
KPV (Stage 1)
{w; 2t, Je (ūr, 2,, as) mi {1
Uses n4, [t/2] random samples.
KPV (Stage 2)
Uses my = - n4.- Samples are disjoint from Stage 1.
KPV (Stage 3)
Uses the full set: l
KAP (Stage 1)
(wi, Zi, di )h (w), 0g )iq1
Uses nip [1/2] random samples.
KAP (Stage 2)
Uses mip = - nhp- Samples are disjoint from Stage
KAP (Stage 3)
(Ui, 2,, à. d
Uses the full set: t.p t.
While we use this structured splitting in our implementation, our consistency proofs in Section (E) do assume disjoint splits, similar [1 23 Indeed, using the full dataset per stage still retains theoretical consistency (e.g., Corollary I in 12 for outcome bridge-based methods).
However, structured data splitting is useful in practice for main reasons:
1. Hyperparameter Tuning:
For regression stages lacking closed-fome LOOCV solution (such as the second stage both KPV and KAP), using first-stage data as a held-out set provides clean and non-overlapping validation loss for regularizer tuning. See Section (S.M. F33 for the corresponding discussion.
2. Observation Types:
Splitting naturally accommodates the structure the different marginal distributions required across stages. For instance, the KAP first-stage utilizes the marginal distribution of (W, Z, A), whereas the third-stage requires data from the distribution of (Y, Z, A), facilitating estimation when only partial marginals are available.
F.3 Hyperparameter optimization procedures
We describe the procedures used to tune the regularization parameters λt,1, At, dp.1 Ap.2. A AMMR, and ÅDR- Specifically, A Ap,1, Aa.3. and ÄDR are selected using leave-one-out cross- validation (LOOCV), which admits closed-form solution in the case of kernel ridge regression. In
43

---

contrast, 2,, A..2 and AMMR are tuned using validation loss on a held-out set. To avoid repetition, we first review the LOOCV procedure in the general kernel ridge regression setting, followed by the closed-form expressions for the relevant regression stages of the KPV and KAP algorithms. We then present the validation loss formulas used to tune Aş.2 and AMMR-
Leave one out cross validation in kernel ridge regression: We consider the problem of estimating the conditional mean function fa( B[y X x] from observational data ( M C Rax R denote the dimensions of imputs and outputs, respectively. The kernel ridge regression
(KRR) estimator for fo is given by
arg min + Σ JEHx
(f.x=()) + Alfl|H
(38)
where H. is reproducing kernel Hilbert space (RKHS) on domain X' with the associated canonical feature map (xx()) Hx, and X >> 0 is the regularization parameter. The closed-form solution
to Equation 38 is
Y=Y + tAI) 1xxx -
(39)
where Y = 111 yel Φx (x1) ... (x(x1)), and Xxx is the kernel matrix
puted from the inputs {ax) 1=1- To select an appropriate value for X, we employ leave-one-out cross-validation (LOOCV), which assesses generalization by sequentially excluding each data point and evaluating prediction error. The LOOCV objective is defined as
LOOCV(A) = ΣΣ yo f=)0=
(40)
where f-y KRR estimator trained all data except the j-th observation. In the KRR setting, the LOOCV loss admits a closed- form expression, as given by the following result:
Theorem F.1 (Algorithm (F.1) in (69 Consider the kernel ridge regression setup introduced in Equation 38) where {= denotes the input data and {yi): denotes the corresponding oulputs. Then, the LOOCV loss is given by
LOOCV(A) = " HAYI
1 Tr (i HAYY H
(41)
where
HA = I- Xx- + MAT) E E H, diag H2) E Rext
The regularization parameter can then be selected by minimizing the LOOCV loss over predefined grid A C R:
X- = arg min
- H HAY
ACR
The proof of Theorem F.II can be found Singh al. 669 (see Algorithm (F.1)). We apply Thcorem
F.I tune the regularization parameters A, X 1, Xp.3. and ÅDR
F33.1
Hyperparameter selection for An.1. A A1p.3+ and ADR
Below, we present the application of the LOOCV tuning procedure to the regularization parameters An,1, Ap,1, 2 and ÀDR-
KPV first-stage regression: The first-stage regression the KPV algorithm, given samples
mh (ws, Ei, Te is kernel ridge regression from inputs (dz(2,) @ d to the outcomes (dhy (i. ) Therefore, the LOOCV loss for At,1 is given by
LOOCV(Ar.)) n
1 Tr [
As. HAS K Kw HJ, H ArT.))
44

---

where
H HAR = I - (K 22 3 AAA ( 3 KA + m RR xng R
H Ap,1 diag( H xx.1) "
We use logarithmically spaced grid of values in the range [5 x 10-5, -, and select the value of
1 that minimizes the LOOCV loss.
KAP first-stage regression: The first-stage regression in the KAP algorithm, using samples (W Zi, 04 me) is a kernel ridge regression from inputs (dw(())) @ ( (a)) to the outcomes
{oz (2)) thp Hence, the LOOCV loss for A,,1 is given by
HA,, LO0CV(Ap1)= Tr H
H. Apal ) K K2
7up
where
H Hx I (K Ww 2 AAA) @ KAa + I) R
H = diag H E IR=x*.,
We perform grid search over 25 logarithmically spaced values in the range [5 x 10-5, 55, 1] and select the value of Ap.1 that minimizes the LOOCV loss.
KAP thard regression: With the third-stage data {úr, 24 àz J the KAP algorithm performs kernel ridge regression from the inputs {d^(á:)) to the outcomes ( Thus, the LOOCV loss for Ap.3 is given by
LOOCV(Apa)
1 to
Tr (H HAps (K2 o rYT) HJ, ATTT))
where
Hx -- K KAkKAN A + 11dhp.1) Rlxxt,
H Hx,. diag(Hx,,.) Xty
We perform grid search over 25 logarithmically spaced values in the range [5 x 10-5, 1] and select the value of A that minimizes the LOOCV loss.
Slack term estimation: Our doubly robust estimator includes term of the form A = a). which requires estimating the conditional mean embed- ding (Zz(Z) & (w)(W) A a] as we derive in Section (S.M. C.3). Using training data (zi, ws, ws,a; 1f we fit a kernel ridge regression from the inputs { {Pa(a))) to the outputs (0,));) @ (w(w; with the regularization parameter ADR. As result, the LOOCV for ADR is given by
E (p Z, Za aW.
LOOCV(POR) = [ HAm (Kzz 2 Kww) H H H
1
where
Xne KAs(KAa + +ADDR)) Rext
H diag (Aa E RRext
For ÅDR, we use grid search over 25 logarithmically spaced values in the range (5 10-5, 1] minimize the LOOCV loss.
E33.2 Hyperparameter selection for An.2. Ap.2. and AMMR
In this section, provide details of the tuning procedures of the other regularization parameters An,2, A and AMMR- Specifically, we leverage the validation loss on a held-out set to tune these parameters.
KPV second-stage regression: The second-stage regression KPV approximates the outcome
bridge function via optimizing the objective given in Equation 20 To tune the regularization
45

---

parameter 2,, we treat the strst samples {0,, 24, a as held-out validation set and minimize the validation loss:
2 = arg min
y1
74A
where is the solution to the optimization problem in Equation 20 as described in Algorithm and Aw|z,a denotes the estimated conditional mean embedding from the first stage.
In our experiments, we used a grid of 25 logarithmically spaced values in the range (5 10-5, .].
KAP second-stage regression: The second-stage regression of KAP estimates the treatment bridge function via optimizing the objective given in Equation (26). To tune the regularization parameter 2 treat the first-stage samples { (wi, 2i, diji=1 held-out validation and minimize the validation loss:
1 Thp
((P)
(p, Pz W. (m,, a,) @d((a))
i=1
(.. 0z) @ @a(āi))
Hzd
This objective admits a closed-form expression:
AVal (x))(10)=
1 011mp
KzzC o KAA
CzKzzB
anmp am++1
(1L)) LL K2zC
KzzB Kz + K AA)
amm+1
K2 K 1
a1:mp amp+1
( LB KrzC KūA)
Here, the matrices C' and C are defined as:
C (www O Kâs + noda.11) (KwW @ Kaa),
C.j
K kww- Kūat mpdp11) (K (w Kaay);
Vj.
Thp
For full derivation, see Bozkurt et al. (23 Section (13.2.2)]. To avoid overfitting, Bozkurt et al. 23 additionally propose minimizing the validation loss augmented by model complexity penalty:
20 Tr (L" + m) LLT
Xa.2 = arg min 2Val (x) +
mip
for some fixed o > 0.
our experiments, we used a grid of 25 logarithmically spaced values the range [5 x 10 to tune A Following the complexity regularization parameters in Bozkurt et al. 23] we set og = or the synthetic low-dimensional setting as well as the legalized abortion and crime dataset, and 0 for dSprite and grade retention datasets.
PMMR regularization parameter tuning: To tune the regularization parameter AMMR we follow a procedure similar to that used in the second-stage regression of the KPV algorithm. Specifically, we set aside small validation subset from the training data {ya, w;, 2,, a, Ji and evaluate the validation loss based on the bridge function predictions. We use a grid of 25 logarithmically spaced values the range 10-5 10-3) tune AMMR, with 10% of the training set held out validation set.
F.4 Additional numerical experiments with misspecified bridge functions
higher-resolution version of the misspecification analysis presented in Section (Figure 3 provided below in Figure for enhanced legibility.
46

---

Here, we conduct an additional ablation study to evaluate the robustness of our methods under misspecification of bridge functions, building on the experimental setup from Section (5). In this experiment, we introduce higher levels of noise into the bridge function coefficients. Recall that, by the representer theorem, the bridge functions admit the following forms:
KPV: h Em شز (wy @ o사(ā,),
PMMR: = Ei= (wo)) 8 da(a;), KAP: 3 = Ei
; E mdZz,, A (wi, āt) @ (A()))
To simulate misspecification, we first train DRKPV (or DRPMMR) in the synthetic low-dimensional setting and then perturb either the outcome (KPV/PMMR) or treatment (KAP) bridge coefficients by injecting Gaussian noise.
Figures (5ah and (b)) display DRKPV results, averaged over five independent runs with standard deviation bands. In Figure (5)) outcome bridge coefficients are perturbed via aij <- Qij E where Eij (0, 0.5). Similarly, in Figure (5b), treatment bridge coefficients are jittered as Yij f- "Yy + Eij with Eij ~ (0, 0.5). Despite the misspecification. DRKPV successfully recovers the true causal function, with the slack term compensating for the introduced error.
Figures (5c) and (ch present analogous results for DRPMMR, where Gaussian noise (O, 0.5) is added to one of the learned bridge functions. As with DRKPV, DRPMMR remains robust, accurately recovering the true causal effect even under significant misspecification.
DRKPV Outcome Bridge Misspecified (0
DRKPV- Treatment Bridge Misspecifled [o]
Greand Yask
(are(o)
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-10
-0.5
0.0
0.5
2.0
(a) DRKPV, Outcome Bridge Misspecitied
(b) DRKPV, Treatment Bridge Misspecified
DRPMMR Outcome Bridge Misspecified (o² 0.2)
DRPMMR Treatment Bridge Misspecified (or 0.21
Saxad Faah
fare
(ax(a)
-05
0.5
1.0
2.0
-0.5
0.0
3.5
(c) DRPMMR, Outcome Bridge Misspecified
(d) DRPMMR, Treatment Bridge Misspecified
Figure 4: (Duplicate of Figure (3) Experimental results in bridge function misspecifications with the synthetic low-dimensional data with jittering noise sampled from (00 0.2): (a) DRKPV estimates under outcome bridge misspecification: (b) DRKPV estimates under treatment bridge misspecifica- tion; (c) DRPMMR estimates under outcome bridge misspecification; and (d) DRPMMR estimates under treatment bridge misspecification.
47

---

ORKPV-Outcome Bridge is Misspecified (o* 0.5)
DRKPV Treatmet Bridge is Misspecified (o 0.5)
Aay
2
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-1.0
-0.5
0.0
0.5
1.0
3.5
2.0
(a) DRKPV, Outcome Bridge Misspecified
(b) DRKPV, Treatment Bridge Misspecified
DRPMMR Outcome Bridge is Misspecified (o" 0.55
DRPMMR Trestment Bridge is Misspecified =0=5)
-10
-70
-25
(c) DRPMMR, Outcome Bridge Misspecified
(d) DRPMMR, Treatment Bridge Misspecified
Figure 5: Experimental results under bridge function misspecifications with the synthetic low dimensional data with jittering noise sampled from (0, 0.5): (a) DRKPV estimation when the outcome bridge is misspecified,
(b) DRKPV estimation when the treatment bridge function is misspecified, (c) DRPMMR estimation when the outcome bridge function is misspecified, (d) DRPMMR estimation when the treatment bridge is misspecified.
E5 Numerical experiments with Job Corps dataset
In this section, we present numerical experiments with the Job Corps dataset [54 55 which we accessed through the public repository provided by Singh et al. [69 //ttps://iite thub. com/ Liyuan9988/Kent (wte999888KKeeeeeala waw.eeeaaaaaaaanttoo/rrrreeee mmmttrr)....r... The U.S. Job Corps Program is an edu- cational initiative designed to support disadvantaged youth. In the corresponding observational dataset, the treatment variable captures the total number of hours participants spent in academic or vocational training, whereas the outcome variable Y measures the proportion of weeks the participant was employed during the program's second year. The covariate vector U E R2 includes demographic and sociocconomic features such gender, ethnicity, age, language proficiency, education level, marital status, houschold size, among others.
To adapt this dataset to the proximal causal learning (PCL) framework, Bozkurt et al. 23 proposed synthetic proxy generation schemes specifically crafted to test the limits of the completeness con- ditions stated in Assumptions (2.3) and (2.5). We adopt the six experimental settings introduced in their work to evaluate the performance of our proposed doubly robust estimators.
Let W = U + Gu and Z g (U71:20)) / max (T7(1201)) + where g(x) 0.8 exp(x) 1+exp(#) 0.1 is applied elementwise, and both the division and maximum operations are performed elementwise. Here, the notation U (11220)) refers to the first 20 components of the vector U. Noise terms are sampled as _O( 1) = = 1,..., 65, e ~ 1[11, Vi = 1,...,20.. 20.
Setting 1:
Setting 2: Let Z U c and W g (T(1:20) / max (U(I-20)) cw with e0 (0, Vi,
(4) ~ ~ 1] Vi.
Setting 3: Let W = U + 6 and Z 9 (T7(20:0)) / max (U7(2-10))) + 62 with ( (0, Vi, (0) ~ 1 Vi.
48

---

Setting 4: Let Z = U + and W = (U(20:10) max (U(20-101)) + fw with c ~ (0,,0 vi,
(i) ~ ~1[-1,1] Vi.
Setting 5: Let W = U + cu and Z 9 (U7(40.00) / max (U7(4060))) c2 with co ~ (0,,1) vi,
& 11-1,11 1 Vi.
Setting 6: Let Z = U + 6; and W = g (U(10:60) max (U((0:60))) + 6u with (0,, Vi, 10 ( vi.
Notice that the odd numbered settings are constructed introduce an incomplete link between the treatment proxy Z and the confounder U, combined with a nonlinear transformation. These configurations are prone to violating Assumption (2.3) that hinges the identifiability of dose-response with outcome bridge function (see Theorem (2.4)) On the other hand, even numbered settings are constructed introduce an incomplete link between the outcome proxy W and the confounder U, combined with a nonlinear transformation. These setups are likely to violate Assumption (2.5)
which underpins the identifiability guarantees of dose response with treatment bridge function (see Theorem (..6)).
Bozkurt et al. 23 showed that their treatment bridge-based method, KAP, produces estimates that are more closely aligned with those of the oracle method Kernel-ATE [69]. which assumes access to the true confounder U, compared to outcome bridge-based approaches such as KPV and KNC. They attribute this to KAP's robustness when the existence of the bridge function is violated, rather than when the assumption underlying the identifiability of the causal function is compromised. Conversely, their results suggest that KPV and KNC tend to yield estimates closer to the oracle method in scenarios where the existence of the bridge function is violated but the assumption ensuring the causal function's s identifiability remains potentially valid- highhightin their strength under different failure modes. While these findings emphasize the complementary robustness of the two classes of methods under varying conditions, determining which condition holds in practice is often difficult. In this work, we demonstrate that our proposed doubly robust estimators effectively unify the strengths of both approaches, yielding consistently strong performance across diverse scenarios.
Figure (6) presents the estimated dose-response curves produced by our proposed methods across all six experimental settings averaged over 5 different realizations with standard deviation envelopes. For comparison, we include estimates from methods based solely on outcome or treatment bridge functions, along with the oracle method Kernel-ATE from Singh et al. [69], which has access to the true confounder U. For numerical comparison, Table (T) reports the mean squared distance between each algorithm's 's estimate and that of the oracle method. Across all settings. DRKPV consistently outperforms its non-doubly robust counterpart KPV, as well as the baselines KNC and KAP. In Settings 1, 3, 4, and 5, DRPMMR also outperforms its non-doubly robust variant PMMR and the other baselines. The only exceptions are Settings 2 and 6, where PMMR performs better than DRPMMR, based on the oracle method's s dose-response estimates the ground truth.
Table 1: Mean squared distance between each algorithm's estimated dose-response curve and that of the oracle method across the six experimental settings for the Job Corps dataset described in Section (S.M.) F.5
DRKPV
DRPMMR
KPV
PMMR
KNC
KAP
Sct. 1
0.87 0.30
1.17 0.49
15.58 + 2.38
1.55 +
4.45 + 1.37
1.95 + 0.21
Set. 2
1.78
0.83
2.02
4.30 2.05
0.24
2.74 1.12
9.67
Sct. 3
2.58 0.65
0.20
1.63 0.35
3.15 6.59
Set. 4
2.29
0.43
+.31
5.25 4.201
1.94
2.11
4.12
Set. 5
+.47
0.89
10.20 + 5.32
1.26
4.45 2.29
1.83 0.71
Set. 6
0.99 0.42
+ +
4.99 2.59
0.62 0.19
2.87
2.70 + 1.48
F.6
Ablation study:
hyperparameter sensitivity on the length scale
Here, we investigate the sensitivity of our methods to the selection of the kernel length-scale hy- perparameter, I, for the Gaussian kernel used throughout our experiments. In the main paper (and detailed in Section (S.M. F.I)), we employed the median heuristic to set the length scale. This process
49

---

ATE Predictions Setting
ATE Predictions Setting 2
65.0
62.5
62
60.0
60
57.5
55.0
58
52.5
7M9
Kamel Abumacee
56
79961
Karral Abamat
50.0
54
500
750
1250
1750
a
a
(b)
ATE Predictions Setting
ATE Predictions Setting 4
65.0
62.5
60.0
57.5
10
58
35.00
Daybl Rodual
Dasbly Ratsol
52.5
Kartel Ahemacat
Karnal Abamar
FMRR
52-
FMMR BP
50.0
750
730
1250
1750
3
(d)
ATE Predictions Setting
ATE Predictions Setting 6
64
65.0
62.5
60.0
(ATE
57.5
58
55.0
Daubly Rabuat
56
Doubly Fakuai
cortel Alwoor FUMR
52.5
KPV
gV Karrel NezaiNn
50.0
Serral Megzzis
500
750
1500 1750
750
1250
1750
a
(c)
(f)
Figure 6: Dose-response estimation curves for the Job Corps experimental settings described in Section (S.M. F.5). Figures (6a)- 61 display the dose-response estimates from our proposed methods, DRKPV and DRPMMR, compared against KPV, PMMR, KNC, KAP, and the oracle method, Kernel- ATE.
involves setting equal the 0.5-quantile of the computed pairwise distances across all training data. To demonstrate the robustness of our proposed DRKPV and DRPMMR methods, we vary the quantile value used for setting for each of the involved kernels. Specifically, we test the algorithms on the synthetic oow-dimensional setting (with t = 2000) and the Legalized Abortion & Crime dataset by varying the quantile within the set (0.25, 0.4, 0.5, 0.6, 0.75, 0.9}. The corresponding results are presented in Table (2). The table demonstrates that the performance of both DRKPV and DRPMMR stable across a broad range of length-scale quantile selections, and confirms that the 0.50-quantile (median heuristic) is robust choice. Degradation in performance is only significant at the extreme high end of the scale (e.g., the 0.90-qua the synthetic dataset), indicating that our methods are generally robust to the specific selection of l.
50

---

Table 2: Ablation study hyperparameter sensitivity: Mean Squared Error (MSE) across different length-scale (4) selections for the RBF kernel. The length scale set to the specified quantile of pairwise distances.
MSE (mean + std)
Algorithm
Quantile
Synthetic Low
Abortion Crime
0.25
0.055 0.035
0.024 0.007
0,40
0.028 = 0.016
0.020 = 0.014
0.50 (Median)
0.026 = 0.015
0.018 = 0.010
DRKPV
0.60
0.027 = 0.018
0.019 = 0.012
0.75
0.023 + 0.016
0.90
0.22 = 0.096
0.021 0.018
0.25
0.026 + 0.011
0.005
0.40
0.020 + 0.011
0.015 +.0000
0.50 (Median)
0.019
0.016 =
DRPMMR
0.60
0.018 0.009
0.75
0.034 0.033
0.025 0.013
0.90
0.22 + 0.130
0.022 + 0.015
E.7 Ablation study: performance with the Matérn kernel class
We now demonstrate the performance of our proposed methods across different kernel choices by conducting experiment with the Matérn kernel class in the synthetic low-dimensional setting (with t = 2000) and the Legalized Abortion & Crime dataset.
As detailed in Section (S.M. F.I. the Matérn kernel is controlled by the smoothness parameter In this study, we focus on cases where v is a half-integer, i.e., v = + 1/2 where p is an integer. The closed form polynomial expression is given in Equation (37).
We report the performance of our proposed methods, DRKPV and DRPMMR, for varying values of the integer within the set (0. 1, 2, 3, 10). This range allows us to explore kernel smoothness from the least smooth Exponential kernel (y 1/2, p = 0) up to an approximation of the highly smooth Gaussian kernel (v 10.5, p 10).
Table 3: Ablation Study: Mean Squared Error (MSE) for various Matérn kernel selections, controlled by the smoothness parameter p + 1/2.
MSE (mean std)
Algorithm
Smoothness Parameter (p)
Synthetic Low-Dim
Abortion & Crime
( 0.5)
0.032 + 0.016
0.019 0.007
1 (v 1.5)
0.022 + 0.013
0.022 + 0.015
DRKPV
(v 2.5)
0.022 + 0.013
0.020 + 0.012
3 (v =3.5)
0.022 + 0.013
0.023 + 0.017
(v 10.5)
0.024 + 0.015
0.023 0.014
0 (v = 0.5)
0.029 + 0.011
0.016 + 0.003
1 (v 1.5)
0.019 + 0.013
0.019 + 0.009
DRPMMR
2 (v = 2.5)
0.019 + 0.013
0.018 + 0.010
(v 3.5)
0.018 + 0.013
0.018 + 0.010
10 (v = 10.5)
0.022 + 0.016
0.021 + 0.011
F88 Scalability analysis: impact of Nyström approximation on performance
We now analyze the computational complexity of our proposed DRKPV method and introduce the Nyström approximation as technique to enhance its scalability for large datasets. While we focus on DRKPV for simplicity, this discussion extends directly to DRPMMR.
51

---

The DRKPV framework comprises three primary computational components: 1) the KPV algorithm (Algorithm (2, 2) the KAP algorithm (Algorithm 4 and 3) the estimation of the slack correction term (step 3 in Algorithm (T). We break down the complexity analysis as follows:
KPV Algorithm (Stage I): Utilizes Th data samples and involves the inversion of ny nh
Gram matrix, yielding complexity dominated by the matrix inversion, O(n)))
KPV Algorithm (Stage 2): Uses my data samples and requires the inversion of m, X my Gram
matrix. Its complexity is therefore ((m)).
KAP Algorithm (Stage I): Uses nip data samples and involves the inverting an n x Rep Gram
matrix, with complexity of O
KAP Algorithm (Stage 2): The complexity is dominated by the inversion of (m, + x (mp + 1)
matrix, scaling with (m,)), where "ipp is the number of samples in the second stage.
KAP Algorithm (Stage 3): Leverages fip data samples, and its complexity is governed by
inversion of Lp hp Gram matrix, scaling with O
DRKPV Slack correction term: As outlined in Equation (6). this step requires inverting & x t Gram matrix, where / is the total number of data samples. Hence, the computational complexity is O(s)).
In conclusion, the overall computational complexity of our method, DRKPV, dominated by the inversion of the largest matrix, scaling O(t"). This is similar to the complexity of standard kernel methods like kernel ridge regression.
In general, the O(t²) computational complexity is a significant bottleneck for applications involving large datasets. To address this, we propose the Nyström approximation concrete step towards making our method more scalable.
Nyström approximation: Recall from Equations
38
39
the kernel ridge regression (KRR) uses observational data {a; ys Vi
estimator for the conditional mean function fo = E( and is given by:
- Y + LAT)
As detailed in the complexity breakdown above, the required matrix inversion of the / x matrix (Kxx + tXI) makes this solution expensive to compute, as the complexity scales cubically with the number of data points. The Nyström approximation technique (73) tackles this bottleneck by relying on a low-rank approximation to the Gram matrix Kxx. This approximation is constructed by using smaller subset of landmark points {s:),= {a:)i= where p < These landmarks typically selected using various sampling schemes 74). The closed form expression with this approximation
is given by (See Eqatuion (2) in (75)):
f(x) = Y T *
KS
(42)
where [Kxskj = (xx 83), [Kssk = (x), and [Ksa]i = kx(s;, The pri- mary computational advantage of this Nyström form is that the matrix requiring inversion, [K] KKs xs + tAKss) is now only p x p. Assuming that the number of landmarks is sub- stantially smaller than the total number of samples (p < t), the computational cost of kernel ridge regression solution is reduced from (t to O(p²t).
We apply this low-rank approximation the computationally intensive kernel regression stages in our DRKPV algorithm:
KPV (Stage 1): Kernel ridge regression from the input features {0z(Z, @ (a())) to the observations (d) ""
KPV (Stage 2): Kernel ridge regression from the input features (iw1zz a.) @ TT to
the observations (ū/ ma (-
KAP (Stage 1): Kernel ridge regression from the input features ( (W; @ a (a.)) (I
to the KAP (Stage 3): Kernel ridge regression from the input features (d:) 1i= to the observations
observations
(+))
Doubly robust correction term (Equation 29 c Kernel ridge regression from the input features
( to the observations (0z(2z) 2 ow(w)))
52

---

The second-stage regression of KAP involves a more specialized setup that deviates from standard KRR. but equally amenable to approximation techniques. Table presents synthetic low- dimensional experiments comparing the full-kemnll DRKPV to its Nyström--aprrximmate version (p 500 and p = 1000), where all regularization parameters are uniformly set 10-3, Results clearly demonstrate that the Nyström approximation effectively reduces the algorithm's runtime while maintaining or even improving prediction accuracy. Future work will also explore stronger scalable alternatives such scalable kernel ridge regression [75 and the neural network approaches [13 35 including for the second-stage regression KAP.
Table 4: Performance and scalability comparison of DRKPV with Nyström approximation in synthetic
low-dimensional experiments.
Algorithm
Data Size (t)
MSE (mean + std)
Algo Run Time (s)
DRKPY (Original)
5000
0.0097 = 0.0046
61.1 + 0.21
DRKPV (Nyström, = 500)
5000
0.0082 = 0.0033
42.5 + 0.19
DRKPV (Nyström, p 100)
5000
0.0071 = 0.0028
39.2 0.22
DRKPV (Original)
7500
0.0070 = 0.0033
194.5 + 0.35
DRKPV (Nyström, p 500)
7500
0.0085 = 0.0044
132.7 0.04
DRKPY (Nyström, p = 100)
7500
0.0091 = 0.0046
127.6 = 0.00
DRKPY (Original)
10000
0.0065 = 0.0026
453.9 + 0.51
DRKPY (Nyström, p 500)
10000
0.0059 = 0.0014
306.7 0.02
DRKPY (Nyström, 100)
10000
0.0062 = 0.0022
299.5 = 0.00
53

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction clearly state the paper's contributions and scope, namely the development of doubly robust estimators for proxy causal learning using kernel mean embeddings.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should elearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Section (6) discusses the limitations of our work and outlines future directions.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed the paper.
The authors are encouraged to create a separate "Limitations" " section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
54

---

Answer: [Yes]
Justification: We provide complete proofs and all necessary assumptions for our novel results, referencing the Supplementary Material when needed. Assumptions and results from prior work are properly cited and referenced.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We include our Python code with instructions and dependencies in the sup- plementary zip. Experimental setups, hyperparameters, and implementation details are provided in Sections (5) and S. S.M. F). with derivations and pseudocode in S.M. (C).
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibilityy can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to hosted model (e.g., in the case of a large language model), releasing of model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
55

---

The authors should provide scripts to reproduce all experimental results for the new reproducible, they 5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We include our Python code with the supplementary zip file, along with instructions to reproduce and required packages. Dataset sources and URLs are cited in Sections (5) and S.M. (F). The synthetic data can be reproduced using our code or the descriptions provided in the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://niis. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized ver- sions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide full details of our experimental setup in the main text (Section (5) and additional training details, including hyperparameter selection procedures, in the supple- mentary material (S.M.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper a level of detail that necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We report results multiple realizations and include boxplots standard deviation bars to illustrate the variability and statistical significance of our methods.
Guidelines:
The answer NA means that the paper does not include experiments.
56

---

The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1- I-sigma error bars, but one should state it. The authors should preferably report 2-sigma error bar than state that they have 96% CI. if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [No]
Justification: We did report compute resources, as our methods lightweight and all experiments can be run on standard personal computer with Python on it.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't tt make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
Our research fully complies with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consideration due laws regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
57

---

Justification: Our paper focuses on methodological and theoretical contributions and does not include discussion of broader societal impacts. We do not foresce any negative societal impacts from the research presented.
Guidelines:
The answer NA means that there is societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimesss considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications. let alone deployments. However, if there is a direct path any negative applications, the authors should point it out, For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mecha- nisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g.. pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Not applicable to the research presented in this paper.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We properly credit all datasets used: Legalized Abortion and Crime [51] 52].
Grade Retention 10, and Job Corps [54] 55). We also cite the GitHub repository from Mastouri et al. where some datasets were sourced.
Guidelines:
The answer NA means that the paper does not use existing assets.
58

---

The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provide well-structured code for our methods in the supplementary material zip file, including a README with instructions to reproduce the experiments.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their sub- missions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Our work does not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve research with human subjects or crowdsourcing.
Guidelines:
59

---

The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs were only used for grammar and writing assistance, not for core methodology or scientific content.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy htt.ge//pprrioeccc///otrrreeee/ /LLM) for what should should not be described.
60


Reviews:
[{"review_id": "brEkqJX6jI", "reviewer": "Reviewer_tteA", "summary": "When hidden factors influence a treatment's outcome, researchers can use proxy variables to estimate the true causal effect. However, standard methods require estimating complex density ratios, which can be unstable.\nThe authors developed a doubly robust method that avoids this step, leading to more stable and effective causal effect estimation with proxy variables.\n\nSpecifically, they introduces a density ratio-free, doubly robust kernel-based estimator for treatment effect estimation in the presence of unmeasured confounding, leveraging proxy variables.\n\n- Their estimator does not rely on estimating density ratios, improving robustness and practicality.\n- They employ RKHS, supporting nonparametric inference.\n- Through synthetic and real-world datasets, the proposed method outperforms standard and proximal baselines, particularly in misspecified or high-dimensional scenarios.", "questions": "See my questions in the above “Weakness” part and below: \n\n- How does the computational cost of the proposed kernel method scale with the number of samples and proxies?\n- How sensitive is the estimator’s performance to the choice of kernel function and its parameters?", "limitations": "See the limitations also in the above “Weakness” and “Questions” part.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Assumptions 2.3 and 2.5 require that the proxies are \"sufficiently informative\" relative to the unobserved confounder, which seems resasonable to me. Though in practice it might be hard to verify. \n\n**Strengths:**\n- They developed a density ratio-free approach, which is more scalable and stable for continuous or high-dimensional treatments.\n- Theorem 2.7 shows that the proposed estimators are designed to be consistent if either the outcome bridge model or the treatment bridge model is correctly specified.\n- They also demonstrate Theorem 2.7 empirically through misspecification experiments, where the estimators recover the true causal effect even when one of the bridge functions is intentionally perturbed.\n- The paper establishes uniform consistency for both proposed estimators (Theorems 4.1 and 4.2).\n\n**Weakness:**\n- The computational cost of kernel methods may make the method scale poorly with the number of samples. Thus might be impractical for large datasets.\n- RKHS methods typically rely on kernel selection and tuning. The paper does not explore how robust the method is to kernel and hyperparameter choices.\n- The DRKPV estimator, which builds on the KPV and KAP algorithms, may use data splitting for its different regression stages, but seems less discussed? Or did I missed it?", "was_revised": true, "final_justification": "My main concern was satisfactorily addressed. I remain my overall assessment.", "rebuttal": "We thank you for your detailed and constructive feedback. Below, we address your concerns point by point. If any issues remain unclear, we are glad to provide additional clarification. Should our responses resolve your concerns, we would appreciate it if you would consider adjusting your score.\n> Assumptions 2.3 and 2.5 require that the proxies are \"sufficiently informative\" ... Though in practice it might be hard to verify.\n\nCompleteness is a ubiquitous tool in statistics and has recently been employed in causal learning to establish identification. The challenge of verifying the completeness assumption is not unique to our work but is shared across proxy causal learning (PCL) and instrumental variable (IV) regression. Appendix Section 2 of [R1] discusses completeness in the context of PCL and explicitly builds on earlier foundational work such as [R2] and [R3], and many more. These references show how completeness has been formulated and justified in general statistical inference and IV problems, and [R1] uses them to argue that the completeness assumption in PCL can be similarly framed.\n\nThe general takeaway, as emphasized in [R1], is that completeness holds for a wide range of semiparametric and nonparametric models, although it is typically not testable—even when all relevant variables are observed [R2]. In addition, [R1] provides concrete examples where completeness is satisfied. For instance, [R3] shows that when proxies and confounders are continuously distributed and the proxy dimension exceeds the confounder dimension, completeness is ensured under mild regularity conditions. Thus, in PCL it is standard to assume completeness for a broad class of models. Intuitively, completeness formalizes the requirement that proxies be sufficiently informative about the confounders, highlighting the importance of collecting rich sets of proxy variables in observational studies to mitigate unobserved confounding.\n\nWe will incorporate an extended discussion in our updated manuscript to explicitly address this point and strengthen clarity.\n\n[R1] Miao, W., Hu, W., Ogburn, E. L., and Zhou, X.-H. (2022), “Identifying Effects of Multiple Treatments in the Presence of Unmeasured Confounding,” Journal of the American Statistical Association.\n\n[R2] Canay, I. A., Santos, A., and Shaikh, A. M. (2013) “On the testability of identification in some nonparametric models with endogeneity.” Econometrica.\n\n[R3] Andrews, D. W. (2017) “Examples of L2-complete and boundedly-complete distributions.” Journal of Econometrics.\n> **Weaknesses and Questions**: computational cost of the proposed kernel method and its scalability to larger sets\"\n\nWe will add a detailed discussion in Appendix F about the computational complexity and scalability of our methods. Briefly, our method builds on multi- and single-stage kernel regressions from KPV/PMMR and KAP, each of which involves kernel matrix inversions. This results in cubic complexity with respect to the sample size. Specifically, KAP stages have complexities: (i) first-stage $O(n_\\varphi^3)$, (ii) second-stage $O((m_\\varphi + 1)^3)$, (iii) third-stage $O(t_\\varphi^3)$. Similarly, KPV stages have complexities: (i) $O(n_h^3)$, (ii) second-stage $O(m_h^3)$. Here, $n_\\varphi, m_\\varphi, t_\\varphi, n_h, m_h$ denote the number of samples used in different stages of KAP and KPV algorithms. Our estimator also requires computing the slack prediction (Appendix C.3), which adds a $O(t^3)$ cost, where $t$ is the sample size of the whole dataset. Overall, the complexity of our proposed method is $O(t^3)$, dominated by the largest kernel inversion.\n\nTo address scalability, we provide an additional set of numerical experiments where we incorporate the Nyström approximation [R4] in all kernel regressions of DRKPV (KPV first/second stage, KAP first/third stage, slack term estimation). The second-stage regression of KAP involves a more specialized setup that differs from standard kernel ridge regression, but it is equally amenable to approximation techniques. Here, we demonstrate substantial scalability gains by applying Nyström to the kernel regression stages; extending it to the KAP second stage is feasible and will be included in future versions. We present synthetic low-dimensional experiments (Section 5) comparing DRKPV with Nyström (500/100 landmarks) to the full kernel version, where each regularization parameter is set to $10^{-3}$. Results (to be added in Appendix F) show Nyström makes DRKPV far more scalable. Future work will also explore stronger approximations such as scalable kernel ridge regression [R5] and neural mean embedding [R6], including for the second-stage regression of KAP. In this paper, we focus on providing a robust and uniformly consistent kernel method, leaving scalable variants for future work.\n\n|Algorithm| Data Size|MSE (mean±std)|Algo Run Time (mean±std)|\n|-|-|-|-| \n|DRKPV (original)| 5000|0.0097±0.0046|61.1±0.21|\n|DRKPV (Nyström, 500 landmarks)|5000|0.0082±0.0033|42.5±0.19|\n|DRKPV (Nyström, 100 landmarks)|5000|0.0071±0.0028|39.2±0.22| \n|DRKPV (original)|7500|0.0070±0.0033|194.5±0.35|\n|DRKPV (Nyström, 500 landmarks)|7500|0.0085±0.0044|132.7±0.04|\n|DRKPV (Nyström, 100 landmarks)|7500|0.0091±0.0046|127.6±0.00|\n|DRKPV (original)|10000 |0.0065±0.0026|453.9±0.51|\n|DRKPV (Nyström, 500 landmarks)|10000|0.0059±0.0014|306.7±0.02|\n|DRKPV (Nyström, 100 landmarks)|10000|0.0062±0.0022|299.5±0.00|\n\n[R4] Williams, C. K. I. and M. Seeger (2001). “Using the Nyström Method to Speed Up Kernel Machines”. NeurIPS 13.\n\n[R5] Meanti, G., L. Carratino, L. Rosasco, and A. Rudi (2020). “Kernel methods through the roof: handling billions of points efficiently”. NeurIPS 34.\n\n[R6] Liyuan Xu and Heishiro Kanagawa and Arthur Gretton, (2021). \"Deep Proxy Causal Learning and its Application to Confounded Bandit Policy Evaluation\". NeurIPS 35.\n> **Weaknesses and Questions**: \"sensitivity of the estimator’s performance to the choice of kernel function and its parameters\"\n\nTo evaluate robustness, we conducted further experiments across (i) different kernel length scales (quantiles of pairwise distances) and (ii) different Matern kernels (with varying smoothness parameter $p$).  For Gaussian kernels, our main results used the common median heuristic (0.5 quantile). Tables below report DRKPV and DRPMMR results on synthetic low-dimensional and legalized abortion/crime datasets. We will add these ablation studies in Section F.\n*   Across quantiles, performance remains stable for a broad range, with degradation only at extreme settings (e.g., 0.9 quantile). \n*   Across Matern kernels, our methods maintain robust performance for varying $p$.\n\nTable illustrating performance over different length-scale selections of the RBF kernel in the synthetic low-dimensional dataset.\n| Algorithm |$l$-quantile | MSE ± Std|\n|-|-|-|\n|DRKPV|0.25|0.055±0.035|\n|DRKPV|0.40|0.028±0.016|\n|DRKPV|0.50|0.026±0.015|\n|DRKPV|0.60|0.027±0.018|\n|DRKPV|0.75|0.034±0.019|\n|DRKPV|0.90|0.22±0.096|\n|DRPMMR|0.25|0.026±0.011|\n|DRPMMR|0.40|0.020±0.011|\n|DRPMMR|0.50|0.019±0.012|\n|DRPMMR|0.60|0.024±0.018|\n|DRPMMR|0.75|0.034±0.033|\n|DRPMMR|0.90|0.22±0.130|\n\nTable illustrating performance over different length-scale selections of the RBF kernel in Legalized Abortion and Crime dataset.\n|Algorithm|$l$-quantile|MSE ± Std|\n|-|-|-|\n|DRKPV|0.25|0.024±0.007|\n|DRKPV|0.40|0.020±0.014|\n|DRKPV|0.50|0.018±0.010|\n|DRKPV|0.60|0.019±0.012|\n|DRKPV|0.75|0.023±0.016|\n|DRKPV|0.90|0.021±0.018|\n|DRPMMR|0.25|0.016±0.005|\n|DRPMMR|0.40|0.015±0.008|\n|DRPMMR|0.50|0.016±0.009|\n|DRPMMR|0.60|0.018±0.009|\n|DRPMMR|0.75|0.025±0.013|\n|DRPMMR|0.90|0.022±0.015|\n\nTable illustrating performance over different Matern kernel selections in synthetic low dimensional dataset.\n| Algorithm | $p$| MSE ± Std|\n|-|-|-|\n|DRKPV|0|0.032±0.016|\n|DRKPV|1|0.022±0.013|\n|DRKPV|2|0.022±0.013|\n|DRKPV|3|0.022±0.013|\n|DRKPV|10|0.024±0.015|\n|DRPMMR|0|0.029±0.011|\n|DRPMMR|1|0.019±0.013|\n|DRPMMR|2|0.019±0.013|\n|DRPMMR|3|0.018±0.013|\n|DRPMMR|10|0.022±0.016|\n\nTable illustrating performance over different Matern kernel selections in Legalized Abortion and Crime dataset.\n|Algorithm|$p$|MSE ± Std|\n|-|-|-|\n|DRKPV|0|0.019±0.007|\n|DRKPV|1|0.022±0.015|\n|DRKPV|2|0.020±0.012|\n|DRKPV|3|0.023±0.017|\n|DRKPV|10|0.023±0.014|\n|DRPMMR|0|0.016±0.003|\n|DRPMMR|1|0.019±0.009|\n|DRPMMR|2|0.018±0.010|\n|DRPMMR|3|0.018±0.010|\n|DRPMMR|10|0.021±0.011|\n> **Weaknesses**: The DRKPV estimator, ..., may use data splitting for its different regression stages, but seems less discussed? Or did I miss it?\n\nWe agree that adding more discussion would help readers better understand the implementation, and we will devote a subsection in Appendix F to this. Specifically, our estimator uses the same data splitting setup as KPV/KAP. For a dataset of $t$ samples $\\{a_i, y_i, z_i, w_i\\}_{i= 1}^t$, the methods split data across stages:\n*    KPV: $n_h$ and $m_h$ for its first- and second-stage, respectively,\n*    KAP: $n_\\varphi$, $m_\\varphi$, and $t_\\varphi$ for its first-, second-, and third-stage, respectively.\n\nUsing the full dataset per stage still retains consistency (see Corollary 1 in [11] for outcome bridge-based method). Our consistency proof (Appendix E) does not assume disjoint splits, nor do [10] and [19]. However, data splitting becomes useful in practice due to two reasons:\n*    For KAP, tuning the second-stage regularizer (no closed-form LOOCV) requires held-out data. In implementation, we utilize first-stage data as a held-out set to evaluate the validation loss and tune the regularizer in second-stage regression.\n*    Splitting accommodates different observation types: e.g., KPV's first-stage uses $\\{A, W, Z\\}$, while the third-stage uses $\\{A, Y, Z\\}$, enabling estimation when only partial marginals are available rather than the full joint distribution.\n> **Paper Formatting Concerns**: Figures 2, 3 may be difficult to read when printed or viewed on small screens.\n\nWe will enlarge the font sizes in the figures to ensure readability."}, {"review_id": "cu82W4oGJ0", "reviewer": "Reviewer_j1Ba", "summary": "This paper addresses causal function estimation in the Proxy Causal Learning (PCL) framework, where confounders are unobserved but proxy variables are available. The authors propose two novel kernel-based doubly robust estimators that integrate outcome and treatment bridge methods. These estimators avoid kernel smoothing over the treatment and density ratio estimation, making them especially effective for continuous or high-dimensional treatments. Using kernel mean embeddings, the method achieves closed-form solutions and strong consistency. Empirical results on PCL benchmarks show superior performance over existing methods, including previous doubly robust approaches.", "questions": "1.\tIn Supp lines 939-945, the author mention that the method of utilizing IF is locally effective. In Theorem 2.7, the author derives an identification formula that inherit the dr property. Will the proposed estimators also inherit local efficiency, or at least satisfy asymptotic normality? Given the nonregularity of continuous treatment estimators, should we expect slower-than-$\\sqrt{n}$ asymptotic normal.\n2.\tGiven that kernel-based methods are involved in parts of the estimation, can the authors clarify the computational efficiency of the proposed approach?", "limitations": "The work does not have potential negative societal issues.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1.\tThe proposed estimators avoid kernel smoothing over the treatment or density ratio estimation, making them particularly suitable for continuous or high-dimensional treatment settings.\n2.\tBoth estimators are shown to possess the doubly robust (DR) property, which enhances their robustness to model misspecification.\n3.\tThe paper provides uniform consistency guarantees and demonstrates empirical performance, offering a compelling combination of theoretical justification and practical utility.\n\nWeaknesses:\n1.\tThe idea lacks novelty. The core idea of combining estimators from [10] and [19] to construct a new DR estimator is similar to the classical causal inference literature. Notably, the claimed benefits—such as avoiding kernel smoothing and handling high-dimensional treatments—are inherited from [19].\n2.\tBoth [10] and [19] derive convergence rates for their respective estimators. Does the proposed DR estimator—formed via their combination—achieve a similar rate? If so, is it the minimum of the two, or can it attain a faster rate under certain conditions?\n3.\tOne key advantage of DR estimators is inference. However, the paper does not analyze the asymptotic distribution of the proposed estimator, and thus it is unclear whether valid confidence intervals can be constructed.\n4.\tIn Lines 276–278, the authors mention that different stages require different samples, obtained via either data splitting or reuse. If data splitting is applied, sample inefficiency might be a concern. If reuse is employed, what are the implications for theoretical guarantees (e.g., bias or overfitting)? This trade-off requires further discussion.", "was_revised": true, "final_justification": "Most of my concerns have been addressed by the author's rebuttal. This is a nice work, and I would like to see it in the conference.", "rebuttal": "We greatly appreciate your valuable and detailed feedback. Below, we respond to your questions and concerns point by point. Should any points remain unclear, we are glad to provide additional explanation. If our clarifications have resolved your concerns, we would be grateful if you would consider increasing your score.\n> **Weaknesses**: ... the claimed benefits—such as avoiding kernel smoothing and handling high-dimensional treatments—are inherited from [19].\n\nWhile we acknowledge that advantages such as avoiding kernel smoothing and accommodating high-dimensional treatments are shared across kernel-based PCL frameworks [19, 10, 11], these stem from conditional mean embeddings (CMEs) in reproducing kernel Hilbert spaces (RKHSs).\n\nOur contribution goes substantially beyond this: we introduce the **first doubly robust, density ratio-free PCL framework** using kernel machines. Specifically:\n* **Novelty of our approach**: Unlike [19], which estimates the causal function solely via integration of the treatment bridge function, our framework leverages the double robustness property. This formulation avoids density ratio estimation and ensures **robustness to misspecification of either bridge function**, a point substantiated by our empirical results.\n* **Handling high-dimensional treatments**: Because our identification relies on conditional means rather than joint expectations with indicator functions or kernel smoothing, our methods naturally extend to higher-dimensional treatments. Other kernel-based frameworks [10, 11, 19] share this advantage, but unlike ours, they are not doubly robust. \n* **Empirical validation**: Our experiments show that DRKPV and DRPMMR consistently outperform existing kernel-based methods including KAP [19], KPV/PMMR [10] across diverse PCL benchmarks including high dimensional treatment settings. Furthermore, doubly robust method PKDR [18], which incorporates kernel smoothing, fails to operate in high-dimensional treatment settings since its implementation is limited to univariate treatments. This limitation reflects the instability of ratio estimation and poor approximation of indicator functions in higher dimensions.  In contrast, our methods successfully estimate causal functions in challenging settings such as the dSprite dataset, which features high-dimensional ($64 \\times 64$ image) treatments.\n> **Weaknesses**: Both [10] and [19] derive convergence rates for their respective estimators. Does the proposed DR estimator—formed via their combination—achieve a similar rate? If so, is it the minimum of the two, or can it attain a faster rate under certain conditions?\n\nOur method achieves strong uniform consistency and, under some smoothness conditions, can converge faster than both KAP and KPV individually. In Appendix F we derive rates for both DRKPV and DRPMMR. For instance, under condition (iv) of Theorem E.27, DRKPV converges with $$O_p(t^{-\\frac{1}{2}\\frac{\\beta_{\\varphi, 3}-1}{\\beta_{\\varphi, 3} + p_{\\varphi, 3}}} + m_h^{-\\frac{1}{2} \\frac{\\beta_{h,2} - 1}{\\beta_{h,2} + p_{h,2}}} m_{\\varphi}^{-\\frac{1}{2}\\frac{\\beta_{\\varphi, 2}-1}{\\beta_{\\varphi, 2} + p_{\\varphi, 2}}}).$$\nHere, $\\beta_{\\cdot}$ denote smoothness parameters and $p_{\\cdot}$ the effective RKHS dimensions in the corresponding regression stages. By comparison, KAP achieves $O_p(t^{-\\frac{1}{2}\\frac{\\beta_{\\varphi, 3}-1}{\\beta_{\\varphi, 3} + p_{\\varphi, 3}}} + m_{\\varphi}^{-\\frac{1}{2}\\frac{\\beta_{\\varphi, 2}-1}{\\beta_{\\varphi, 2} + p_{\\varphi, 2}}})$ (see Theorem 5.4 in [19]), while KPV achieves $O_p(m_h^{-\\frac{1}{2} \\frac{\\beta_{h,2} - 1}{\\beta_{h,2} + p_{h,2}}})$ (see Proposition (1) in [10]). This indicates that our method converges faster than the KAP algorithm. Furthermore, our method converges faster than KPV if $t^{-\\frac{1}{2}\\frac{\\beta_{\\varphi, 3}-1}{\\beta_{\\varphi, 3} + p_{\\varphi, 3}}} < m_h^{-\\frac{1}{2} \\frac{\\beta_{h,2} - 1}{\\beta_{h,2} + p_{h,2}}}$. The best achievable rate in this setting is $O(t^{-1/4})$, a typical nonparametric rate. This arises because our identification is based on a conditional expectation involving the treatment bridge and outcome bridge functions, rather than the expectation of a non-centered influence function. While the latter may theoretically yield faster rates, in practice such approaches suffer from density ratio estimation and kernel smoothing.  The advantages of our approach are twofold:\n* **Uniform convergence**: Unlike most doubly robust frameworks, which prove only pointwise consistency, our results (DRKPV and DRPMMR) establish uniform convergence. Depending on regression smoothness, this can yield faster rates than both outcome- and treatment-bridge methods. \n*   **Empirical performance**: As shown in Section 5, DRKPV and DRPMMR consistently outperform prior methods across benchmarks.\n> **Weaknesses**: If data splitting is applied, sample inefficiency might be a concern. If reuse is employed, what are the implications for theoretical guarantees (e.g., bias or overfitting)? This trade-off requires further discussion.\n\nIn the camera-ready version, we will include a detailed discussion of sample splitting in Appendix F. Our estimator uses the same splitting setup as KPV and KAP:\n*    KPV: $n_h$ and $m_h$ for its first- and second-stage, respectively,\n*    KAP: $n_\\varphi$, $m_\\varphi$, and $t_\\varphi$ for its first-, second-, and third-stage, respectively.\n\nUsing the full dataset in each stage also preserves consistency (see Corollary 1 in [11] for an outcome bridge–based method). Our consistency proofs in Appendix E never require disjoint sample splits; under Theorem E.27 andE.28 with optimal regularization, our methods converge to the true causal structural function. Splitting is nevertheless useful in practice for two reasons:\n*    Regularization tuning: In KAP, the second-stage regularizer (no closed-form LOOCV) requires held-out data. We use first-stage data for validation and tuning.\n*    Different observation types:For example, KPV’s first stage uses $\\{A, W, Z\\}$, while the second-stage uses $\\{A, Y, Z\\}$. Splitting thus enables estimation when only partial marginals are available.\n> **Weaknesses and Questions**: ... the paper does not analyze the asymptotic distribution of the proposed estimator. Will the proposed estimators also inherit local efficiency, or at least satisfy asymptotic normality? Given the nonregularity of continuous treatment estimators, should we expect slower-than-$\\sqrt{n}$ asymptotic normal.\n\nWhile Theorem 2.7 provides an identification formula derived from the efficient influence function (EIF), and our estimator resembles a one-step correction, we emphasize that it is not a classical one-step estimator and therefore does not automatically inherit local efficiency. Unlike standard EIF-based estimators that leverage first-order orthogonality between two coupled nuisance components (e.g., outcome regression and propensity score), our method involves **three components**—the outcome bridge, treatment bridge, and a correction term—each estimated separately. As a result, the bias terms accumulate additively, rather than canceling multiplicatively, as in traditional doubly robust estimators. Establishing local efficiency would thus require further analysis, including proving asymptotic linearity, implementing cross-fitting, and verifying attainment of the semiparametric variance bound—steps we leave for future work.\n\nSimilarly, asymptotic normality is not guaranteed. Although the estimator is consistent and based on an EIF-derived formula, a central limit theorem would require controlling higher-order remainder terms across all three nuisance components—particularly challenging in our kernel-based framework. We therefore refrain from claiming asymptotic normality in the present version, though we believe it could be established in future works.\n\nDespite this, our estimators offer several practical advantages: they yield closed-form solutions, support uniform consistency guarantees, exhibit double robustness to outcome or treatment bridge misspecification, and handle continuous or high-dimensional treatments more naturally than previous methods. Empirically, they consistently outperform existing baselines.\n\nFinally, we agree with the reviewer that in the continuous treatment setting, one should generally expect slower-than-$\\sqrt{n}$ convergence due to the nonregularity of the parameter—a point supported by both theory and our convergence results.\n\nWe thank the reviewer for these questions and we will subsequently add a discussion paragraph in the paper.\n> **Questions**: ... computational efficiency of the proposed approach?\n\nThe main computational bottleneck is kernel matrix inversion, so—as with other kernel ridge regression–based methods—our estimators scale cubically in the sample size ($O(t^3)$). As detailed in our response to Reviewer tteA (see Complexity and Scalability), both DRKPV and DRPMMR fall into this regime; we will expand this discussion and include breakdowns in Appendix F. To mitigate scalability, we applied the Nyström approximation [R1] across all DRKPV kernel regressions (KPV first/second, KAP first/third, and slack prediction), achieving substantial runtime reductions while maintaining accuracy. Detailed results are provided in our response to Reviewer tteA, and we will add these ablations in Appendix F. We plan to explore stronger approximations such as scalable kernel ridge regression [R2] and neural mean embedding methods [R3].\n\n[R1] Williams, C. K. I. and M. Seeger (2001). “Using the Nyström Method to Speed Up Kernel Machines”. NeurIPS 13.\n\n[R2] Meanti, G., L. Carratino, L. Rosasco, and A. Rudi (2020). “Kernel methods through the roof: handling billions of points efficiently”. In: NeurIPS 34.\n\n[R3] Liyuan Xu and Heishiro Kanagawa and Arthur Gretton, (2021). \"Deep Proxy Causal Learning and its Application to Confounded Bandit Policy Evaluation\". In NeurIPS 35."}, {"review_id": "CBLJKtrm9p", "reviewer": "Reviewer_R2kv", "summary": "The work  presents a way of estimating causal functions under the proxy causal learning (PCL) framework. They account for confounders by proxy of variables for treatment (Z) and outcome (W). The two estimators presented are both of the doubly  robust (DR) kernel-based type, and combine the benefits of  both Z and W bridge-based estimattors. \nThe density-ratio freedom implies they they circumvent the need of the indicator finction or kernel smoothing over treatment variables.The experiments show improved performance over sysnthetic high-dimensional data.", "questions": "Between Z/W, one may sometimes be information-poor, and the other the more informative proxy. How does that play into the calculus?", "limitations": "I'm not able to determine those critically!", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Kernel mean embeddings help with the closed-form and strongly-consistent guarantees. \n\nIt appears that the work is an important step in measuring confounders.\n\nExperiments are on high dim (quasi continuous)  data, but not on copious volumes.  Kernel computation may be a bottleneck in that.", "was_revised": true, "final_justification": "I cnceded to the arguments in the comments to my review.", "rebuttal": "We sincerely appreciate your thoughtful and valuable feedback. Below, we respond to your comments point by point. If anything remains unclear, we would be happy to elaborate further. We hope our responses address your concerns and reinforce your positive evaluation of our work.\n> **Questions:** Between Z/W, one may sometimes be information-poor, and the other the more informative proxy. How does that play into the calculus?\n\nIn Appendix F.4, we conduct simulations designed precisely for such scenarios.\n\n**Setup**: The JobCorps dataset (an observational study) is adapted to the PCL framework by generating synthetic proxies with varying informativeness. The treatment variable $A$ measures hours spent in training, while the outcome $Y$  is the proportion of weeks employed in the second year. Confounders $U\\in \\mathbb{R}^{65}$ include demographic and socioeconomic features. [19] proposed proxy-generation schemes to test the completeness assumptions (2.3 and 2.5), which we adopt to evaluate our doubly robust estimators. Based on [19],we expect the behavior:\n*    Outcome-bridge methods perform better when $W$ is less informative, since the identification via treatment-bridge is challenged (Assumption 2.5).\n*    Treatment-bridge methods perform better when $Z$ is less informative, since identification via outcome-bridge is challenged (Assumption 2.3).\n\nAcross six proxy-informativeness settings, we compare DRKPV and DRPMMR with outcome-bridge and treatment-bridge baselines. Results (Table 1, Figure 5) in Appendix F.4 show that even when one class of methods underperforms, our doubly robust estimators consistently provide strong performance by effectively combining the strengths of both approaches.\n\nThis robustness arises from the double robustness property: our methods remain valid as long as either the outcome bridge or the treatment bridge is well specified.\n\n> **Strengths And Weaknesses:** Experiments are on high dim (quasi continuous) data, but not on copious volumes. Kernel computation may be a bottleneck in that.\n\n**Complexity Analysis**: The main computational bottleneck is kernel matrix inversion. As detailed in our response to Reviewer tteA (see Complexity and Scalability discussion), our method scales cubically with sample size due to multi- and single-stage kernel regressions in KPV, KAP, and PMMR. Specifically, both of our methods have complexity $O(t^3)$, where $t$ is the number of observations. In Appendix F, we will expand this discussion and include complexity breakdowns.\n\n**Scalability**: To mitigate this, we conducted additional experiments using the Nyström approximation [R1]. We discuss this extensively in our response to Reviewer tteA, which provides detailed Nyström-based results. Briefly, we applied it across all DRKPV kernel regressions (KPV first/second, KAP first/third, slack prediction), showing substantial runtime reductions while maintaining accuracy. DRPMMR can be approximated similarly. We will include these ablation studies in Appendix F. For full results and explanations, please see our response to Reviewer tteA. In future work, we aim to apply more advanced techniques, including scalable kernel ridge regression [R2] and neural mean embedding framework [R3].\n\nReferences\n\n[R1] Williams, C. K. I. and M. Seeger (2001). “Using the Nystrom Method to Speed Up Kernel Machines”. In: NeurIPS 13.\n\n[R2] Meanti, G., L. Carratino, L. Rosasco, and A. Rudi (2020). “Kernel methods through the roof: handling billions of points efficiently”. In: NeurIPS 34.\n\n[R3] Liyuan Xu and Heishiro Kanagawa and Arthur Gretton, (2021). \"Deep Proxy Causal Learning and its Application to Confounded Bandit Policy Evaluation\". In NeurIPS 35."}, {"review_id": "KTc9mvX47Y", "reviewer": "Reviewer_oWY2", "summary": "This paper proposes two kernel-based doubly robust (DR) estimators that leverage both outcome and treatment bridge functions, enabling them to handle continuous and high-dimensional covariates. Following a framework similar to Cui et al. [17], the estimators integrate the KPV/PMMR method [15] for estimating $h$ and KAP [19] for estimating $\\varphi$, thereby constructing the DR estimator. The authors establish consistency for the proposed estimator and support their approach with extensive simulation studies.", "questions": "1. **Figure 1**\n\nThe arrow from U to W should be solid, and the arrow from Y to W should be dashed. \nConsider an example where $U \\sim N(0,1)$, $W \\sim N(0,1)$, $A \\sim N(U,1)$, and $Y \\sim N(A+W,1)$. \nThis example is compatible with Figure 1, but we find $E[U|W,A=0] = E[U|A=0]=0$, indicating that $\\ell(U) = U$ serves a counterexample of Assumption 2.5.\nThus, if $U$ and $W$ are not associated, Assumption 2.5 cannot be satisfied, which requires a solid arrow between $U$ and $W$. \nFor further justification, see Table A.1 in Tchetgen Tchetgen et al. (2024, Statistical Science), which also specifies the requirement for a direct connection between $U$ and $W$.\n\n\n2. **Theorem**\n\nWhat notion of convergence is used in the theorems?\n\nThe paper includes high-dimensional simulations, but it’s unclear whether Theorems 4.1 and 4.2 are valid in such settings. How is the curse of dimensionality addressed? Please clarify the assumptions under which these theoretical results extend to high-dimensional cases.\n\n\n3. **Reference**\n\nIn explaining Equation (2), it would be helpful to cite foundational work on stabilized weights for continuous treatments (e.g., Robins et al., 2000, Epidemiology). This would situate your methods more clearly in the literature and justify the weighting strategy.\n\n4. **Contribution**\n\nThe fundamental concept of integrating KPV/PMMR with KAP is sound. Nevertheless, the manuscript fails to clearly delineate whether this integration yields significant theoretical or methodological innovations. Notably, [19] also asserts applicability to continuous treatment and high-dimensional settings, suggesting that these aspects are not exclusive contributions of the current paper. Consequently, the paper's contribution, as articulated, appears to be primarily confined to the development of the DR property. Consequently, the paper's contribution, as articulated, appears to be primarily confined to the development of the DR property. The paper should place clearer emphasis on why this DR development is a nontrivial task.\n\n5. **Typo**\n\nLine 320, Page 8: $f_{ATE} \\rightarrow \\theta_{ATE}$", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Quality**\n\nThe submission demonstrates technical soundness, presenting well-substantiated claims grounded in rigorous theoretical analysis and simulation-based evidence.\n\n**Clarity**\n\nThe paper is clearly written and well-organized.  \n\n**Significance**\n\nThe work addresses a topic of considerable interest, as many practical causal inference problems involve continuous treatments. The high-dimensional extension is also promising, though the theoretical justification in this context requires further clarification.\n\n**Originality**\n\nWhile the paper builds on ideas from [15], [17], and [19], the novelty of its contributions is not clearly delineated. The main methodological development appears between Line 250 (Page 6) and Line 270 (Page 7), but this section seems to offer only incremental innovation. Furthermore, some relevant literature is omitted, reducing the clarity of the paper's positioning relative to existing work.", "was_revised": true, "final_justification": "I appreciate the authors’ engaging discussion. The work makes a valuable contribution to the literature, and I have accordingly raised my evaluation score.", "rebuttal": "We are grateful for your careful and detailed review. In what follows, we address your concerns point by point. If any issues remain open, we will gladly provide further clarification. If our responses satisfactorily resolve your concerns, we would be thankful if you would consider adjusting your score.\n> **Originality**: The main methodological development appears between Line 250 (Page 6) and Line 270 (Page 7), but this section seems to offer only incremental innovation. \n> **Questions**: 4-) The fundamental concept of integrating KPV/PMMR with KAP is sound. Nevertheless, the manuscript fails to clearly delineate whether this integration yields significant theoretical or methodological innovations. ... The paper should place clearer emphasis on why this DR development is a nontrivial task.\n\nAlthough our work builds on prior literature, our contributions extend substantially beyond the section you referenced, and we will revise the manuscript to make this clearer.\n*    **Novel doubly robust (DR) identification in PCL**: Theorem 2.7 provides the first DR identification result in the PCL setting that avoids indicator functions and kernel smoothing. Prior DR methods [17, 18] are restricted to discrete or univariate continuous treatments, whereas our approach naturally extends to high-dimensional and continuous treatments. We highlight this advancement in Remarks 2.8 and 2.9.\n*    **Density-ratio-free DR formulation**: While avoiding density ratio estimation was introduced in [19], our integration of treatment bridge–based methods into a DR framework introduces substantive methodological differences. For example, [17] employs correction terms of the form $\\mathbb{E}[q(Z,A)h(W,A) I(A - a)]$, whereas our formulation uses $\\mathbb{E}[\\varphi(Z,A)h(W,A)| A = a],$ reflecting a distinct treatment bridge function $\\varphi$. .This yields a structurally new identification formula and results in the **first density-ratio-free DR PCL estimators**. Importantly, our kernel-based construction provides easily implementable closed-form solutions, which are nontrivial to derive as they require careful algebra leveraging properties of reproducing kernel Hilbert spaces (Appendix C.3).\n*    **Strong uniform convergence guarantees**: Our results (Theorems 4.1 and 4.2; details in Appendix E) establish **uniform consistency**, which is **strictly stronger than the pointwise convergence** typically used in earlier works. Achieving such rates in a doubly robust setting is technically demanding, requiring precise error control across multiple regression stages, including the novel correction term estimation.\n*   **Practical impact**: These methodological contributions translate into tangible performance gains. As shown in Section 5, our DR estimators (DRKPV and DRPMMR) substantially outperform outcome-only and treatment-only bridge methods across diverse PCL benchmarks, including challenging high-dimensional setups such as dSprite and the semi-synthetic JobCorps dataset. Furthermore, our method outperforms the previous doubly robust approach PKDR [18], which relies on kernel smoothing.\n\nWe will revise our manuscript to more clearly articulate our contributions, ensuring that the nontrivial nature of our DR development is evident. In particular, by utilizing the extra content page, some of the revisions we incorporate include:\n*    Abstract line 11 will be revised as: \"By using kernel mean embeddings, we propose the first density ratio-free doubly robust estimators for proxy causal learning (PCL), which have closed-form solutions and strong consistency guarantees.\"\n*    Line 102 will be revised as: \"We establish uniform consistency of our proposed estimators which is stronger than the pointwise convergence typical in doubly robust causal learning\"\n*    Revision to Remark (2.9): \"This structural difference enables applicability to high-dimensional treatments, where kernel-smoothing–based DR approach fails.\"\n\nTogether, these points demonstrate that our work provides a substantial and nontrivial advancement of the PCL literature. For further revisions on consistency section, see our corresponding response below. \n\n> **Questions**: 1-) The arrow from U to W should be solid, and the arrow from Y to W should be dashed.\n\nWe appreciate the opportunity to clarify our notation. While the counterexample you provided indeed violates Assumption (2.5), we kindly disagree that it aligns with the structure of Figure (1). In our notation, a dashed arrow indicates the possibility that either: (i) $U \\rightarrow W$, (ii) $W \\rightarrow U$, or (iii) $U$ and $W$ share a common ancestor. The example you provided implies no dependence between $U$ and $W$, which differs from our intended interpretation. We acknowledge that our explanation may not be sufficiently clear, and we will accordingly update our manuscript to address this. Specifically, we will add the clarification after line 28 in our updated manuscript.\n\nFurthermore, this convention is consistent with Table A.1 in Tchetgen Tchetgen et al. (2024, Statistical Science), which shows that both $U \\rightarrow W$ and $W \\rightarrow U$ are possible. \n\n> **Significance**: The high-dimensional extension is also promising, though the theoretical justification in this context requires further clarification.\n> **Questions**: 2-) What notion of convergence is used in the theorems? The paper includes high-dimensional simulations, but it’s unclear whether Theorems 4.1 and 4.2 are valid in such settings. How is the curse of dimensionality addressed? Please clarify the assumptions under which these theoretical results extend to high-dimensional cases.\n\nIn Theorems 4.1 and 4.2, we establish **uniform consistency** of our estimators DRKPV and DRPMMR. In fact, we obtain high-probability finite-sample bounds on the error in the supremum norm (see Appendix E.3, Theorems E.27 and E.28). Specifically, under smoothness assumptions on the regression stages (e.g. Assumptions E.5, E.6, E.11, E.12, E.17), our estimators converge to the true causal function at the rates we derive in Appendix E.3. These results rely on conditional mean embedding and kernel ridge regression (KRR) learning rates (see, e.g., [37, 52, 53]), which provide rigorous convergence guarantees.\n\nWe agree with the reviewer that the curse of dimensionality is a real concern. Kernel Ridge Regression (KRR) is known to achieve minimax-optimal rates in moderate dimension, but in very high-dimensional regimes, particularly when the input dimension $d$ grows with the sample size $n$, i.e. $d/n^\\beta \\rightarrow c$ for $\\beta \\in(0,1)$, the situation becomes more subtle [R1]. In particular:\n\n*    When functions lie in Sobolev classes, the achievable rate explicitly depends on both smoothness and dimension. Intuitively, the regression function must be “smooth enough” relative to the ambient dimension $d$ for KRR to remain consistent. More specifically, S. Fischer and I. Steinwart [51, Corollary 5 & 6] (for standard scalar KRR) and Li et al. [R2 Corollary 1 & 2] (for vector valued KRR) show that if the target function has smoothness $s$, then the optimal rate of convergence in $L_2$ norm is $O(n^{-\\frac{s}{s + d/2}})$ and that this rate is achieved by KRR. This rate clearly shows the curse of dimensionality where for a fixed $s$, the bound becomes vacuous as $d \\to +\\infty$.\n*    Donhauser et al. [R1] show that for rotationally invariant kernels, a polynomial approximation barrier arises: the learned function is effectively restricted to low-degree polynomials as $d$ grows, regardless of eigenvalue decay. This implies that consistency in high dimensions is limited unless additional structural assumptions are imposed.\n\nIn our work, we do not claim to resolve the curse of dimensionality. Instead, our contribution is to show that doubly robust PCL estimators can be constructed without density ratio estimation and kernel-smoothing, thereby extending applicability to continuous and high-dimensional treatments where prior DR methods [17, 18] fail. Our theorems remain valid in high-dimensional settings provided the smoothness and effective RKHS dimension assumptions hold, but we acknowledge that when $d$ increases with $n$, methods based on standard RBF or Matern kernels may indeed fail.\n\nWe will clarify this point in the updated version by extending the discussion in our consistency section.\n\nReferences\n\n[R1] Konstantin Donhauser, Mingqi Wu, Fanny Yang, (2021). \"How rotational invariance of common kernels prevents generalization in high dimensions\". PMLR\n\n[R2] Zhu Li, Dimitri Meunier, Mattes Mollenhauer, and Arthur Gretton. (2024) \"Towards optimal sobolev norm rates for the vector-valued regularized least-squares algorithm\". JMLR\n\n> **Questions**: 3-) In explaining Equation (2), it would be helpful to cite foundational work on stabilized weights for continuous treatments (e.g., Robins et al., 2000, Epidemiology). This would situate your methods more clearly in the literature and justify the weighting strategy.\n> **Originality**: Furthermore, some relevant literature is omitted, reducing the clarity of the paper's positioning relative to existing work.\n\nWe will incorporate the recommended literature in the updated manuscript. In particular, we will clarify the connection between the treatment bridge–based approach and inverse-probability-weighted causal learning methods. We would also be glad to consider any further references you believe would help strengthen the positioning of our work.\n> **Questions**: 5-) Typo in line 320. \n\nThank you for pointing out this typo; we have corrected it."}]

Author rebuttal:

Meta-review: This elegant paper proposed two doubly-robust kernel-based estimators for causal inference in the presence of unobserved confounders.  It builds on two different proxy approaches, when proxies are available for unobserved confounders, focusing on either the influence on treatment (propensity) or influence on outcome.  The paper does a nice job of reviewing the relevant recent literature, situating its contributions in the recent literature, and showing its advantages and limitations empirically and theoretically.

The paper should be proofread again carefully; for example, on line 159 the actual heading for section 2.3 says "Doubly bobust..."

The reviewers and authors had a robust discussion, so the authors should also please incorporate items from that discussion to improve clarity of the paper.
Decision: Accept (poster)