Paper deep dive
Robust Weighted Triangulation of Causal Effects Under Model Uncertainty
Rohit Bhattacharya, Ina Ocelli, Ted Westling
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 85%
Last extracted: 7/20/2026, 2:39:35 AM
Summary
The paper introduces a framework for robust causal effect triangulation that combines multiple candidate causal models under conditions of model uncertainty. It proposes a weighted triangulation functional using data-driven measures of model validity (testable implications) to assign smooth weights to model-specific estimates, avoiding explicit model selection and post-selection inference problems. Theoretical bounds are provided for the distance between the functional and the true causal effect, demonstrating robustness when at least one model is correct and testable.
Entities (8)
Relation Signals (7)
Triangulation Functional → addresses → Model Uncertainty
confidence 90% · When there is model uncertainty, analysts may seek to use estimates from multiple candidate models... Here, we develop a framework for causal effect triangulation
Triangulation Functional → combines → identified functionals
confidence 90% · We propose a triangulation functional that combines identified functionals from each model with data-driven measures of model validity.
Causal DAG → implies → Verma Constraint
confidence 85% · the observed distribution P(V) consists of not only conditional independences from d-separation, but also generalized equality constraints, known as Verma constraints
Triangulation Functional → uses → Gaussian Kernel
confidence 85% · We propose a smooth approximation... by replacing I(βk=0) with Gaussian kernels
Backdoor Adjustment → istypeof → Causal Model
confidence 80% · backdoor adjustment model... rely on qualitatively distinct assumptions
Frontdoor Model → istypeof → Causal Model
confidence 80% · One such alternative is the frontdoor model
Instrumental Variable → istypeof → Causal Model
confidence 80% · Alternatively, one may use instrumental variable (IV) models
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:A fundamental challenge in causal inference with observational data is correct specification of a causal model. When there is model uncertainty, analysts may seek to use estimates from multiple candidate models that rely on distinct, and possibly partially overlapping, sets of identifying assumptions to infer the causal effect, a process known as triangulation. Principled methods for triangulation, however, remain underdeveloped. Here, we develop a framework for causal effect triangulation that combines model testability methods from causal discovery with statistical inference methods from semiparametric theory, while avoiding explicit model selection and post-selection inference problems. We propose a triangulation functional that combines identified functionals from each model with data-driven measures of model validity. We provide a bound on the distance of the functional from the true causal effect along with conditions under which this distance can be taken to zero. Finally, we derive valid statistical inference for this functional. Our framework formalizes robustness under causal pluralism without requiring agreement across models or commitment to a single specification. We demonstrate its performance through simulations and an empirical application.
Tags
Links
- Source: https://arxiv.org/abs/2603.01119v2
- Canonical: https://arxiv.org/abs/2603.01119v2
Trouble viewing inline? Open PDF directly →
Full Text
76,038 characters extracted from source content.
Expand or collapse full text
Robust Weighted Triangulation of Causal Effects Under Model Uncertainty Rohit Bhattacharya Ina Ocelli Ted Westling Dept. of Mathematics and Statistics University of Massachusetts Amherst Abstract A fundamental challenge in causal inference with observational data is correct specification of a causal model. When there is model uncertainty, analysts may seek to use estimates from multiple candidate models that rely on distinct, and possibly partially overlapping, sets of identifying assumptions to infer the causal effect, a process known as triangulation. Principled methods for triangulation, however, remain underdeveloped. Here, we develop a framework for causal effect triangulation that combines model testability methods from causal discovery with statistical inference methods from semiparametric theory, while avoiding explicit model selection and post-selection inference problems. We propose a triangulation functional that combines identified functionals from each model with data-driven measures of model validity. We provide a bound on the distance of the functional from the true causal effect along with conditions under which this distance can be taken to zero. Finally, we derive valid statistical inference for this functional. Our framework formalizes robustness under causal pluralism without requiring agreement across models or commitment to a single specification. We demonstrate its performance through simulations and an empirical application. 1 Introduction Causal inference with observational data is rarely conducted under a single, indisputable set of assumptions. For a single observed dataset, there may be multiple plausible causal models that identify the target causal parameter under different sets of assumptions, some of which may be untestable. A common response to model uncertainty is to combine evidence from each of these models with the hope of reaching more robust conclusions than reliance on just a single one, a process commonly referred to as triangulation [thurmond2001point, farmer2006developing, lawlor2016triangulation]. Despite its intuitive appeal, triangulated effect estimation remains underdeveloped. A prominent line of work formalizes triangulation through evidence factors, which combine p-values across multiple analyses to test causal hypotheses [rosenbaum2010evidence, rosenbaum2011some, karmakar2019integrating]. While powerful, these approaches are fundamentally geared toward hypothesis testing rather than estimation, and rely on assumptions that are often difficult to justify in observational settings. In particular, they require that different analyses do not share sources of bias and that the joint distribution of p-values satisfies certain stochastic dominance properties under the null. yangstatistical leveraged the joint convergence properties of semiparametric estimators to design a robust test of the causal null hypothesis that remains valid as long as at least one model is correct. However, yangstatistical also focused on hypothesis testing rather than effect estimation. In this paper, we develop a general framework for triangulating causal effect estimates across multiple candidate models without requiring explicit model selection or correctness of a plurality of models, as required by voting-based triangulation procedures. Our approach instead assigns data-driven weights to each model based on testable implications of its identifying assumptions, such as conditional independence or generalized equality constraints. These weights are used to form a smooth aggregation of model-specific effect estimates, yielding an estimator that is consistent for a weighted triangulation functional. We show that the absolute difference between this functional and the true causal parameter can be bounded as a function of (i) the maximal bias among incorrect models and (i) our ability to separate correct and incorrect models using observed data. We use this bound to provide conditions under which the bias of the triangulation functional is small, yielding a form of robustness to misspecification of some causal models through triangulation. In particular, we show the bias can be controlled if at least one candidate model is correct and testable from observed data. We demonstrate the effectiveness of the proposed method through simulations and an empirical application. Contributions. Our contributions can be viewed from two complementary perspectives. First, we advance the triangulation literature by providing a principled and quantitative method for combining causal effect estimates that achieves robustness to model misspecification by leveraging testable implications in the observed data. Second, we contribute to the literature on post-selection inference, which studies valid inference after data-driven model selection. Rather than selecting a single model, our approach avoids explicit selection by using smooth weights, thereby mitigating post-selection bias while still incorporating information from model diagnostics. Unlike recent work in causal discovery [gradu2025valid, chang2026post], which focuses on learning the full causal graph, our framework targets a more modest but practically important goal: testing a minimal set of assumptions required to identify a causal parameter and combining estimators of the identified functionals. Related work. The methods in rakshit2025adaptive, kang2016instrumental, sun2021multiply, yao2024deciphering allow for some models to be misspecified, but require that a plurality of models be correctly specified. This can be a strong assumption in observational settings, where several models may fail for similar reasons. Moreover, these methods are typically tailored to specific classes of causal models, such as proxy-variable approaches or instrumental variables. Bayesian model averaging approaches [horii2021bayesian, steiner2025bayesian] provide another avenue for combining estimates, but often rely on parametric assumptions such as linearity or are similarly restricted to specific model classes. 2 Causal Graph Preliminaries ZZAAMMYYCCUU(a) (V∪U)G(V∪ U)ZZAAmmYYCCUU(b) (V∪U)do(m)G(V∪ U)_do(m) Figure 1: (a) A hidden variable causal DAG. (b) Graph representing intervention on M from which we can read the Verma constraint Z⟂Y∣CZ \!\!\! Y C in P(V)/P(M∣A,Z,C)P(V)/P(M A,Z,C). AAYYC4,C5\C_4,C_5\C1,C2,C3\C_1,C_2,C_3\U1U_1U2U_2(a) M-bias exampleAAYYCCUU(b) Backdoor modelAAMMYYUU(c) Frontdoor modelZZAAYYUU(d) IV model Figure 2: Causal DAGs used for motivating robust triangulation. (a) A causal DAG where uncertainty about what variables to adjust for may yield M-bias. (b, c, d) Under model uncertainty, analysts may wish to triangulate effects from each of these models that rely on qualitatively distinct assumptions—violations of these assumptions are shown via blue dashed edges. Although the candidate causal models need not be graphical, we frame our discussion with causal directed acyclic graphs (causal DAGs), as they provide a transparent way to represent identifying assumptions and their testable implications. The causal model of a DAG (V)G(V) can be understood as distributions generated by a system of structural equations equipped with the do(⋅)do(·) operator [pearl2009causality]. Specifically, for each variable Vi∈V_i∈ V, there is an equation of the form Vi←fi(pa(Vi),ϵi),V_i← f_i (pa_G(V_i), _i ), where pa(Vi)pa_G(V_i) denotes a set of values for the parents of ViV_i in G and ϵi _i is an exogenous (latent) error term. A joint distribution P(V)P(V) induced by such a system is said to be Markov with respect to G, meaning that it factorizes as P(V)=∏Vi∈VP(Vi∣Pa(Vi))P(V)= _V_i∈ VP (V_i _G(V_i) ). Equivalently, P satisfies the global Markov property with respect to G stated in terms of the well-known d-separation criterion: X⟂d-sepY∣Z⟹X⟂Y∣Z in PX \!\!\! _d-sepY Z X \!\!\! Y Z in P [pearl2009causality]. The distribution P is considered faithful to G if the converse is also true, i.e., X⟂d-sepY∣Z⇔X⟂Y∣Z in PX \!\!\! _d-sepY Z X \!\!\! Y Z in P. In fully observed causal DAG models, counterfactual distributions arising from interventions on a set A⊂VA⊂ V, written as P(V∖A∣do(a))P(V A (a)), are identified via a truncated factorization known as the g-formula [robins1986new, pearl2009causality]: P(V∖A∣do(a))=∏Vi∈V∖AP(Vi∣Pa(Vi))|A=a. P(V A (a))= _V_i∈ V AP (V_i _G(V_i) ) |_A=a. (1) When some variables U are unobserved, identification theory becomes more complex (as in Scenario 2 of Section 3). In addition, the observed distribution P(V)P(V) consists of not only conditional independences from d-separation, but also generalized equality constraints, known as Verma constraints [robins1986new, verma1990equivalence]. Verma constraints are obtained via d-separation in conditional DAGs corresponding to identifiable post-intervention distributions. As an example, consider the DAG (V∪U)G(V∪ U) in Figure 1(a). By d-separation, G implies Z⟂M∣A,CZ \!\!\! M A,C. However, there is no set X⊂V∖Z,YX⊂ V \Z,Y\ such that Z⟂d-sepY∣XZ \!\!\! _d-sepY X. That is, there is no ordinary conditional independence between Z and Y in the observed joint P(V)P(V). However, it is well known that there exists a Verma constraint Z⟂Y∣CZ \!\!\! Y C in the Markov kernel P(V)/P(M∣A,Z,C)P(V)/P(M A,Z,C) [verma1990equivalence]. Under a causal interpretation of (V∪U)G(V∪ U), for any fixed value of m, P(V)/P(M=m∣A,Z,C)|M=m=P(Z,A,C,Y∣do(m))P(V)/P(M=m A,Z,C)|_M=m=P(Z,A,C,Y (m)) by the g-formula in (1). That is, the kernel P(V)/P(M∣A,Z,C)P(V)/P(M A,Z,C) factorizes according to the conditional DAG in Figure 1(b), where M is now a fixed node with all incoming edges removed, and from which the Verma constraint can be read using d-separation—notice that Z and Y are d-separated given C in (V∪U)do(m)G(V∪ U)_do(m). When the measures of model validity used rely on Verma constraints, we will require P(V)P(V) to be Verma constraint faithful to (V∪U)G(V∪ U). That is, X⟂Y∣ZX \!\!\! Y Z in P(V)/P(A|B)P(V)/P(A|B) implies that edges in G are such that P(V∖A∣do(a))=P(V)/P(A|B)|A=aP(V A (a))=P(V)/P(A|B)|_A=a and X⟂d-sepY∣ZX \!\!\! _d-sepY Z in do(a)G_do(a). 3 Motivating Scenarios Suppose an analyst is interested in estimating the average causal effect θ≡[Y∣do(A=1)]−[Y∣do(A=0)]θ [Y do(A=1)]-E[Y do(A=0)]. In this section, we present two classic scenarios where proper triangulation of effect estimates from multiple plausible causal models would lead to more robust causal inference. Scenario 1. Adjust for all pre-treatment covariates or not? When a set of pre-treatment covariates L blocks all backdoor paths between A and Y we obtain a well-known instance of the g-formula, known as the backdoor formula: θ=∑lP(l)×([Y∣A=1,l]−[Y∣A=0,l]). \!\!θ= _lP(l)× (E[Y A=1,l]-E[Y A=0,l] ). (2) A longstanding debate in causal inference is whether analysts should include all observed pre-treatment covariates in the backdoor formula (2). Some researchers [rosenbaum2002observational, ding2015adjust, rubin2009should] argue that conditional ignorability (blockage of all backdoor paths) is likely to hold only when adjusting for a large set of variables, and thus advocate adjusting for all pre-treatment covariates C associated with the treatment and outcome. Other researchers emphasize that adjustment for all such variables can induce collider bias, such as M-bias [pearl2009remarks, shrier2009remarks, sjolander2009propensity]. Thus, when the true DAG is unknown—as is typical in observational studies—the analyst may face significant uncertainty about whether excluding certain covariates risks confounding bias or including them risks collider bias. This setting naturally motivates triangulation across multiple adjustment strategies. As a concrete example, let Figure 2(a) be the true (but unknown) DAG. Here, adjusting for all observed pre-treatment covariates C yields invalid effect estimates due to the colliders at C4C_4 and C5C_5. Due to model uncertainty, the analyst might partition the observed covariates into two groups: those that certainly do not give M-bias and are essential for adjustment, and those about which there is some uncertainty. Say they correctly deem C1,C2,C3\C_1,C_2,C_3\ to be essential, but are uncertain about C4,C5\C_4,C_5\. They may then wish to triangulate estimates from multiple overlapping adjustment sets, e.g., C1,C2,C3\C_1,C_2,C_3\ that excludes all uncertain covariates, C1,C2,C3,C4\C_1,C_2,C_3,C_4\ that includes some but not others, and C1,C2,C3,C4,C5\C_1,C_2,C_3,C_4,C_5\ that includes all pre-treatment covariates. In this example, only the first set leads to valid effect estimates, underscoring the importance of a triangulation procedure that is robust to incorrect identifying assumptions in multiple models as long as at least one is correct. Scenario 2. Use backdoor, frontdoor, or an instrument? In many observational settings, there is uncertainty about whether unmeasured confounding between A and Y is present. A natural response to this is to triangulate estimates obtained from a backdoor adjustment model (Figure 2(b)) with those from alternative identification strategies that explicitly allow for unmeasured A–Y confounding but impose different structural assumptions. One such alternative is the frontdoor model (Figure 2(c)). The frontdoor model permits unmeasured confounding between A and Y, but assumes (i) that A affects Y only through a mediator set M (i.e., no direct A→YA→ Y edge), and (i) that the unmeasured variable U does not confound the A–M or M–Y relationships (i.e., no U→MU→ M edge). Under these assumptions, the causal effect is identified by the frontdoor formula [pearl1995causal]. A conditional version that adjusts for baseline covariates L is θ=∑m,l θ= _m,l \ (∑a′P(a′,l)⋅[Y∣a′,m,l]) ( _a P(a ,l)·E[Y a ,m,l] ) ⋅(P(m∣A=1,l)−P(m∣A=0,l)). · (P(m A=1,l)-P(m A=0,l) ) \. (3) Alternatively, one may use instrumental variable (IV) models (Figure 2(d)). An instrument Z typically satisfies three structural assumptions: (i) conditional independence from the unmeasured confounder U (no U→ZU→ Z edge), (i) an exclusion restriction (no Z→YZ→ Y edge), and (i) instrument relevance (Z→AZ→ A exists). Under additional non-graphical assumptions, such as effect homogeneity, the causal effect is point-identified as [angrist1996identification]: θ=∑lP(l)⋅([Y∣Z=1,l]−[Y∣Z=0,l])∑lP(l)⋅([A∣Z=1,l]−[A∣Z=0,l]). θ= _lP(l)· (E[Y Z=1,l]-E[Y Z=0,l]) _lP(l)· (E[A Z=1,l]-E[A Z=0,l] ). (4) Unlike Scenario 1, this setting involves three qualitatively distinct causal models. Each has their own unique drawbacks, but may imply testable restrictions in P(V)P(V) under some assumptions of faithfulness and causal ordering of variables [entner2013data, bhattacharya2022testability]. A natural approach then would be to test such implications, retain non-rejected models, and aggregate their estimates. However, this two-step strategy introduces a post-selection inference problem, particularly when models share data, have overlapping adjustment sets, or shared sources of bias. Further, standard remedies for this, such as sample splitting [hansen2000sample, newey2018cross], can substantially reduce effective sample size when multiple models are considered. These scenarios and associated challenges motivate our new triangulation functional and inference procedure. 4 New Triangulation Method We now present the general setup that we consider. Let θ denote the causal parameter of interest, such as the average or conditional average causal effect. We assume the observed data consists of n IID realizations O1,…,OnO_1,…,O_n drawn from some unknown distribution P. For a finite integer K>1K>1, let ℳ1…,ℳK M_1…, M_K denote the different candidate causal models, and let ψk _k (we suppress the dependence on P for these parameters and ones that follow for notational brevity) denote the identifying functional of ℳk M_k. That is, if the assumptions of ℳk M_k are true, then θ=ψkθ= _k. Let β1,…,βK _1,…, _K denote observed data parameters and A be a set of assumptions such that if A is true and βk=0 _k=0, then the identifying assumptions of ℳk M_k are true. We do not require, however, that if A holds and βk≠0 _k =0, then the assumptions of ℳk M_k must be incorrect; we allow this scenario to also correspond to the model ℳk M_k being untestable using the observed data. The role of each βk _k is to encode a testable implication of the identifying assumptions of model ℳk M_k. Examples of such testable implications from the causal discovery literature include conditional independence constraints for backdoor models [entner2013data], Verma constraints for frontdoor models [bhattacharya2022testability], and tetrad constraints for proxy-based models [xie2024automating]. That is, each βk _k is an observed data parameter constructed so that it is zero if and only if a particular equality constraint, such as the aforementioned ones, holds in the observed distribution P. More specifically, βk _k may be a regression coefficient in a parametric model, a log-odds ratio in a semiparametric model [chen2007semiparametric], or a generalized covariance measure [shah2020hardness, he2025on, bergen2026the]. Under assumptions A (typically including a faithfulness condition and partial knowledge of temporal ordering of variables), βk=0 _k=0 implies the absence of specific edges in G, and absence of these edges is sufficient for the identifying assumptions of ℳk M_k to hold. Thus, whether βk=0 _k=0 provides empirical evidence about the validity of ℳk M_k. 4.1 Triangulation functional The above setup motivates the following naive triangulation functional: ψnaive=∑k=1K(βk=0)⋅ψk∑j=1K(βj=0) _naive= _k=1^KI( _k=0)· _k _j=1^KI( _j=0), where (⋅)I(·) is the indicator function. That is, ψnaive _naive is an average over all ψk _k derived from models ℳk M_k for which βk=0 _k=0. If assumptions A used to test model correctness hold, ψnaive _naive is causally robust: ψnaive=θ _naive=θ if at least one model is correct and testable, as it reduces to an average of only those ψk _k such that ψk=θ _k=θ. We consider this to be naive as filtering based on (βk=0)I( _k=0) is unlikely to succeed for any model ℳk M_k when βk _k is unknown and must be estimated from finite samples. Instead, we propose a smooth approximation of the naive functional by replacing (βk=0)I( _k=0) with Gaussian kernels of the form δa(βk)=1|a|πe−(βk/a)2 _a( _k)= 1|a| πe^- ( _k/a )^2, where a>0a>0 is a constant that controls the sharpness of the approximation. As a→0a→ 0, the function becomes concentrated around βk=0 _k=0. Our proposed triangulation functional is then given by ψ=∑k=1Kwkψk,wherewk=δa(βk)∑j=1Kδa(βj). ψ= _k=1^Kw_k _k, w_k= _a( _k) _j=1^K _a( _j). (5) This functional exhibits an approximate causal robustness property due to smoothing via the kernels, as stated below. Theorem 1. Suppose assumptions A used to justify the tests of model correctness hold, and let ⊆1,…,KC \1,…c,K\ and ℐ=1,…,K∖I=\1,…c,K\ be the subsets of indices k such that model ℳkM_k is correct and incorrect, respectively. Then |ψ−θ|≤maxk|ψk−θ|1+Da|ψ-θ|≤ _k| _k-θ|1+D_a (6) where Da=(∑k∈δa(βk))/(∑k∈ℐδa(βk))D_a= ( _k∈ C _a( _k) )/ ( _k∈ I _a( _k) ). Further, if at least one model ℳk M_k is correct and testable using the observed data, then Da≥eε2/a2/|ℐ|D_a≥ e ^2/a^2/| I|, where ε=mink∈ℐ|βk| = _k∈ I| _k|. The proof is in Appendix A. Theorem 1 demonstrates that the absolute difference between the triangulation functional ψ and the true causal parameter of interest θ is bounded by the maximal bias of the functionals ψk _k from the incorrect causal models divided by a “discrimination factor" 1+Da1+D_a. The discrimination factor depends both on the weighting function chosen and the true values of the testing functionals β1,…,βK _1,…c, _K. Roughly speaking, when the βk _k’s from correct causal models are small in magnitude compared to those from incorrect causal models on the scale determined by δa _a, then the discrimination factor is large, and the absolute difference between ψ and θ is reduced. We note that the second statement of Theorem 1 implies that if at least one model ℳk M_k is both correct and testable using the observed data and all testing functionals βk _k in incorrect models are non-zero, then lima→0Da=∞ _a→ 0D_a=∞, so that lima→0ψ=θ _a→ 0ψ=θ. To make the robustness property more concrete, consider setting the Gaussian kernel parameter to a=0.1a=0.1 and suppose there is one correct model =1C=\1\ and two incorrect models ℐ=2,3I=\2,3\. Let ε=mink∈ℐ|βk|=0.2 = _k | _k|=0.2, indicating that violations of the incorrect models are only weakly detectable since the value of βk _k is close to 0 even when ℳk M_k is incorrect. From Theorem 1, this yields 1+Da≥1+e(0.2/0.1)2/2≈301+D_a≥ 1+e^(0.2/0.1)^2/2≈ 30, implying that the bias of the triangulation functional is at most the maximal bias among incorrect models divided by 3030. While smaller values of a reduce the absolute difference between the triangulated functional ψ and the target causal effect θ, in practice we cannot set a too small because we only have estimates of ψk _k and βk _k. Thus, we must balance the variability of these estimators against the bias induced by a>0a>0. We return to this discussion and provide concrete recommendations for setting a in Section 4.3. We also note that if no candidate model is both correct and testable, then the triangulation functional does not converge to θ as a→0a→ 0. Fortunately, this behavior is straightforward to diagnose in practice: the estimated values of ψ often diverge due to ∑j=1Kδa(βj)≈0 _j=1^K _a( _j)≈ 0. The normalizing term, however, also poses challenges in finite samples even when a correct and testable model ℳk M_k exists in theory, e.g., due to sampling variability. We briefly address this issue of numerical stability in the following subsection on inference. 4.2 Inference Procedure We now discuss our approach to estimation and inference for the triangulation functional ψ. Suppose that for each k, we have estimators ψk,n _k,n and βk,n _k,n of each identifying functional ψk _k and each measure of model correctness βk _k, respectively. Our triangulation estimator is then ψn=∑k=1Kwk,nψk,n,for wk,n=δa(βk,n)λn+∑jδa(βj,n), \!\! _n= _k=1^Kw_k,n _k,n,\ for \ w_k,n= _a( _k,n) _n+ _j _a( _j,n), (7) where λn>0 _n>0 is an additional term used to ensure numerical stability of the weights wk,nw_k,n even when the normalizing function in the denominator is close to zero. To avoid affecting first-order asymptotics of the estimator, we require λn=o(n−1/2) _n=o(n^-1/2). A simple choice is to fix λn=1/n _n=1/n. Note that λn _n in the triangulation functional is optional. In practice, we find in our simulations and data application that we do not need it and that using the well-known log-sum-exp trick [blanchard2021accurately] to compute the normalizing function is sufficient to prevent numerical instability. However, we keep it in our exposition for the sake of completeness. Since ψ1,n,…,ψK,n,β1,n,…,βK,n _1,n,…, _K,n, _1,n,…, _K,n are constructed using a common observed dataset, they are typically dependent, even asymptotically. We suppose that under a set of statistical regularity conditions S (e.g., rates of convergence and complexity constraints on the nuisance estimators), we can construct asymptotically linear estimators of both βk _k and ψk _k with influence functions ϕβk _ _k and ϕψk _ _k, respectively, for each k. That is, under S we have βk,n−βk=1n∑i=1nϕβk(Oi)+op(n−1/2) _k,n- _k= 1n _i=1^n _ _k(O_i)+o_p(n^-1/2) and ψk,n−ψk=1n∑i=1nϕψk(Oi)+op(n−1/2) _k,n- _k= 1n _i=1^n _ _k(O_i)+o_p(n^-1/2). Here ϕβk _ _k and ϕψk _ _k are assumed to satisfy [ϕβk]=[ϕψk]=0E[ _ _k]=E[ _ _k]=0, [ϕβk2]<∞E[φ^2_ _k]<∞, and [ϕψk2]<∞E[φ^2_ _k]<∞. These estimators could be parametric or semiparametric in nature; we will discuss specific approaches to estimation more below. Asymptotic linearity of any finite collection of estimators implies joint convergence to a multivariate normal distribution [van2000asymptotic]. Let κ=[β1,…,βK,ψ1,…,ψK]κ=[ _1,…, _K, _1,…, _K] and κn _n denote the corresponding vector of estimates. Then, n1/2(κn−κ)→N(0,Σ)n^1/2( _n-κ) dN(0, ), where Σ is a 2K×2K2K× 2K covariance matrix of the influence functions of each ψk _k and βk _k. Each entry of Σ is of the form [ϕβj,ϕβk],[ϕψj,ϕψk]E[ _ _j, _ _k],E[ _ _j, _ _k], or [ϕβj,ϕψk]E[ _ _j, _ _k]. Applying the delta method then gives the asymptotic distribution of the triangulation estimator as, n1/2(ψn−ψ)→N(0,γTΣγ), n^1/2( _n-ψ) dN(0,γ^T γ), (8) where γ is a vector of partial derivatives defined as, γ=[∂ψ∂β1,…,∂ψ∂βK,∂ψ∂ψ1,…,∂ψ∂ψK]. γ= [\ ∂ψ∂ _1,…, ∂ψ∂ _K, ∂ψ∂ _1,…, ∂ψ∂ _K\ ]. (9) The following lemma also gives us closed form expressions for computing the vector γ. The proof is in Appendix B. Lemma 1. The partial derivatives of ψ are ∂ψ∂ψk=wk and ∂ψ∂βk=2βkwka2(ψ−ψk). ∂ψ∂ _k=w_k and ∂ψ∂ _k= 2 _kw_ka^2 (ψ- _k ). We now suggest three approaches to constructing asymptotically linear estimators and associated approaches to valid inference. Our preferred strategy is influence function-based estimation. In this approach, the influence functions ϕβk _ _k and ϕψk _ _k are derived explicitly and used to construct βk,n _k,n and ψk,n _k,n using, e.g., the one-step construction [bickel1982adaptive] or estimating equations [chernozhukov2018double]. A benefit to this approach is that machine learning estimators can be used for nuisance functions, and asymptotic linearity follows under sufficient rates of convergence of these nuisance estimators. When the influence functions are available in closed form—several of which have been derived in the context of causal graphs jung2021estimating, bhattacharya2022semiparametric, guo2024average—we can construct consistent estimators ϕβk,n _ _k,n and ϕψk,n _ _k,n of the influence functions. A consistent estimator Σn _n of Σ is then simply the sample covariance matrix of the influence function estimators. An estimator γn _n of γ can be obtained by plugging in ψk,n _k,n and βk,n _k,n into the form of γ provided in Lemma 1. A variance estimator of ψn _n is then given by var^(ψn)=γnTΣnγn/n var( _n)= _n^T _n _n/n, which can be used to construct Wald-type confidence intervals for ψ. For some functionals ψk _k or βk _k, explicit influence functions might not be available in the current literature. In this case, one option is to use plug-in estimators ψk,n _k,n and βk,n _k,n. However, plug-in estimators are often only asymptotically linear when correctly-specified parametric models are employed for nuisance estimators. In this case, the empirical bootstrap yields asymptotically valid inference under mild smoothness conditions on the parametric models [efron1994introduction]. Specifically, on each bootstrap sample, the parametric nuisance estimators are re-estimated and plugged in to obtain bootstrap estimates ψk,n∗ _k,n^* and βk,n∗ _k,n^*. These estimates are then substituted as in (7) to obtain a bootstrap estimate ψn∗ _n^*. The distribution of bootstrap estimates can then be used to construct a confidence interval in the standard ways. Finally, if plug-in estimators with nuisance estimators based on machine learning are used, the resulting estimators may not be asymptotically linear due to excess bias, and hence the empirical bootstrap is not guaranteed to work. In this case, subsampling may still be asymptotically valid [politis2001asymptotic]. Subsampling involves drawing b random subsamples of size m<nm<n from the full dataset, computing estimates ψm(1),…,ψm(b) _m^(1),…, _m^(b) for each subsample, and constructing a (1−α)(1-α) confidence interval based on the α/2α/2 and 1−α/21-α/2 quantiles of these estimates. A sufficient condition for this procedure to yield valid inference is that τn(ψn−ψ) _n( _n-ψ) converges to a non-degenerate distribution, which can be true even when the estimators are not asymptotically linear. We summarize our recommendations below. The order in the if-elif-else logic reflects our preferred inference strategies. Inference Procedure for Triangulation 1. Construct estimators ψk,n _k,n and βk,n _k,n for each model ℳk M_k and obtain the point estimate ψn _n using (7). 2. If every ψk,n _k,n and βk,n _k,n is influence-function-based, let SE^(ψn)=γn⊤Σnγn/n SE( _n)= _n _n _n/n and construct a (1−α)(1-α) CI as ψn±z1−α/2SE^(ψn) _n± z_1-α/2\, SE( _n). 3. Else, if ψk,n _k,n and βk,n _k,n are plug-in estimators based on parametric nuisance models: Construct a (1−α)(1-α) CI as [2ψn−q1−α/2, 2ψn−qα/2][2 _n-q_1-α/2,\,2 _n-q_α/2], where qpq_p is the pthp^th empirical quantile of estimates ψn(1),…,ψn(b)\ _n^(1),…, _n^(b)\ obtained via empirical bootstrap and re-fitting of nuisance models. 4. Else, construct a (1−α)(1-α) CI as [2ψn−q1−α/2, 2ψn−qα/2][2 _n-q_1-α/2,\,2 _n-q_α/2], based on the empirical quantiles of estimates ψm(1),…,ψm(b)\ _m^(1),…, _m^(b)\ obtained via subsampling with m=n4/5m=n^4/5 and re-fitting of nuisance models. 4.3 Setting the kernel bandwidth Finally, we address setting the kernel bandwidth parameter a. Our recommendations in this section apply when estimators that converge at the rate n−1/2n^-1/2, which includes influence function-based estimators and plug-in estimators based on parametric nuisance models (options 2 and 3 of the display in Section 4.2). Setting a when using estimators with slowers rates of convergence, such as many plug-in estimators based on flexible nuisance models, is an interesting problem left to future work. We recommend choosing a such that (1) n1/2a→∞n^1/2a→∞ and (2) alog(n)→0a (n)→ 0. Examples of this include setting a=n−1/3a=n^-1/3 and a=1/log(n)a=1/ (n). In our numerical experiments and data application we choose a=n−1/3a=n^-1/3, and also test this setting against other kernel bandwidths that lie outside of our recommended range. Our recommendation is based on the following theoretical analysis. From Theorem 1 and Section 4.2, the bias of the triangulation estimator ψn _n for the true parameter θ decays at a rate proportional to e−1/a2+oP(n−1/2)e^-1/a^2+o_P(n^-1/2), while its standard deviation (SD) goes to zero at rate n−1/2n^-1/2 when using influence function-based or parametric estimators. If alog(n)→0a (n)→ 0 (i.e., a goes to zero faster than 1/log(n)1/ (n)), then e−1/a2=o(1/n)e^-1/a^2=o(1/ n), so that the bias goes to zero faster than the SD. Hence, as long as at least one model is correct and testable and regularity conditions for nuisance estimation hold, then our methods produce asymptotically valid inference for the true causal effect if alog(n)→0a (n)→ 0. Now let ′ C denote the set of indices of all models ℳk M_k that are both correct and testable. When the estimators βk,n _k,n are asymptotically linear, then nβk,n n _k,n is asymptotically normal if model k is correct and testable, and cnβk,n→p∞c_n _k,n _p∞ for any sequence cn→∞c_n→∞ otherwise. This implies the following: (i) If a goes to zero faster than than n−1/2n^-1/2, then asymptotically ψn _n is equal to ψk,n _k,n for some k randomly selected from ′ C ; (i) If a goes to zero at the rate n−1/2n^-1/2, then asymptotically ψn _n is a random weighted combination of all ψk,n _k,n for k∈′k∈ C , with the joint distribution of weights depending on the asymptotic covariance of nβk,n n _k,n for k∈′k∈ C ; (i) If a goes to zero slower than n−1/2n^-1/2, then asymptotically ψn _n is an equally-weighted average of all ψk,n _k,n for k∈′k∈ C . In our view, scenario (i) is the most desirable behavior. Asymptotically, scenario (i) will result in focusing exclusively on the model whose testing functional estimator is closest to zero in any particular sample, which is overly narrow. The asymptotic covariance of the testing functionals is not necessarily related to the asymptotic variance of ψk,n _k,n, so there is not a good reason to weight estimators as in scenario (i). Thus, (i) arguably reflects the most reasonable behavior. It follows that we recommend choosing a such that n1/2a→∞n^1/2a→∞ (i.e., a goes to zero slower than n−1/2n^-1/2). 5 Applying The Triangulation Procedure In Practice AAZZYYC45C_45C123C_123U1U_1U2U_2AAZZMMYYCCUU Figure 3: Causal DAG used in our simulations to demonstrate robustness to (a) M-bias in some adjustment sets, and (b) misspecification of backdoor, frontdoor, or IV models. We now illustrate practical applications of our general method, along with numerical experiments and an empirical data application. The first subsection outlines a general framework for triangulating causal effect estimates derived from backdoor models with different candidate adjustment sets. The second subsection describes triangulation with backdoor, frontdoor model, and IV models. Explicit descriptions of the data generating processes and estimators are in Appendix D and E respectively, with proofs in Appendix C. 5.1 Multiple Candidate Adjustment Sets We first establish a set of assumptions A under which we can use observed data parameters βk _k to test the validity of each proposed adjustment set corresponding to different causal models ℳk M_k. Let C denote all pre-treatment covariates under consideration for adjustment and Z denote an “anchor” variable, such that the following assumptions hold: 1:P is faithful wrt a causal DAG (V∪U) A_1:P is faithful wrt a causal DAG G(V∪ U) 2: satisfies the causal ordering Z,C<A<Y A_2:G satisfies the causal ordering \Z,C\<A<Y 3:Z→A→Y exists in A_3:Z→ A→ Y exists in G (10) Note that 1 A_1 above is ordinary faithfulness as we do not rely on Verma constraints for this particular application, 2 A_2 simply imposes some weak background knowledge on the causal ordering of the variables, and 3 A_3 is a relevance assumption that ensures the testable implications are not trivial and actually rule out non-identifying edges in G.111Since A→YA→ Y exists, the sharp causal null is assumed false. For robust causal null hypothesis tests, see yangstatistical. Proposition 1. Let ℳk M_k be a backdoor model that proposes W⊆CW C as an adjustment set and βk _k be any observed data parameter such that βk=0⇔Y⟂Z∣A,W _k=0 Y \!\!\! Z A,W. Under assumptions A in (10), βk=0 _k=0 implies ℳk M_k is correct, i.e., ψk=θ _k=θ, where ψk _k is the backdoor formula (2) with L=WL=W. A natural choice for βk _k above is the log-odds ratio. For distinct sets of variables A,B,CA,B,C and a choice of reference values a0,b0a_0,b_0, the odds ratio function OR(A,B∣C)OR(A,B C) is a non-parametric measure of association given by [chen2007semiparametric], OR(a,b∣c)=P(a∣b,c)P(a0∣b,c)×P(a0∣b0,c)P(a∣b0,c), (a,b c)= P(a b,c)P(a_0 b,c)× P(a_0 b_0,c)P(a b_0,c), and log(OR(A,B∣C))=0⇔A⟂B∣C (OR(A,B C))=0 A \!\!\! B C. Though the log-odds ratio is a function, it is most commonly treated as a single scalar parameter in a semiparametric model; see for e.g., [chen2007semiparametric, tchetgen2010doubly, malinsky2019potential]. We adopt this setup, computing log-odds ratios of the form log(OR(Y,Z∣A,W)) (OR(Y,Z A,W)) for different choices of adjustment sets W as given by ℳ1,…,ℳk M_1,…, M_k, and use these as our βk _k in the triangulation functional (5). We return to the M-bias scenario in Section 3 as a concrete application. The candidate adjustment sets are ℳ1:C1,C2,C3 M_1:\C_1,C_2,C_3\, ℳ2:C1,C2,C3,C4 M_2:\C_1,C_2,C_3,C_4\, and ℳ3:C1,C2,C3,C4,C5 M_3:\C_1,C_2,C_3,C_4,C_5\, where only ℳ1 M_1 is correct. Figure 3(a) shows a DAG with an anchor variable Z satisfying assumptions A in (10) with or without the red dashed edge. Note that choices for anchor variables in the causal discovery literature are often candidate IVs, such as Z in Figure 3(a). While one could include IV-based estimates in the triangulation functional, when a valid adjustment set exists, backdoor estimators are more efficient and avoid homogeneity assumptions. By Proposition 1, β1=log(OR(Y,Z∣A,C123))=0 _1= (OR(Y,Z A,C_123))=0, while β2 _2 and β3 _3 are non-zero. We evaluate two settings, one where ε=mink∈ℐ=|βk|=0.71 = _k∈ I=| _k|=0.71, and another where ε=0.36 =0.36. We achieve the latter by excluding U1→ZU_1→ Z in Figure 3(a). Per Theorem 1, the second setting is more challenging for our triangulation estimator, as the absolute bias |ψ−θ||ψ-θ| is higher, since it is harder to detect the incorrect models. Here we use influence function-based estimators of the triangulation functional. As shown in Figure 4 (top row), the estimator remains causally robust in both cases, despite the plurality of models being incorrect. We require more samples, however, to obtain reliable estimates when ε=0.36 =0.36, as predicted by our theory. To test our variance computation proposal and Wald-type CI construction, we also compute the coverage of the estimator for both ψ and the target parameter θ. We report coverage for both, as these parameters are different in general, albeit with a small bounded difference. At n=5000n=5000, the estimator achieves the nominal coverage of 95%95\% for both ψ and θ in both settings. 5.2 Backdoor, frontdoor, and IV Figure 4: Point estimates averaged over 200200 trials; shaded bands correspond to 2.52.5 and 97.597.5 percentiles of the estimates. Similar to the previous subsection, we first establish assumptions under which we can use observed data to test each ℳk M_k, where ℳk M_k could be a backdoor, frontdoor, or IV model. As before, let C denote pre-treatment covariates used for adjustment and Z an anchor variable. As stated earlier, the anchor Z is often a candidate IV, and we will treat it as such in this subsection. Further, let M be a mediator set such that: 1:P is Verma faithful wrt a causal DAG (V∪U) A_1:P is Verma faithful wrt a causal DAG G(V∪ U) 2: satisfies the ordering Z,C<A<M<Y A_2:G satisfies the ordering \Z,C\<A<M<Y 3:Z→A→M→Y exist in A_3:Z→ A→ M→ Y exist in G (11) 4:If Ui∈U causes M it also causes Y A_4:If U_i∈ U causes M it also causes Y Note that 1 A_1 includes ordinary faithfulness as a special case when the intervention set is empty, 2 A_2 imposes a causal ordering of the variables as before, and 3,4 A_3, A_4 are relevance assumptions that ensure the testable implications are not trivial and actually rule out non-identifying edges in G. Proposition 2. Let ℳ1,ℳ2,ℳ3 M_1, M_2, M_3 be backdoor, frontdoor, and IV models respectively and β1,β2,β3 _1, _2, _3 be observed data parameters that are zero iff the independence Y⟂Z∣A,CY \!\!\! Z A,C, the Verma constraint Y⟂Z∣C in p(V)/p(M|A,Z,C)Y \!\!\! Z C in p(V)/p(M|A,Z,C), and the independence M⟂Z∣A,CM \!\!\! Z A,C hold in P respectively. Under assumptions A in (11), β1=0⟹θ=ψ1 in (2) with L=C, _1=0 θ= _1 in eq:backdoor with L=C, β2=0⟹θ=ψ2 in (3) with M=M,L=C∪Z, _2=0 θ= _2 in eq:frontdoor with M=M,L=C∪\Z\, β2=β3=0⟹θ=ψ3 in (4) with Z=Z,L=C. _2= _3=0 θ= _3 in eq:iv with Z=Z,L=C. We will use log(OR(Y,Z∣A,C)) (OR(Y,Z A,C)) and log(OR(M,Z∣A,C)) (OR(M,Z A,C)) as β1 _1 and β3 _3 respectively. For β2 _2, we use a reweighted log-odds log(ORP~(Y,Z∣C)) (OR_ P(Y,Z C)) defined with respect to the distribution P~=p(V)/P(M∣A,Z,C) P=p(V)/P(M A,Z,C). This is estimated using a reweighted regression procedure, similar to those described in robins1997estimation, bhattacharya2022testability, robins2000marginal. The validity of the IV model also relies on assumptions of homogeneity of the treatment effect—we assume such conditions hold apriori and focus only on structural assumptions of the causal models. Further, per Proposition 2, testability of the IV model relies on testability of the frontdoor model, since β2 _2 must be zero as well—other known tests of the IV assumptions are also in over-identified models [kitagawa2015test]. Thus, to test the IV model we use β3~=β2+β3 _3= _2+ _3 under the assumption that these parameters do not exactly cancel out, which can be taken to be a form of faithfulness. We return to Scenario 2 in Section 3. Figure 3(b) is an example of a DAG that satisfies assumptions A in (11). Since influence function-based estimators of parameters encoding Verma constraints, such as log(ORP~(Y,Z∣C))) (OR_ P(Y,Z C))), are underdeveloped, we rely on plug-in estimators for each ψk,βk _k, _k with parametric nuisance estimators. This also serves as a test of the bootstrapping branch of our triangulation procedure. We test two scenarios, one in which frontdoor and IV are correct and testable while backdoor is incorrect, and the other in which only frontdoor is correct and testable (using the blue dashed edge in Figure 3(b)). Figure 4 shows that our estimator maintains robustness in both scenarios. When both frontdoor and IV are correct and testable, our triangulation estimator also exhibits lower variance than relying solely on the IV model. Thus, the estimator has benefits even when an analyst may be certain of some identifying assumptions if these do not yield the most efficient estimator. At n=5000n=5000 and when both frontdoor and IV are correct, we achieve coverage of 97%97\% and 96%96\% for ψ and θ respectively; when only the frontdoor model is correct, the coverage is just below nominal coverage—93%93\% for both ψ and θ. We also use Scenario 2 to test our recommendations in Section 4.3 for setting a. In particular, we use the scenario where the frontdoor and IV models are correct, and try four different settings of the bandwidth parameter: a=n−1/3a=n^-1/3, 1/log(n)1/ (n), 1/log(n)1/ (n), and 1/n1/n. The results are shown in Figure 5. The choices a=n−1/3a=n^-1/3 and a=1/log(n)a=1/ (n) lie within our recommended range and have roughly similar performance. Setting a=1/log(n)a=1/ (n) takes a to zero too slowly, resulting in excess bias and invalid inference for the true parameter θ. Finally, setting a=1/na=1/n takes a to zero faster than our recommendation and results in larger variance. Figure 5: Triangulated point estimates using different settings of a over 200 trials where frontdoor and IV are correct. 5.3 Framingham Data Application Method ACE βk,n,wk,n _k,n,w_k,n Backdoor 0.086(0.048,0.115)0.086\ (0.048,0.115) −0.04,0.45-0.04,0.45 Frontdoor 0.011(0.006,0.016)0.011\ (0.006,0.016) −0.02,0.55-0.02,0.55 IV −3.12(−14.7,11.09)-3.12\ (-14.7,11.09) −0.41,≈0-0.41,≈ 0 Triangulation 0.044(0.026,0.062)0.044(0.026,0.062) – Table 1: Results for the empirical data application. We use data from the Framingham Heart Study [kannel1968framingham] to estimate the effect of blood glucose levels on coronary heart disease. We treat blood glucose levels as a binary treatment variable, where A=1A=1 corresponds to higher than average levels. The outcome Y is a binary indicator of coronary heart disease. Similar to our setup in Section 5.2, we consider backdoor, frontdoor, and IV models, where Z=educational attainmentZ=educational attainment, C=sexC=sex, and M=hypertensionM=hypertension. We intentionally adjust for only one confounder C, so that the assumptions of backdoor in particular are difficult to justify, and so triangulation with other estimates from frontdoor or IV may be important. Results are shown in Table 1, where the last column reports model weights. The frontdoor model receives the highest weight and the IV model the lowest. The backdoor model also retains substantial weight—likely because sex is one of the strongest confounders for both heart disease and blood glucose, making its assumptions approximately valid. The triangulated estimate lies between the backdoor and frontdoor estimates, largely discounting the IV estimate based on diagnostics, and indicates a 4.4% increase in coronary heart disease under intervention to raise blood glucose. As a check, we confirm including additional confounders (e.g., age) does indeed increase the weight applied to the backdoor model, while its point estimate drops to about 0.05. 6 Discussion and Conclusion In this work we proposed a general framework for triangulating causal effects that uses data-driven weights based on measures of model validity. We showed that our triangulation functional trades robustness to causal model misspecification in exchange for some modest bias, as well as additional assumptions like faithfulness and background knowledge. That is, while robustness through triangulation is desirable, our framework is not without limitations. Faithfulness, in particular, is an assumption that is often contested in the causal discovery literature. While violations of faithfulness are rare in a measure-theoretic sense [meek1995strong, boeken2024bayesian], near violations of faithfulness can be fairly frequent [uhler2013geometry]. In our framework, near violations of faithfulness can result in small ε in Theorem 1 and thus large bias. Analysis for cases where ε shrinks with n would provide deeper understanding of the impact of near violations of faithfulness. This is an area of potential future research. We also developed inference strategies with frequentist guarantees for our proposed triangulation functional, and demonstrated their performance through numerical studies and a data application. This included challenging scenarios in which the plurality of models vote similarly on an incorrect causal effect. Extensions to this may focus on sensitivity analysis, deriving other practical scenarios in which the framework can be used, and data-driven selection of the kernel parameter a. Another interesting avenue for future research is to incorporate ideas from Bayesian paradigms for model averaging into our framework, particularly those which provide robustness to unfaithfulness and the need for point identification [silva2016causal]. Acknowledgements.The authors would like to thank Daniel Malinsky for helpful discussions on influence function-based estimation of odds ratios, Hyunseung Kang and Oliver Dukes for helpful discussions on subsampling, and five anonymous reviewers for their helpful peer review. RB would like to thank the Isaac Newton Institute for Mathematical Sciences for the support and hospitality during the programme Causal inference: From theory to practice and back again when some of the work on this paper was undertaken. This work was supported by: EPSRC grant EP/Z000580/1 (RB), NSF CRII grant 2348287 (RB), and NSF DMS 2113171 (TW). References Robust Weighted Triangulation of Causal Effects Under Model Uncertainty (Supplementary Material) In this supplement we provide proofs that were omitted from the main paper for space, as well as details of data generating processes and the specific estimators used in our numerical experiments and data application. Appendix A Proof Of Theorem 1 Proof. Since ∑kwk=1 _kw_k=1, wk≥0w_k≥ 0, and ψk=θ _k=θ for any k∈k∈ C, we have |ψ−θ| |ψ-θ| =|∑k=1Kwkψk−θ| = | _k=1^Kw_k _k-θ | =|∑k=1Kδa(βk)[ψk−θ]|∑k=1Kδa(βk) = | _k=1^K _a( _k) [ _k-θ ] | _k=1^K _a( _k) =|∑k∈ℐδa(βk)[ψk−θ]|∑k∈ℐδa(βk)+∑k∈δa(βk) = | _k∈ I _a( _k) [ _k-θ ] | _k∈ I _a( _k)+ _k∈ C _a( _k) ≤∑k∈ℐδa(βk)maxk|ψk−θ|∑k∈ℐδa(βk)+∑k∈δa(βk) ≤ _k∈ I _a( _k) _k | _k-θ | _k∈ I _a( _k)+ _k∈ C _a( _k) =maxk|ψk−θ|1+Da. = _k | _k-θ |1+D_a. Now by the definition of δa _a, Da D_a =∑k∈δa(βk)∑k∈ℐδa(βk)=∑k∈e−βk2/a2∑k∈ℐe−βk2/a2. = _k∈ C _a( _k) _k∈ I _a( _k)= _k∈ Ce^- _k^2/a^2 _k∈ Ie^- _k^2/a^2. Under our assumptions, βk=0 _k=0 for any model that is correct and testable. If there is at least one correct and testable model, then e−βk2/a2=1e^- _k^2/a^2=1 for this model, and so ∑k∈e−βk2/a2≥1 _k∈ Ce^- _k^2/a^2≥ 1. For the denominator, by monotonicity of x↦e−(x/a)2x e^-(x/a)^2 for x≥0x≥ 0, we get ∑k∈ℐe−βk2/a2≤|ℐ|e−ε2/a2 _k∈ Ie^- _k^2/a^2≤|I|e^- ^2/a^2. The result follows. ∎ Appendix B Deriving The Partial Derivative Vector In Lemma 1 To obtain the variance of the combined estimator using the delta method, we require the following vector of partial derivatives, γ=[∂ψ∂β1,…,∂ψ∂βK,∂ψ∂ψ1,…,∂ψ∂ψK]. γ= [\ ∂ψ∂ _1,…, ∂ψ∂ _K, ∂ψ∂ _1,…, ∂ψ∂ _K\ ]. (12) We divide the task of deriving the partial derivatives ∂ψ∂βk ∂ψ∂ _k and ∂ψ∂ψk ∂ψ∂ _k for all k∈1,…,Kk∈\1,…,K\ into two subsections as follows. B.1 Partial derivative of ψ with respect to βk _k First note that by plugging in the definition of ψ and by linearity of differentiation we have, ∂ψ∂βk=∂βk(∑i=1Kwiψi)=∑i=1K(∂(wiψi)∂βk). ∂ψ∂ _k= ∂ _k ( _i=1^Kw_i _i )= _i=1^K ( ∂(w_i _i)∂ _k ). (13) For each term in (15), we can apply the product rule to get, ∂ψ∂βk=∑i=1K(wi∂ψi∂βk+ψi∂wi∂βk). ∂ψ∂ _k= _i=1^K (w_i ∂ _i∂ _k+ _i ∂ w_i∂ _k ). (14) The estimators ψi _i are not functions of βk _k (even when i=ki=k). Thus, ∂ψi∂βk ∂ _i∂ _k is always 0. So the derivative simplifies to, ∂ψ∂βk=∑i=1Kψi∂wi∂βk. ∂ψ∂ _k= _i=1^K _i ∂ w_i∂ _k. (15) Plugging in the definition for weights wiw_i we get, ∂ψ∂βk=∑i=1Kψi×∂βk(δ(βi)/(λn+∑j=1Kδ(βj))). ∂ψ∂ _k= _i=1^K _i× ∂ _k (δ( _i)/ ( _n+ _j=1^Kδ( _j) ) ). (16) We now separate the above terms in the summation based on whether i=ki=k or not, as these are the only two cases that lead to substantively different partial derivates with respect to βk _k: ∂ψ∂βk=ψk×∂βk(δ(βk)/(λn+∑j=1Kδ(βj)))+∑i≠kψi×∂βk(δ(βi)/(λn+∑j=1Kδ(βj))). ∂ψ∂ _k= _k× ∂ _k (δ( _k)/ ( _n+ _j=1^Kδ( _j) ) )+ _i =k _i× ∂ _k (δ( _i)/ ( _n+ _j=1^Kδ( _j) ) ). (17) We first deal with, ∂βk(δ(βk)/(λn+∑j=1Kδ(βj))). ∂ _k (δ( _k)/ ( _n+ _j=1^Kδ( _j) ) ). By the quotient rule, ∂βk(δ(βk)/(λn+∑j=1Kδ(βj)))=∂δ(βk)∂βk(λn+∑j=1Kδ(βj))−δ(βk)∂∑j=1Kδ(βj)∂βk(λn+∑j=1Kδ(βj))2. ∂ _k (δ( _k)/ ( _n+ _j=1^Kδ( _j) ) )= ∂δ( _k)∂ _k ( _n+ _j=1^Kδ( _j) )-δ( _k) ∂ _j=1^Kδ( _j)∂ _k ( _n+ _j=1^Kδ( _j) )^2. (18) We introduce two helper derivatives of the Dirac delta function to simplify (18): ∂δ(βj)∂βk=∂βk(1|a|πe−(βja)2)=−2βka2δ(βk) when j=k0 when j≠k. ∂\ δ( _j)∂ _k= ∂ _k ( 1|a| πe^- ( _ja )^2 )= cases- 2 _ka^2δ( _k)\ when \ j=k\\ 0 when \ j =k. cases (19) Plugging these into (18) and simplifying gives, ∂βk(δ(βk)/(λn+∑j=1Kδ(βj))) ∂ _k (δ( _k)/ ( _n+ _j=1^Kδ( _j) ) ) =(−2βka2δ(βk))(λn+∑j=1Kδ(βj))−δ(βk)(−2βka2δ(βk))(λn+∑j=1Kδ(βj))2 = (- 2 _ka^2δ( _k) ) ( _n+ _j=1^Kδ( _j) )-δ( _k) (- 2 _ka^2δ( _k) ) ( _n+ _j=1^Kδ( _j) )^2 (20) =(−2βka2δ(βk))(λn+∑j≠kδ(βj))(λn+∑j=1Kδ(βj))2 = (- 2 _ka^2δ( _k) ) ( _n+ _j =kδ( _j) ) ( _n+ _j=1^Kδ( _j) )^2 (21) =−2βkwka2(λn+∑j≠kδ(βj))(λn+∑j=1Kδ(βj)). =- 2 _kw_ka^2 ( _n+ _j =kδ( _j) ) ( _n+ _j=1^Kδ( _j) ). (22) In the above, the first equality comes from plugging in the helper derivatives, the second follows from factorizing a common term, and the third comes from merging terms that correspond to wkw_k. Now we tackle the second set of terms in (17) involving, ∂βk(δ(βi)/(λn+∑j=1Kδ(βj))), ∂ _k (δ( _i)/ ( _n+ _j=1^Kδ( _j) ) ), where i≠ki =k. By the quotient rule again, ∂βk(δ(βi)/(λn+∑j=1Kδ(βj)))=∂δ(βi)∂βk(λn+∑j=1Kδ(βj))−δ(βi)∂∑j=1Kδ(βj)∂βk(λn+∑j=1Kδ(βj))2. ∂ _k (δ( _i)/ ( _n+ _j=1^Kδ( _j) ) )= ∂δ( _i)∂ _k ( _n+ _j=1^Kδ( _j) )-δ( _i) ∂ _j=1^Kδ( _j)∂ _k ( _n+ _j=1^Kδ( _j) )^2. (23) Plugging in the helper derivatives in (19) and simplifying gives, ∂βk(δ(βi)/(λn+∑j=1Kδ(βj))) ∂ _k (δ( _i)/ ( _n+ _j=1^Kδ( _j) ) ) =0−δ(βi)(−2βka2δ(βk))(λn+∑j=1Kδ(βj))2 = 0-δ( _i) (- 2 _ka^2δ( _k) ) ( _n+ _j=1^Kδ( _j) )^2 (24) =2βka2δ(βk)δ(βi)(λn+∑j=1Kδ(βj))2 = 2 _ka^2 δ( _k)δ( _i) ( _n+ _j=1^Kδ( _j) )^2 (25) =2βkwkwia2. = 2 _kw_kw_ia^2. (26) Plugging (22) and (26) back into (17) gives us the following expression for the partial derivative of the combined estimator ψ with respect to βk _k: ∂ψ∂βk=−ψk(2βkwka2(λn+∑j≠kδ(βj))(λn+∑j=1Kδ(βj)))+∑i≠kψi(2βkwkwia2). ∂ψ∂ _k=- _k ( 2 _kw_ka^2 ( _n+ _j =kδ( _j) ) ( _n+ _j=1^Kδ( _j) ) )+ _i =k _i ( 2 _kw_kw_ia^2 ). (27) This can be further simplified to yield the final expression as, ∂ψ∂βk ∂ψ∂ _k =2βkwka2(∑i≠kwiψi−ψk(λn+∑j≠kδ(βj))(λn+∑j=1Kδ(βj))+wkψk−wkψk) = 2 _kw_ka^2 ( _i =kw_i _i- _k ( _n+ _j =kδ( _j) ) ( _n+ _j=1^Kδ( _j) )+w_k _k-w_k _k ) (28) =2βkwka2(ψ−ψk[(λn+∑j≠kδ(βj))(λn+∑j=1Kδ(βj))+wk]) = 2 _kw_ka^2 (ψ- _k [ ( _n+ _j =kδ( _j) ) ( _n+ _j=1^Kδ( _j) )+w_k ] ) (29) =2βkwka2(ψ−ψk). = 2 _kw_ka^2 (ψ- _k ). (30) B.2 Partial derivative of ψ with respect to ψk _k Here we again apply the linearity of differentiation to get, ∂ψ∂ψk=∂ψk(∑i=1Kwiψi)=∑i=1K(∂(wiψi)∂ψk). ∂ψ∂ _k= ∂ _k ( _i=1^Kw_i _i )= _i=1^K ( ∂(w_i _i)∂ _k ). (31) Applying the product rule gives us, ∂ψ∂ψk=∑i=1K(wi∂ψi∂ψk+ψi∂wi∂ψk). ∂ψ∂ _k= _i=1^K (w_i ∂ _i∂ _k+ _i ∂ w_i∂ _k ). (32) Similar to the previous subsection, notice that none of the weights wiw_i are a function of any of the estimators ψk _k (even when i=ki=k). Thus, ∂wi∂ψk ∂ w_i∂ _k is always 0. So the derivative simplifies to, ∂ψ∂ψk=∑i=1Kwi∂ψi∂ψk. ∂ψ∂ _k= _i=1^Kw_i ∂ _i∂ _k. (33) Finally, this simplifies nicely as, ∂ψi∂ψk=1 when i=k0 when i≠k. ∂ _i∂ _k= cases1\ when i=k\\ 0\ when i =k. cases (34) Plugging these in gives us the final expression for the partial derivative of ψ with respect to ψi _i, ∂ψ∂ψk=wk. ∂ψ∂ _k=w_k. (35) Appendix C Proofs of Propositions 1 and 2 Proof of Proposition 1 Proof. entner2013data show under assumptions 1 A_1 and 2 A_2 in (10) that W⊆CW C is a valid backdoor adjustment if Y⟂̸⟂Z∣WY \!\!\! Z W and Y⟂Z∣A,WY \!\!\! Z A,W. The additional assumption 3 A_3 we make in (10) ensures that Y⟂̸⟂Z∣WY \!\!\! Z W is already true due to the existence of the path Z→A→YZ→ A→ Y. Thus, under assumptions 1,2,3 A_1, A_2, A_3, the independence Y⟂Z∣A,WY \!\!\! Z A,W alone is sufficient to ensure that model ℳk M_k is correct (i.e., W is a valid backdoor adjustment set). The conclusion then follows, as we suppose that βk _k is an observed data parameter such that βk=0⇔Y⟂Z∣A,W _k=0 Y \!\!\! Z A,W. ∎ Proof of Proposition 2 Proof. The argument for β1=0 _1=0 implying correctness of the backdoor model with adjustment set C is essentially the same as the proof of Proposition 1 with W=CW=C, since the assumptions used in Proposition 1 are a superset of those in (10). For the second implication, we make an argument similar to the one used in bhattacharya2022testability. Under the Verma faithfulness assumption 1 A_1, we know that if β2=0 _2=0, the causal DAG (V∪U)G(V∪ U) must support identification of P(A,Z,C,Y∣do(m))P(A,Z,C,Y (m)) by the g-formula and Y⟂d-sepZ∣CY \!\!\! _d-sepZ C in do(m)G_do(m). By assumption 3 A_3 we know that Z→AZ→ A exists in G. We now show that the existence of A→YA→ Y in G contradicts existence of the Verma constraint under faithfulness. Suppose A→YA→ Y does exist in G. Then, Y⟂̸⟂d-sepZ∣CY \!\!\! _d-sepZ C in do(m)G_do(m) due to the open path Z→A→YZ→ A→ Y, which is a contradiction. Thus, the frontdoor exclusion restriction of no A→YA→ Y is satisfied. Per tian2002general, the distribution P(V∖A∣do(a))P(V \A\ (a)) where A is a single treatment variable is identified if and only if there is no path from A to any child X of A of the form A←⋯→XA←·s→ X such that every collider on the path is an observed variable Vi∈V_i∈ V and every non-collider on the path is an unmeasured variable in Ui∈U_i∈ U. By the causal ordering assumption 2 A_2 and the previous argument ruling out the existence of A→YA→ Y, the only child of A is M, so we only need to show that no such path exists from A to M. The existence of such a path between A and M also contradicts the presence of the Verma constraint, as it would imply the existence of some Ui∈U_i∈ U that causes M and thus also causes Y (by assumption 4 A_4), preventing identification of p(A,Z,C,Y∣do(m))p(A,Z,C,Y (m)). Thus, the second frontdoor restriction is satisfied and P(Z,C,M,Y∣do(a))P(Z,C,M,Y (a)) is identified as [tian2002general] (∑a′P(a′,Z,C)⋅[Y∣a′,M,C,Z])⋅p(M∣A=a,Z,C)⋅ ( _a P(a ,Z,C)·E[Y a ,M,C,Z] )· p(M A=a,Z,C)· It is straightforward then to obtain the frontdoor functional in (3) for θ by summing over Z,CZ,C to obtain P(Y∣do(a))P(Y (a)). Finally for the third implication, Z already satisfies the IV relevance condition by assumption 3 A_3. The second condition for IV validity is that Z should be d-separated from Y given C in a graph where we delete the outgoing edges from A. This is satisfied when M⟂Z∣A,CM \!\!\! Z A,C and the Verma constraint holds per Corollary 1.1 in bhattacharya2022testability. ∎ Appendix D Details Of Data Generating Processes In this subsection we describe the data generating processes (DGPs) underlying each of our numerical experiments. In our descriptions, we define expit(x)≔1/(1+exp(−x))expit(x) 1/(1+exp(-x)) for x∈ℝx . We use Bern(p)Bern(p) as shorthand for the Bernoulli distribution with probability p and N(μ,σ2)N(μ,σ^2) as shorthand for the normal distribution with mean μ and variance σ2σ^2. D.1 DGPs for Section 5.1 When the U1→ZU_1→ Z edge is absent, the data are generated according to Figure 3(a) as, Z Z ∼Bern(0.5) (0.5) U1,U2,C1,C2,C3 U_1,U_2,C_1,C_2,C_3 ∼N(0,1) N(0,1) C4 C_4 ∼N(−2.5⋅U1+2⋅U2, 1) N (-2.5· U_1+2· U_2,\ 1 ) C5 C_5 ∼N(−2.5⋅U1+2⋅U2, 1) N (-2.5· U_1+2· U_2,\ 1 ) A A ∼Bern(expit(2.75⋅Z−3⋅U1+C1+C2+C3)) (expit(2.75· Z-3· U_1+C_1+C_2+C_3) ) Y Y ∼Bern(expit(1.5⋅A+2⋅U2+C3+C4+C5)) (expit(1.5· A+2· U2+C_3+C_4+C_5) ) When the U1→ZU_1→ Z edge is present, the data are generated according to Figure 3(a) as, U1,U2,C1,C2,C3 U_1,U_2,C_1,C_2,C_3 ∼N(0,1) N(0,1) Z Z ∼Bern(expit(U1)) (expit(U_1) ) C4 C_4 ∼N(−2.5⋅U1+2⋅U2, 1) N (-2.5· U_1+2· U_2,\ 1 ) C5 C_5 ∼N(−2.5⋅U1+2⋅U2, 1) N (-2.5· U_1+2· U_2,\ 1 ) A A ∼Bern(expit(2⋅Z−3⋅U1+C1+C2+C3)) (expit(2· Z-3· U_1+C_1+C_2+C_3) ) Y Y ∼Bern(expit(A+2⋅U2+C3+C4+C5)) (expit(A+2· U2+C_3+C_4+C_5) ) D.2 DGPs for Section 5.2 When both frontdoor and IV models are correct, data are generated according to Figure 3(b) without the blue dashed edge as, Z Z ∼Bern(0.5) (0.5) C C ∼N(0,1) N(0,1) U U ∼N(0,1) N(0,1) A A ∼Bern(expit(2⋅Z+2⋅C+2⋅U)) (expit(2· Z+2· C+2· U) ) M M ∼Bern(expit(−2+4⋅A−0.5⋅C)) (expit(-2+4· A-0.5· C) ) Y Y ∼N(2⋅M+2⋅C+2⋅U, 1) N (2· M+2· C+2· U,\ 1 ) When only the frontdoor model is correct, data are generated according to Figure 3(b) with the blue dashed edge as, Z Z ∼Bern(0.5) (0.5) C C ∼N(0,1) N(0,1) U U ∼N(0,1) N(0,1) A A ∼Bern(expit(Z+C−0.5⋅U)) (expit(Z+C-0.5· U) ) M M ∼Bern(expit(−1+2⋅A−Z+C)) (expit(-1+2· A-Z+C) ) Y Y ∼N(2⋅M−0.75⋅C−2⋅U, 1) N (2· M-0.75· C-2· U,\ 1 ) Appendix E Estimators Used In Numerical Experiments And Data Application Below we describe the specific estimators used in our numerical experiments and data application in more detail. E.1 Estimators Used In Section 5.1 Estimators for βk _k For each log-odds ratio βk _k, we use an influence function-based estimator of log(OR(Y,Z|A,W)) (OR(Y,Z|A,W)) from tchetgen2010doubly and tan2019doubly. An R implementation of their method is publicly available from Wu and Malinsky at https://github.com/chaoqiw0324/ortest; we translate this to Python for our purposes. First, let ζ(a,w)≔[Y∣z0,a,w]ζ(a,w) [Y z_0,a,w] and η(a,w)≔[Z∣y0,a,w]η(a,w) [Z y_0,a,w]. Let O=Y,A,Z,WO=\Y,A,Z,W\. Then, an unbiased estimating function of βk _k is g(o,ζ,η,βk)=(y−ζ(a,w))⋅(z−η(a,w))⋅e−βk⋅y⋅z. g(o,ζ,η, _k)=(y-ζ(a,w))·(z-η(a,w))· e^- _k· y· z. Since g is an estimating function, we can construct a point estimate βk,n _k,n of the log-odds ratio by first constructing estimators ζn _n and ηn _n of ζ and η respectively, and using these to obtain the value of βk _k that solves the equation ∑i=1ng(oi,ζn,ηn,βk)=0 _i=1^ng(o_i, _n, _n, _k)=0. This estimator is asymptotically linear under doubly robust conditions on the nuisance estimators ζn _n and ηn _n. Further, by standard Z-estimator theory, the influence function ϕβk=B−1g(o,ζ,η,βk) _ _k=B^-1g(o,ζ,η, _k), where B−1=[−∂g/∂βk]B^-1=E[-∂ g/∂ _k]; see for example the review in cole2025five. We can then construct an estimator of the influence function ϕβk,n _ _k,n as 1n∑i=1n(−∂g(oi,ζn,ηn,βk)∂βk)−1⋅g(oi,ζn,ηn,βk)|βk=βk,n. 1n _i=1^n ( -∂ g(o_i, _n, _n, _k)∂ _k )^-1· g(o_i, _n, _n, _k)\ |_ _k= _k,n. Estimators for ψk _k For each ℳk M_k that uses an adjustment set W, we use the AIPW estimator [bang2005doubly] of the backdoor formula (2) with L=WL=W. Let π(w)≔P(A=1∣W=w)π(w) P(A=1 W=w) and μ(a,w)≔[Y∣A=a,W=w]μ(a,w) [Y A=a,W=w]. Then the AIPW estimator is an asymptotically linear estimator of ψk _k under doubly robust conditions on the nuisance estimators πn _n and μn _n of π and μ respectively, with influence function ϕψk _ _k given by ϕψk=y−μ(a,w)a−π(w)π(w)(1−π(w))+μ(1,w)−μ(0,w)−ψk. _ _k=\y-μ(a,w)\ \ a-π(w)π(w)(1-π(w)) \+\μ(1,w)-μ(0,w)\- _k. E.2 Estimators Used In Section 5.2 and Data Application For these experiments we construct plug-in estimators of βk,n _k,n and ψk,n _k,n based on parametric nuisance estimation. Estimators for βk _k For β1 _1 and β3 _3, we construct parametric nuisance estimators ζn _n and ηn _n of ζ(y,c,a,z)≔P(Y∣A,C,Z)ζ(y,c,a,z) P(Y A,C,Z) and η(m,a,z,c)≔P(M∣A,Z,C)η(m,a,z,c) P(M A,Z,C) respectively. In the parametric models we consider, estimates of the log-odds ratio log(OR(Y,Z∣A,C)) (OR(Y,Z A,C)) and log(OR(M,Z∣A,C)) (OR(M,Z A,C)) are given by the coefficients of Z in ζn _n and ηn _n respectively. For β2 _2, the parameter used to test the Verma constraint, let α≔P(Y∣Z,C)α P(Y Z,C) and g(y,z,c,α)g(y,z,c,α) be any unbiased estimating function used to construct a parametric nuisance estimator αn _n of α. The reweighted log-odds ratio log(ORP~(Y,Z∣C)) (OR_ P(Y,Z C)) is obtained as the coefficient of Z, one of the parameters estimated in αn _n, where αn _n is the solution to a set of reweighted estimating equations ∑i=1ng(yi,zi,ci,α)/ηn(mi,ai,zi,ci)=0 _i=1^ng(y_i,z_i,c_i,α)/ _n(m_i,a_i,z_i,c_i)=0 and ηn _n is a parametric nuisance estimator of η as before. In terms of practical implementation, log(ORP~(Y,Z∣C)) (OR_ P(Y,Z C)) can be treated as the coefficient of Z in a reweighted (linear or logistic) regression of the outcome on the covariates C and anchor variable Z, where the weights are given by 1/η(m,a,z,c)1/η(m,a,z,c). Estimators for ψk _k For ψ1 _1, the backdoor functional (2) with L=CL=C, we construct a plug-in estimator under parametric specification of the outcome regression model μ≔[Y∣A,C]μ [Y A,C]. For ψ2 _2, the frontdoor functional (3) with L=C∪ZL=C∪\Z\, we use a plug-in estimator of the more convenient dual IPW functional instead under parametric specification of the mediator model η≔P(M∣A,Z,C)η P(M A,Z,C) [fulcher2020robust, bhattacharya2022semiparametric] ψ2,Dual IPW=[P(M∣A=1,Z,C)P(M∣A,Z,C)×Y]−[P(M∣A=0,Z,C)P(M∣A,Z,C)×Y]. _2,Dual IPW=E [ P(M A=1,Z,C)P(M A,Z,C)× Y ]-E [ P(M A=0,Z,C)P(M A,Z,C)× Y ]. Finally, for ψ3 _3, the IV functional (4) with L=CL=C, we construct a plug-in estimator under parametric specification of the nuisance function in the numerator ν≔[Y∣Z,C]ν [Y Z,C] and denominator ξ≔[A∣Z,C]ξ [A Z,C].