A Free Lunch with Influence Functions? An Empirical Evaluation of Influence Functions for Average Treatment Effect Estimation
Authors: Matthew James Vowels
Sina Akbari
Necati Cihan Camgoz
Richard Bowden
Venue: Accepted by TMLR ()

Abstract: The applications of causal inference may be life-critical, including the evaluation of vaccinations, medicine, and social policy. However, when undertaking estimation for causal inference, practitioners rarely have access to what might be called `ground-truth' in a supervised learning setting, meaning the chosen estimation methods cannot be evaluated and must be assumed to be reliable. It is therefore crucial that we have a good understanding of the performance consistency of typical methods available to practitioners. In this work we provide a comprehensive evaluation of recent semiparametric methods (including neural network approaches) for average treatment effect estimation. Such methods have been proposed as a means to derive unbiased causal effect estimates and statistically valid confidence intervals, even when using otherwise non-parametric, data-adaptive machine learning techniques. We also propose a new estimator `MultiNet', and a variation on the semiparametric update step `MultiStep', which we evaluate alongside existing approaches. The performance of both semiparametric and `regular' methods are found to be dataset dependent, indicating an interaction between the methods used, the sample size, and nature of the data generating process. Our experiments highlight the need for practitioners to check the consistency of their findings, potentially by undertaking multiple analyses with different combinations of estimators. 

Published in Transactions on Machine Learning Research (03/2023)
A Free Lunch with Influence Functions? An Empirical Eval- uation of Influence Functions for Average Treatment Effect Estimation
Matthew J. Vowels
natthew.wowellsee www.wlelsswwie.ocmm
Institute of Psychology University of Lausanne Switzerland
Sina Akbari BAN
wha.akbarr@eepf..cm
EPFL
Switzerland
Necati Cihan Camgoz
w.oomgoz:@surrey.cc.uu
CVSSP
University of Surrey U.K.
Richard Bowden CVSSP
r.bonden@ssurrey.comm
University of Surrey U.K. Reviewed on OpenReview: https: // openreview. net/ forum? id= dQuBRgCiLr
Abstract
The applications of causal inference may be life-critical, including the evaluation of vac- cinations, medicine, and social policy. However, when undertaking cstimation for causal inference, practitioners rarely have access to what might be called 'ground-truth' in su- pervised learning setting, meaning the chosen estimation met cannot be evaluated and must be assumed to be reliable. is therefore crucial that we have a good understanding of the performance consistency of typical methods available to practitioners. we provide comprehensive evaluation of recent semiparametric methods (including neural network approaches) for average treatment effect estimation. Such methods have been pro- posed means to derive unbiased causal effect estimates and statistically valid confidence intervals, even when using otherwise -oramarme,, data-adaptive machine learning tech- niques. We also propose a new estimator 'MultiNet', and a variation on the semiparametric update step 'MultiStop', which we evaluate alongside existing approaches. The performance of both semiparametric and 'regular' methods are found an interaction between the methods used, the sample size, and nature of the data generating process. Our experiments highlight the need for practitioners to check the consistency of their findings, potentially by undertaking multiple analyses with different combinations of estimators.
this work
dataset dependent, indicating
1 Introduction
Causal inference and causal effect estimation are of the upmost importance for policy making (Kreif & DiazOrdaz 2019). the development of medical treatments (Petersen et al. 2017), the evaluation of evidence

---

Published in Transactions on Machine Learning Research (03/2023)
within legal frameworks (Pearl 2009 Siegerink et al., 2016), and others. Unfortunately, causal inference is different to the typical supervised learning paradigm in machine learning insofar as ground-truth for the causal effect is not available. Practitioners are therefore forced to take a leap of faith, hoping that their estimation methods are reliable and yield meaningful estimates. It is therefore crucial that we have a good understanding of the consistency of the methods we use, before deploying them in practice.
In restricted circumstances we can use well-studied techniques which have known performance guarantees (e.g., linear models for linear data generating processes and Gaussian errors). Indeed, most methods being utilized in empirical fields such psychology epidemiology parametric models (van der Laan & Rose 2011 Blanca ct al., 2018). which are convenient because they facilitate closed-form statistical inference and confidence intervals (e.g. for the purpose of null hypothesis testing). However, there now exist numerous powerful, non- paraammme,r,,, data-adaptive machine learning techniques which do not require strong assump- tions about the functional form (e.g., such as the assumption of linearity). These methods can be "augmented" with semiparametric techniques, which can be usod to reduce bias and derive confidence intervals for per- forming statistical inference. Practitioners may also find themselves choosing between the alluring 'deep learning' based methods and those more conventional machine learning approaches which perhaps, rightly or wrongly, have less associated hype.
Unfortunately, the empirical consistency of the performance of these machine learning approaches for average treatment effect estimation (both with and without the use of semiparametric techniques) is not clear. Indeed, recent work by Curth et al. (2021b) has already highlighted that conventional benchmark datasets for evaluating estimators can be misleading. Furthermore, a theoretical framework for establishing the statistical guarantees of neural-network based approaches is yet elusive (Curth et al.) 2021a), although notable recent contribution is presented by Farrell et al. 20221. In this work we undertake much needed evaluation of machine learning and semiparametric approaches to average treatment effect estimation.
The associated contributions of this paper are:
A new update step method 'MultiStep' (see Section which attempts to improve upon existing update methods by continuously optimizing the solution according to two criteria which characterize the optimum solution (namely, finding the IF with the smallest expectation and variance)
A new method 'MultiNet' (see Section with a single NN
5 which attempts to mimic the performance of an ensemble
An extensive comparison (see Section 6 of the estimation performance of NNs and othe algorithms with and without semiparametrie techniques
We evaluate causal inference task performance in terms of (a) precision in estimation (and the degree to which we can achieve debiasing using semiparametric techniques). (b) double robustness (the extent to which we can improve the robustness of our estimators using semiparametric techniques), and (c) normality of the distribution of estimates (thus, by implication, whether is possible to use closod-form expressions for confidence intervals and statistical inference). In order to evaluate the methods in terms of these three criteria, we 10 different datasets (varying in positive violations, non non-linearity, and sample size). We evaluate across three relevant metrics, and implement all possible combinations of approaches to the estimation. Our findings suggest that evaluation on handful of benchmark datasets is not sufficient for helping us understand the likelihood of causal estimation methods to perform well in real-world applications.
In general, the benefits of influence functions and the performance of a range of contemporary estimators is found to be inconsistent Whilst our MultiNet and MultiStep methods provide competitive performance across datasets, and whilst we observe that initial estimation methods can sometimes benefit. from the application of the semiparametric techniques, our evaluation highlights interactions between the underlying data generating process, sample sizes, and the estimators and update steps used. The conclusion is thus that
II is hopefully clear from the context, but consistency have meanings. Firstly in statistical techhiccal terms, for whether an estimator converges towards value as the sample size increases. Sccondly, in the practical/ cemirical sense, use it to describe whether not particular method, algorithm, estimation approach performs well across range of scenarios, over range of metrics.
2

---

Published in Transactions on Machine Learning Research (03/2023)
practitioners should take care when interpreting their results, and attempt to validate them by undertaking multiple analyses with different estimators. This is particularly important for the task of causal inference where, in real-world applications, ground truth is unlikely to be available.
The paper is structured as follows: We begin by reviewing previous work in Sec. 2 and provide background theory on the motivating case of estimating causal effects from observational data in Sec. 3 Tn this section, we also provide top level introduction to II's (Sec. 3.2 p. We present our own update approach MultiStep
in Sec. [ Our NN method MultiNet is presented in Sec. 5 The evaluation is presented in Sec. and
at the beginming of this section, we summarise the open questions which inform subsequent evaluation design. We present and discuss results in Sec. 6.3 and Sec.
6.4 and we provide a summary of the experiments, Finally, a broader impact statement is included in
conclusions, and opportunities for further work in Sec.
Sec. g
Code for reproducing experiments is provided
here:
httpe://gthhu..cc////tttterrere
FreelunchSemiPPrrmetii
2 Context and Previous Work
Being able to perform statistical tests and reliably quantify uncertainty is especially important when eval- uating the elficacy of treatments or interventions. One approach to perform such tests is by assuming a parametric model for the underlying generating mechanism, and e.g. normally distributed estimates. How- ever, it has been argued that lincar models are incapable of modeling most realistic data generating processes and that we should instead be using modern machine learning techniques (van der Laan & Rose, 2011 van
der Laan & Gruber 2012: van der Laan & Starmans 2014 Vowels, 2021). Unfortunately, most machine learning models are nou-parametric insofar as the estimates derived using such techniques are not directly parameterizable as (e.g.) a Gaussian with a mean and variance. As such, the estimates derived using such non-parametric techniques are not readily amenable to null-hypothesis significance testing other common statistical inference tasks. Furthermore, even though machine learning algorithms are more flexible, they are still likely to be biased because they are not targeted to the specific parameter of interest (van der Laan & Rose 2011). So, what can we do?
By leveraging concopts from the field of semiparametric statistics we can try to address these issucs. Indeod, by combining elements of semiparametrie theory with machine learning methods, we can, at least theoret- ically, enjoy the best of both worlds: We can avoid having to make unreasonably restrictive assumptions about the underlying generative procoss, and can nonetheless undertake valid statistical inference. Further- more, we can also leverage estimator update process to achieve greater precision in existing estimators, without needing additional data (van der Laan & Rose 2011 Tsiatis 2006 Bickel et al.] 1998), an advantage which we might call a 'free lunch? and one we wish to explore and test empirically in this paper.
In terms of the available estimation techniques, one approach which combines machine learning and semi- parametric theory is known as targeted learning (van der Laan & Rose 2001 van der Laan & Starmans 2014 3 This technique, and many related techniques involving influence functions (IFs) and semiparametric theory, have primarily been popularized outside the field of machine learning. In parallel, machine learning
has focused on the development of equivalent methods using deep neural network (NN) methods for causal
inference (sce e.g., BBca et al., 2020 Wu & Fukumizu, 2020: Shalit et al. 2017 Yoon et al.) 2018 Louizos eet al.] 2017), which, owing to their 'untargeted' design (more on this below), may exhibit residual bias. As such, many of the principles and theory associated with semiparametrics and IF's are underused and underappreciated within the machine learning community, and remains unknown to what extent these techniques can be appliod NN based estimators.
The possible applications of semiparametrics in machine learning are broad but under-explored, and IFs in particular have only seen sporadic application in explainable machine learning (Koh & Liing 2017 Sani
2The term 'free lunch" reference to the adage of unknown origin (but probably North American) 'there ain't such thing free lunch'. It was famously used by Wolpert and Macready in the context of optimization (Wolpert & Macready 1997).
For an overview some other related methods see Curth al. 2021a).
3

---

Published in Transactions on Machine Learning Research (03/2023)
et al.
2020), natural language processing (Han et al. 2020) models, causal model selection (Alaa & van der Schaar 2019) and uncertainty quantification for deep learning (Alaa & van der Schaar 2020). Outside of
machine learning, in particular in the fields of epidemiology and econometrics, semiparametric methods are becoming more popular, and include targeted learning (van der Laan & Rose 2011 and the well-known double machine learning approach by Chernozhukov el al. (2018). In stalistics, alternatives have been developed which include doubly robust conditional ATE estimation (Kennody, 2020 and IF-learning (Curth et al. 2021a).
Within the field representing the confluence of causal inference and machine learning, the focus seems to have been on the development of NN methods. See, for example, CEVAE (Louizos et al.) 2017p, CFR-Net (Shalit et al. 2017), GANITE YYoon et al. 2018), Intact-VAB (Wu & Fukumizu 2022) etc. However, these methods have been developed without a consideration for statistical inference or semiparametric theory, and this gap has been noted by Curth et al. (2021b) and Curth & van der Schaar (2021). Indeod, to the best of our knowledge, the application of semiparametric theory to dcbias neural network-based estimators has only be used three times in the field representing the confluence of machine learning and causal inference. Firstly, in DragonNet (Shi et al., 2019), method designed for ATE estimation; secondly in TVAE (Vowels) 2021). a variational, latent variable method for conditional ATE and ATE estimation; and thirdly, by Farrell et al. (2021) where a restricted class of multilayer perceptrons were evaluated for their performance potential as plug-in estimators for semiparameteric estimation of causal effects and shown to yield promising performance. The first two methods incorporate targeted regularization, but do not readily yield statistical inference because to do so requires asymptotic normality (and this is not evaluated in the studies) as well as explicit evaluation of the IF. More broadly, semiparametrics has been discussed in relation to theory in machine learning, for example [Bhattacharya et al. (2020) provides discussion of influence functions in relation to Directed Acyelic Graphs with hidden variables, Rotnitzky & Smucler (2020) and Henckel et al. (020) discuss the application of semiparametric techniques for identifying efficient adjustment sets for causal inference tasks, Zhong & Wang (2021) apply semi-parametrics with deep neural networks to achiove statistical inference in the partially linear quantile regression setting, and Juun et al. (2020) generalize the coverage of work on semiparametric estiraatio to general causal estimands. However, in general the work is quite sparse, particularly in relation to the applicability of the theory to neural networks, and the accessibility of the relevant theory to general practitioners of machine learning.
ct al.
Prior comparisons of the performance of semiparametric approaches and causal inference methods exist. For example, the robustness of targeted learning approaches to causal inference on nutrition trial data was presented by Li et al. (2022) and includes a useful summary table of previous findings and includes its own evaluations. However, il does not include comparisons with NN-based learners, and seeks the answers to different questions relevant to practitioners in the empirical fields. Another example evaluation was undertaken by [LquueFernaaddz et al. (2018) but has a didactic focus. Finally, and more broadly, Cuuth et al. (2221b) evaluate the use of benchmark datasets for the evaluation of treatment effect estimators. They find that certain datasets "..systematically favor some algorithms over others', highlighting need for fairer evaluations of techniques designed for causal inference. This kind of benchmark dependency, referred to as the 'enncmaakk lottery', has been noted before in other supervisede ieeee-eeeeeelinngggggggg domains such as computer vision and natural language processing (Dehghani et al.] arxii preprint). It is worth noting, however, that in contrast to supervisod learning tasks, ground-truth rarely available for causal inference, making it all the more crucial that we undersland the performance of our estimators.
3 Background
Regarding notation, we use upper-case letters e.g. A, B to denote random variables, and bold font, upper- case letters to denole sets of random variables e.g. A, B. Lower-case a and b indicate specific realisations of random variables A and B. Specifically, we use x; ~ P(X) E R" to represent the m-dimensional, pre- treatment covariates (we use bold symbols to signify multi-dimensional variables) for individual assigned factual treatment (T) (0, 1} resulting in outcome y/ ~ P(Y T). Together, these constitute dataset D (Dv, 41 x:] where n is the sample size, sampled from a 'true' population distribution P.
4

---

Published in Transactions on Machine Learning Research (03/2023)
Covariates
do
T
tj
X
X
Treatment
Outcome
T
Y
T
(a)
(b)
Figure 1: Directed Acyelic Graphs (DAGs) for estimating the effect of treatment T t on outcome Y with confounding X.
3.1 Causal Inference
Fig. Da is characteristic of observational data, where the outcome related to the covariates as well as the treatment, and treatment is also related the covariates. For example, if consider age to be typical covariate, young people may opt for surgery, whereas older people may opt for medication. Assuming that age-related risk mechanism exists, then age will confound our estimation of the causal effect of mreatment outcome.
One of the most common causal estimands is the Average Treatment Effect (ATE):
T(x) = Ex P(X) Pyy- 1)x 1)x x/(8l P(
(1)
Here, the use of the do operator [Pear 2009) in do(T = 1) and do(T = 0) simulates interventions, setting treatment to a particular value regardless of what was observed. In Fig. D such an intervention removes the dependence T' on X, and this graph is the same as the one for an RCT, where the treatment unrelated the covariates. Using do-calculus we can establish whether, under number of strong assumptions the desired causal estimand whether the effect is identifiable. For the graph in Fig. [h, the outcome under intervention can be expressed Eyy [do(r=<'))M) ( [v] which is estimable from observational data. Here, t' is the specific intervention of interest (e.g., t' = 1). In particular, it tells us that adjusting for the covariates X is sufficient to remove the bias induced through the "backdoor' path X - - This particular approach is sometimes referred to as backdoor adjustment.
be expressed in terms of function of the observed distribution, and thus
One may use regression to approximate the outcome under intervention, and indeed, plug-in estimators Q be used for estimating the difference between the outcome under two different interventions treatment (the Average Treatment Effect (AE)
+(2:x) - - -, x;)),
(2)
1=1
We use the circumflex/hat (1) notation to designate estimated (rather than true/population) quantity, In the simplest case, may use linear or logistic regression for the estimator Q, depending on whether the outcome is continuous or binary. Unfortunately, if one imagines the true joint distribution to fall somewhere within an infinite set of possible distributions, we deliberately handicap ourselves by using a family of linear models because such family is unlikely to contain the truth. The consequences of such model misspecification can be severe, and results in biased estimates (Vowels) 2021 van der Laan & Rose 2011). In other words, no matter how much data we collect, our estimate will converge to the incorrect value, and this results in false positive rate which converges to 100%. This clearly affects the interpretability and reliability of null-hypthhesi tests. Furthermore, even with correct specification of our plug-in in estimators (i.e., we have achieved identification of the causal effect and our assumptions about the structure and functional relationships hold), our models are unlikely to be 'targeted' to the desired estimand, bocause they often
4These assumptions the Stable Unit Treatment Value Assumption (SUTVA), Positivity, and Ignorabil- ity/ /Unnoonfonndeeeesssss see Section A below more information.
5

---

Published in Transactions on Machine Learning Research (03/2023)
estimate quantities superftuous to the estimand but necessary for the plug-in estimator (e.g., other relevant factors or statistics of the joint distribution). As result, in many cases there exist opportunities to reduce residual bias using what are known as influence functions.
3.2 Influence Functions
Semiparametric theory and, in particular, the concept of Influence Functions (IFs), are known be
challenging to assimilate (Fisher & Kennedy, 2019: Levy 2019 Hines et al.] 2021). Here attempt to
provide a brief, top-level intuition, but detailed exposition lies beyond the scope of this paper. Interested
readers are encouraged to consider work by Kennedy (2016): | Fisher & Kennedy (2019); Hampel (1974):
Ichimura & Newey [2021]; Hines et al. [2021]: 0; Bickel et al. 1998): ; Newey
[1994] 1990 Chernozhukor et al.
2018 van der Laan & Rubin (2006), and Tsiatis 2006
An estimator (Pr)) for an estimand (P) (for example, the ATE) has an IF, 0, if it can be expressed follows:
1 - Σ p(v. Pa) + o,((/Vm) (PP) ((P)=v(PP) ) +
(3)
where ys is a sample from the true distribution P, Pa is the empirical distribution or, alternatively, a model of some part thereof (e.g., a predictive distribution parameterized by a NN, or histogram estimate for density function, etc.), d function (the Influence Function) with mean of zero and finite variance Tsiatis 2006 pp.21), and Op (1) ) is error term that converges in probability to zero.
In many cases Pn is not equivalent to the sample distribution, perhaps because some or all of is being modelled with estirattors. As result, the error does not converge in probability zero and some residual error remains. The IF ó is being used to model the residual bias that stems from the fact that Pr is no longer equivalent to a direct sample from P. We can imagine the sample distribution P,, lies a linear path towards the true distribution P. This linear model can be expressed using what is known as a parametric submodel, which represents a family of distributions indexed by a parameter 6. The direction associated with Pe can then be expressed as a pathwise derivative in terms the function representing our estimand V
The IF itself is function which models how much estimate deviates from the true estimand, up to the error term. If an estimator can be written in terms of its IF, then by central limit theorem and Slutsky's theorem, the estimator converges in distribution to normal distribution with mean zero and variance equal to the variance of the IF. This is a key result that enables us to derive confidence intervals and perform statistical inference (such as mull hypothesis significance testing). Note that the convergence in distribution to a normal distribution does not impact our assumptions about the functional form governing the data generating process itself. In other words, by deriving a normally distributed parameter estimate, we do not sacrifice flexibility in how we model the observed distribution, or flexibility in the functions used to derive the estimates themselves. Appendix B contains information about how to leverage this normality to derive confidence intervals and p-values, and Appendix C derives the IF for general estimands as well as presenting an algorithm for doing so automatically.
4 The MultiStep Update Method
This section presents the first contribution a new process for updating initial estimators using the influence function. The current work concerns the estimation of the Average Treatment Effect (ATE). For the DAG: T Y. T +Y (also see Fig. [b), the IF for the ATE can be expressed (Hines 2021 van der
Laan & Rose 2011

---

Published in Transactions on Machine Learning Research (03/2023)
(1(1) 1- 1-6;(0)) (( - G(x)
(i))(()))) (1, (ATE
PATE
(4)
where = (X), In a slight abuse of notation, dg is the Dirac delta function at the point at which y = D, where y can be a datapoint in our empirical sample (note the shift from specific datapoint y generic empirical samples f). In order to evaluate this we need to evaluate it at Pn, and we also need plug-in estimators ((x) (((x) (propensity score model), and (1, Re P [y] (outcome model).
4.1 Updating Initial Estimators
If we can estimate the IF o then we can update our initial estimator (Pr) according to Eq. in order to reduce the residual bias which the IF is essentially modeling. To be clear, this mcans we can improve our initial NN estimators, without needing additional data. We consider four ways to leverage the IF reduce bias which we refer to as (1) the one-step update, (2) the submodel update (sometimes referred to as targeted update), (3) our own proposed MultiStep procedure (which we present in this section), and (4) targeted regularization. A brief description of methods
(1), (2), and (4) is provided in Appendix D) The first. three approaches can be trivially applied to estimators which have already been trained, making them attractive as post-processing methods for improving estimation across different application arcas. To illustrate these approaches, we consider the ATE be our chosen target estimand, the IF for which is defined in Equation Fig. 2 provides an illustration of the components involved in updating the estimator in the context of ATE estimation.
In order to motivate our own proposed MultiStep method, we begin by noting the limitations of the one-step and submodel update processes. In general, these updates are performed only once (Hines et al. 2021 van der Laan & Rose 2011). In other words, the bias of our initial estimator must be able to be approximated by a linear submodel, such that taking a single step in the direction of the gradient sufficient. We attempt to improve the empirical robustness of the one-slep and submodel update steps by modifying the objective in the update step itself.
4.2 Explicit Optimization Goals
Under the assumptions described above, the one-step and the submodel update approaches yield what is known as the efficient inffuence function. Importantly, the efficient influence function is the one which meets two conditions: (A) the influence function with the smallest variance, and (B) the influence function with a mean of zero, i.e., Σ. p(z;, P.) (Triatis 2006). Our idea is thus explicitly specify these two conditions as goals in continuous optimization problem.
In contrast, the one-step and submodel processes achieve these two conditions indirectly by following the Von Mises gradient step and finding the least-squarcs (or maximum-licclhhood) solution to a supervised learning problem, respectively (see Eq. 15 in the Appendix for more details). In particular, the one-step approach updates an initial estimator Qtt, with the empirical estimate of the influence function (which serves as the gradient), whilst the submodel method updates it by adjusting for a some amount & of a biasing quantity known the clever covariate (z,). We refer to both update processes 'indirect' because, in contrast, the objective used to find the influence function with the smallest. variance and zero mean can be specified explicitly.
We refer to our update variant as MultiStep because whilst it still uses the linear submodel (see Eq. 15 and Appendix D for more details), we optimize the expression 5 below by searching over & E such that the mean and variance of the influence function are minimized:
minger [o [[ P)]1 + az Varlo(z, P)]]
(5)
In words, rather than implicitly finding the solution to the IF indirectly via gradient step (one-step) or a maximum likelihood submodel approach, we explicitly specify that the solution should minimize empirical
7

---

Published in Transactions on Machine Learning Research (03/2023)
X
Q-method
Outcome Model
U-method Update Step
G-method
X
Propensity Score
This figure illustrates the components involved in using IFs to improve our estimates of the Average Treatment Effect (ATE), where the ATE is our target estimand V. We combine the output from an outcome model Q with a propensity score model G and an update step method U. This yields an estimate Vb
Figure 2:
approximations (circumflex/ ceruumflex/haat hat notation) of both (A) the expectation and (B) the variance of the influence function. The degree to which each of the constraints are enforced depends on hyperparanelers an E R+ and a2 R+ which weight the two constraints. In this objective, the parameters y are related the influence function as follows:
θATE(Z, P,) H(z,, y: xtt,,x)) - -x)-((()x)--((11 +ll(z)) +(Q(1,x, + (2, 1)) - (Ô(0) x,) + +HH( 0)) ATE
(6)
In other words, y is the coofficient on the elever covariate which itself represents the quantity which is biasing our principal estimator Q. The objective in expression 5 therefore finds the f which corrects for the biasing action of H by minimizing the mean and the variance of the influence function, which itself is computed according to Eq. 6
One of the consequences of finding the efficient IF is that also achieve improved model robustness. This is because, when the update step is used to derive unbiased estimate, we achieve consistent estimation (ie., we converge in probability to the true parameter as the sample size increases) even if one of the models is misspecified (e.g., the ATE requires both a propensity score model and an outcome model, and thus the facilitates double robustness). Furthermore, in cases where both models are well-specified, we achieve efficient estimation. It is worth noting, however, that this double-robustness property does not apply the limiting distribution of the estimates being Gaussian when data-adaptive plug-in estimators are used [Benkeser et al., 2017 van der Laan 2014). In other words, if only incorrectly specified, the estimates may not be normally distributed, thus invalidating statistical inference. In our later evaluation, we thus might expect models to fail at achieving normally distributed estimates before they fail at yielding unbiased estimates. It is possible to extend the framework such that the double robustness property also applies to the limiting normal distribution the estimates (Benkeser et al. 2017 van der Laan 20144. but we leave this future work. For more technical dotails on the double robustness property see van der Laan & Rose (2011): Hines al. (2021); Benkeser et al. (2017), and Kurz (2021).
both of the two models is/are
5 The MultiNet Estimator
This section concerns our second contribution a proposal for a new estimator for the ATE, 'MultiNet', which aims to compete with existing estimation methods described Scction 2 It represents a combination of ideas from two well-known approaches in causal inference- the Super Learner (van der Laan el. al. 2007 and CFR-Net, both of which already mention above.
5.1 Motivation
One of the primary considerations when choosing estimation algorithns/models is whether the estimator can represent family of distributions which is likely contain the true Data Generating Process (DGP). Indeed, of the motivations for semiparametries is to be able use non-parametric data-driven algorithms
8

---

Published in Transactions on Machine Learning Research (03/2023)
which have the flexibility model complex DGPs, whilst still being able to perform statistical inference. If the estimator is functionally misspecified (i.e., misspecified in terms of the functional form used to model the relationships in the DGP- - linear, quadratic, spline etc.), then we are unlikely to arrive at an estimator which is asymptotically linear and therefore also amenable to the Influence Function update process. This behooves us to seek estimators which 'let the data speak' 4van der Laan & Starmans 2014
Consider the Super Learner (SL) yan der Laan et al. 2007), which is an ensemble method especially designed for parameter estimation in the context of causal inference. The SL process involves taking a weighted average of predictions from each candidate learner, and this quantity is taken the output. The advantage of SL is that the candidate library includes sufficient diversity with respect to functional form and complexity such that the true DGP is likely to fall within the family of statistical models which can be represented by the ensemble. The motivation for reducing bias resulting from functional misspecification is therefore similar to the motivation for influence functions. In both cases, accuracy / precision in estimation is the priority in the domain of causal inference, where the parameters concern critical decision making processes, such as the efficacy of medications.
In contrast with the SL, many methods (including boosted trees, neural networks, linear regressions) are based on single learners with one outcome prediction. As part of the development and evaluation of Influence Function updating methods for reducing bias and deriving consistent. estimators, here we also present a new neural-network based estimator / learning algorithm. Early experimentation highlighted to us that even though NNs are flexible universal function approximators (Hornik) 1993; Hornik et al. 1989) they may nonetheless yield estimators which are not 'good enough' to enable to leverage their asymptotic properties (such as bias reduction with IFs). In such cases, the IF update may actually worsen the initial estimate, pushing us further off course. This problem arose even for simple datasets with only quadratic features. Indeed, the problem with using neural networks for 'tabular' data (as opposed to, say image data) is well known in the machine learning community, and interested readers are directed towards the survey by Kadra et al. (2021). Rescarchers have, in general, noted that gradient boosted trees (Fruund & Schapire
1997) to consistently outperform neural network based learners (Shwartz- & Armon 2021 Kadra el al.
2021 Borisov et al.] 2022). However, Borisov et al. (2022) also found that ensembles of boosted trees and neural networks can nonetheless outperform boosted trees alone, and Kadra et al. (2021) found that sufficiently regularized neural networks could yield competitive performance, even exceed the performance of boosted trees. Thus, in our view the avenues for research into neural network methods for tabular data are still open (and research on the subject continues regardless- - see TVAE, CFR-net, and DragonNet, for example). Furthermore, neural network based methods work well in ensemble combinations with boosted trees, we should attempt to maximise the performance of the neural network learners in order to maximise the performance of the associated ensemble.
5.2 Details
A block diagram for MultiNet is shown in Fig. 3 The method comprises three main elements: (1) a CounterFactual Regression (CFR) network backbone (Shalit et al. 2017). (2) outcome prediction 'taps' at each layer, thereby emulating the outputs from multiple discrete learners in a traditional ensemble and (3) a constrained regression procedure designed to reflect the equivalent 'meta-learner' idea in the Super Learner and which combines the intermediate predictions from each of the layers. We discuss each of these three clements in turn below. The idca therefore represents a combination of two well-known methods in causal inference- CFR-net, and Super Learner.
5.2.1 CFR Backbone
CFR is popular NN method for causal inference tasks comprised of fully connected layers. It includes separate outcome arms depending on the treatment condition, and forms the backbone of MultiNet. Specif- ically, for each consecutive layer in MultiNet, we predict y/4,x for L = {0. 1} and compute the corresponding layer-wise cross-entropy loss (for a binary outcome). This simulates the multiple outputs of a typical en- semble method each layer represents a different (and increasingly complex) function of the input, e.g., for network with L, layers, the predicted outcome 9 under treatment. 1 for the 7t layer of the network
9

---

Published in Transactions on Machine Learning Research (03/2023)
MultiNet
[ ===a. [xx]
x xx) 1'
Constrained Regression
YT
((xxxyr((xx)) = -(NN1))
X
jull
Milla
Figure 3: Block diagram for MultiNet. At each layer 1 = (1,..., L} the network, the outcome y is estimated using covariates x (which can include treatment 1). The treatment is used to select between two estimation arms. Once the network has been trained, the outcomes from each layer are combined and a constrained regression is performed. The weights / in the regression constrained to be positive and sum to 1. An equivalent single-headed network can be used for the treatment model f|x.
is (1) = qf (a/ (a0=(((/22xxt 1)))) where of represents the 7h layer of the network for the CFR arm when treatment is 1. The equivalent for the outcome under treatment is thus @r( of (ff- ( 0)))) In our implementation,
the ReLU activation function.
5.2.2 Ensemble Behavior
Collecting the layer-wise predictions for the outcome vi to UL together form matrix Y which has shape (L x N) where L is the number of layers and N is the number of datapoints. This idea stemss from the one presented as part of the SuperLearner. By consequence of the diversity afforded by multiple learners and the 'meta-learning" ensemble combination of the predictions from these learners, the Super Learner is able achieve efficient rates of convergence and simultaneously reducing the chances of overfitting.
5.2.3 Ensemble Meta-Learner / Constrained Regression
The constrained regression combines the predictions accumulated from each layer of the 'ensemble', and is only applied after MultiNet has been trained. For each observed treatment condition in the training dataset, we train a meta-learner by solving YT g = with layer-wise weights B which are constrained to sum to one and be neg-negative (similar to a non-negative least squares objective). These weights therefore determine how much importance is given to each layer in the network. We expect, for instance, that for linear datasets the early layers will assigned a greater weight, whereas for more complex problems, predictions from deeper layers will become more important. Indeed, we observed this behaviour during development, but leave detailed exploration of this to future work.
For solving the constrained regression we use SciPy (Jones et al. 2001 non-negative least squares solver. The weights are then used to create weighted averages of the predictions from the consecutive network layers for subsequent predictions, emulating the SuperLearner-styv ensemble approach with single network. Note that one the strengths of this approach is that the layer-wise wise outputs and constrained regression tochniques can be flexibly applied to other neural network architectures.
5.2.4 Variations
In order to explore the possibility for increased performance, we explore some additional variations on the idea. Firstly, allow each layer-wise loss gradient to influence all prior network parameters. This
10

---

Published in Transactions on Machine Learning Research (03/2023)
is similar to the implementation of the auxiliary loss idea in the Inception network (Szegedy et al. 2015 and we refer to this variant 'MN-Inc'. The second variant involves only updating the parameters of the corresponding layer, preventing gradients from updating earlier layers. We call this variant the 'cascade" approach, and refer to this variant as 'MN-Casc'.
Finally, order to increase the diversity across the layers and to approximate the diversity of an ensemble, we also explore the use of loss masking. For this, we partition the training data such that each layer has a different 'view' of the observations. The loss is masked such that each layer is trained on different, disjoint subset of the data. We refer to variants of MultiNet with loss masking 'MN+LME
5.2.5 Objective Function
The objective function of MultiNet is:
n
L m
C min n ΣΣ T Σ
m
(7)
where m is the mask for datapoint i in layer I, and C2 is the cross-entropy loss for datapoint i, predicted by layer 1. Note that, according to the CFR backbone, there are two prediction arms for each layer, one for when treatment t = and one for when t = O. In the variants without loss masking, we set m in the objective function to be equal to for all layers and all datapoints.
It is worth noting that other ensemble neural network approaches exist, such as the Bntebbe approach by Wen et al. (2020), and the well-known use of dropout Gy]Ga& & Ghahramani (2016). Whilst incorporate dropout into MultiNet's 's training procedure, we do not use the dropout itself as means to create an ensemble, although this represents an interesting opportunity for further work. We instead explore 'tapping' cach subsequent layer in neural network to collect intermediate outcome predictions, and combining them with constrained regression the manner of the Super Learner. We indeed observe improvements in performance over the CFR-Net backbone, highlighting that this simple modification can provide enough diversity to yield competitive, ensemble-style performance.
6 Evaluation
This section concerns third contribution comprehensive evaluation of the performance of estimators and influence functions for the task of ATE estimation. The main results of this evaluation are given in and the results of a Shapley value analysis (Shapley 1953 Lundberg et al. 2020) are presented in Additional experimental results and discussion can be found in the appendices.
Sec. 6.3
Sec.
6.4
6.1
Open Questions
So far, we have presented the relevant background for causal inference and IF's proposed a new MultiStep update process and proposed a new NN based estimator called MultiNet. The following open questions remain: (1) Can estimation methods be improved using the one-step, submodel, MultiStep (ours), or tar- geted regularization approaches? (2) How do various different outcome, propensity score, and update step methods compare? We aim to answer these questions through extensive evaluation of different met (Sections 6.3 and 6.4h. In particular, (a) precision in estimation, (b) robustness, and (c) statistical inference (normality of the distribution of estimales). We use these open questions to inform the design of our experiments, which are described below.
examine the performance of the different approaches in terms of
6.2 Experimental Setup
6.2.1 Methods, and Evaluation Criteria
We evaluate a number of different methods in terms of their ability to estimate the ATE. summary of the complete set of methods explored as part of the evaluation is shown in Table [ As described above, we
11

---

Published in Transactions on Machine Learning Research (03/2023)
are interested in three properties relating to performance: estimation precision, robustness, and normality, Estimation precision is evaluated using mean squared error (MSE) calculated as r Ef - =| where = 100 is the number of simulations, and the standard error (s.e.) of the ATE estimates is computed as the standard deviation of é. Note that this is not the same as estimating the precision in the Conditional ATE, because the indexing occurs over simulations (í.e. T is an average treatment effect, estimated for an entire empirical sample), not individual participants. The MSE was chosen for its scnsitivity large errors (arguably important in the application domain of causal inference), and because was found that it provided more informative spread of results than the root-MSE. Robustness will be evaluated by comparing initial estimators that fail to exhibit the desired properties, with the results once these estimators have been updated. For normality, we examine the empirical distribution of the estimates. Using these distributions, we provide p-values from Shapiro-Wilk lests for normality (Saapiro & Wilk 1965). Doing so provides an indication of the estimator's asymptotic linearity and whether the IF's are facilitating statistical inference as intended.
6.2.2 Data
Recent work has highlighted the potential for the performance of modern causal inference methods to be heavily
and has recommended the use of bespoke datasets which transparently test specific attributes of the evaluated models across different dimensions (Curth et al.] 2021b). We therefore undertake most of the evaluation using variants of a DGP which we refer to as the LF da-dateset and which has been used for similar evaluations in the literature (Luque-Fernander el al.] 2018). We also evaluate using the well-known IHDP dataset (Hill 2011 Dorie 2016), as well a generalized version of the LF-dataset.
Figure 4: Graph for the 'LF' dataset used by Luqque-Fernanderr et al. 2018
LF Dalasel Variants: The initial and original LF-datasel variant, (v1), models 1-year mortality risk for cancer patients treated with monotherapy or dual therapy. One motivation for starting with this DGP is that its polynomial functional form is not sufficiently complex to unfavourably bias the performance of any method from the start. The dataset also exhibits near-positivity violations, and will therefore highlight problems associated with the propensity score models which are necessary for the update process. We also adjust the level of non-linearity in order to assess the robustness of each method to increased complexity. Accordingly, introduce exponential response into the potential outcome under monotherapy (t [t=]] for the second variant. (v2). Our I.F-datasels comprise 100 samples from set of generating equations. Both variants are designed to highlight problems which may arise due to near positivity violations.
The graph for the synthetic 'LF' dataset used in work by Luque-Fernandez Fernandez et al. (2018) is given in Fig. 4 The DGP is based on a model for cancer patient outcomes for patients treated with monotherapy (t = T) and dual therapy (t = 0). We create two version of this datasct (v1) and (v2), where the data generating process of the second variant includes an exponential non-linearity in the outcome model. The two variants are designed to yield near positivity violations in order to highlight weaknesses in methods which depend on reliable propensity score model.
We also create a generalized version of the LF (v2) datasel which we refer to as 'General', or 'Gen' for short, in the experimental results. The functional form underlying the key relationships in this data generating process are the same LF (v2), except we include more structural complexity. The corresponding graph can be found in Fig. 5 Specifically, we inchude collider structure (C), a mediator (M), two instrumental and two risk variables (T,, 12, R1, and R2, respectively), and new confounder substructure (T Xs X1 + Y).
12

---

Published in Transactions on Machine Learning Research (03/2023)
Q Mothod
Method
U Method
Datasets
Evaluation Criteria
Linowr dinenr/logist Regrosson (QLIR)
Liurar; Regrossion (G-LR)
OoeStpp (U-ones)
[yl] (-00000 10000)
Sepwoe Peror (MSR)
SuperLearuur (Q-SL)
Supeel.ummet (G-SL)
Submodel (U-eab)
(v2) n=(500, 3006, 10000)
Stapiro Wilx Test (p)
CER (Q-CFR)
CPR (G. CFR)
MaltiSepp (U munts))
Gen u- (500, 5000, 10000)
Stondard Error Estunattion (s.c.)
MultiNet MN) vuriante
MulsiNes (G MN) variants
Targeded Regularization
IHDP n-747
(QTVAE)
P-kauuer (G-P)
(U-Base)
DragonNet (Q-D)
Dragie Not (G-D)
Slearne (Q-S:
T-eamer (Q-I)
DML (DML.
Table I: A summary of all variants and metrics explored as part the evaluation. Q-methods concern the model for the outcome, G-methods concern the propensity score models for the treatment assignment, and U-methods concern the manner in which the influence function is used to update the initial estimates. Note that the Double Machine Learning (DML) method is treated its own, without combination with G- and U-Methods.
For estimation our adjustment set comprises all the X variables, even though it is possible (a) to include the risk variables, and (b) lo omit either X5 X1.
For LF (v1), LF (v2), and Gen, we create further variants with different sample sizes n {500, 5000, 10000} in order to explore sensitivity to finite samples. The gencrating equations for these variant datasets (including propensity scores for and v2) are given in Appendix
IHDP: The final dataset comprises 100 simulations from the well-known IHDI 15 dataset. We use the version
corresponding with usual setting A of the NPCI data generating package Dorie 2016 (see Shi et al. 2019
Shalit et al. 2017 and Yao et al. 2018) and comprises 608 untreated and 139 treated samples (747 in total). This variant actually corresponds with variant B from Hill (2011). There are 25 covariates, 19 of which are and the rest are continuous. The outcome generating process is designed such that under treatment, the potential outcome is exponential, whereas under no treatment the outcome a linear function of the covariates (Curth et al. 2021b). This dataset represents a staple benchmark for causal inference in machine learning. However, is worth noting that recent work has shown to preforentially bias certain imaioos (Curth et al. 2021b) so include this dataset for completeness but discount our interpretation of the results accordingly.
X
X2
11
R₁
X
X
Ra
M
Figure 5: A generalized version of the LF dataset which we refer to as the 'Gen' dataset in the experimental results.
6.2.3 Algorithns/Estimatoos
The training details and hyperparameters are presented in Appendix F compare linear/logistic regression (LR); Super Learner (SL) comprising LR, LR with extra quadratic features, a Support Vector classifier, a random forest classifier [Breiman 2001 a nearest neighbours classifier
For the outcome model Q we
Available from http:://owwfrreeo.ccom
13

---

Published in Transactions on Machine Learning Research (03/2023)
AAltman 1992 and an AdaBoost classifier (Freund & Schapire, 1997p; an implementation of the backbone to CounterFacturey Regression network (without the integral probability metric penalty) (Shalit et al.) 2017 (CFR); DragonNet (D) with and without targeted regularization (Shi et al.) 2020): TVAE (Vowels et al,
2021 (which includes targeted regularization); T-lcarner (T) Kunzel et al. 2019) with a gradient boosting
machine (Friedman 200T); S-learner (S) (Kunzel et al., 2019) with a gradient boosting machine (Friedman
2001}; and our MultiNet (MN) variants (MN-Ine, MN- Case, MN--In++ L.M, MN- L.M When estimating
the IF of the ATE, we also need estimators for the propensity score treatment model, which we refer to as G. For this we use LR. and SL, ElasticNet 'P-learner' (Zou & Hastie 2005), DragonNet, as well as CFR and MN. The latter two NN methods must be modified for this task, and for this we simply remove one of the outcome arms, such that we can estimate t/x. This seletion of algoritms was chosen to represent a suitable diverse set of common yet. weak learners (e-g-, logistic regression), modern well-knoww neural network causal inference methods (e.g., CFR-net), neural network methods which already incorporate semi- sami-parameteric techniques (DragonNet and TVAE), and recent and/or or popular methods proposed in the causal inference or biostatistics literature (e.g. Super Learner, T- and S-learners). Finally, we also include the Double Machine Learning (DML) method Chernozhukov et al.. (2018) implemented in the DoubleMLIRM class in the DoubleML python library Bach et al. (2022). The DML estimator uses stacking classifier comprising LR, random forest, nearest neighbours algorithm, support vector machine, and an AdaBoost algorithm (and is thus similar in its adaptability to the Super Learner, as specified above). Note that the DML, method does not need to be 'updated' (it is already a doubly-robust estimator), and as such we consider it separately and provide results for DML, and compare it with the other methods, in the appendices.
6.2.4 Update Steps
We evaluate the onestep (U-ones), submodel (U-sub), MultiStep (U-multi), and targeted regularization (TTeg) approaches to the update process.
The MultiStep update variants are optimized using the Adam Kingma & Ba 2017 optimizer. For small datasets (n < 1000) we undertake full gradient descent (i.e., using the full data), and for larger datasets we use stochastic mini-batch gradient descent. The batch size for datasets with a sample size n > 1000 is set to 500, we undertake 4000 steps of optimization, and the learning rate for the Adain optimizer is set to 5 10-4, The MultiStep objective has hyperparameters aj and a2 which weight the constraints in the objective (expectation and variance of the influence function, respoctively). We set both to onc.
6.3 Results
Given the large number of combinations in a full-factorial design (approximately 5000 results), we undertake an initial set of experiments to narrow down the evaluation space to focus on the most competitive mothods. With this 'shortlist', we investigale the contribution of each Q-, G-, and U-method across the 7 different dataset variants. The results for this initial evaluation are presented in Appendix H] In summary, we identified that CFR did not perform sufficiently well to warrant further investigation (note, as mentioned above, that our implementation of CFR does not utilize the full objective function presented by the original authors, only their proposed architectural foundation). Furthermore, the best performing MN variant was MNN and we use this variant for the subsequent analyses. Finally, targeted regularization was inconchsive. However, previous work has identified its potential to improve DragonNet and TVAE Shi et al. 2020 Vowels et al. 2021) and so we restrict the application of targeted regularization to these methods only, in the main evaluation presented below. Note that even following this initial evaluation / shortlisting process, the number of results remains large. We have allempted to summarize them in Figs. but
include further results in the Appendix in Table 4 Table and Figs. 20
Whilst it is possible and potentially helpful to simply present the full set of results, the overwhelming number of combinations does little to help us understand whether the use of particular Q, G- or U-methods are more or less likely to improve or worsen the performance in any particular combination on any particular dimension/metiio /ectric. Therefore, Figs. K 7 and 8 provide results for P(O)M) = p( (M ())()))p LF and General dataset variants for MSE and s.c. Here, M is the method, and O is the quantile (we split into 5 quantiles) for MSE and s.e., respectively. In words, the associated plots provide estimation for the
across the
14

---

Published in Transactions on Machine Learning Research (03/2023)
probability of achieving performance result in each quantile O, for given method M, thereby providing means to directly assess the relative performance of each Q- G-, and U-method. For instance, we can split the MSE results into equal probability quantiles, and count the number of times the use of each outcome, propensity score, and update method results in a performance which falls into each of these quantiles. Using Bayes rule we gel an estimale for the probability of achieving results in a particular quantile (e.g., the best performing methods fall in the zeroth quantile of MAE results), given a particular choice of method. Using these calculated probabilities, we also select all results from the best quantile, and see how the performance shifts over different sample sizes. Note that because these results are based on rank ordering, it is not possible to judge absolute performance, only relative performance. Indeed, the purpose of the initial results above was to use the absolute performance as a way of shortlisting the methods so that a more comparative evaluation could be undertaken using the more competitive methods.
To evaluate the normality of the estimates, after calculating the p-value from the Shapiro- Wilk test, we calculate the proportion of each Q-, G-, and U-methods which yield normally distributed estimates (p > 0.01). For example, if a particular Q-method has high 'probability of normality' according e.g. Fig. 10 this means that large proportion of the results yielded normally distributed estimates.
In the following Sections 6.3. 1..1 6.3.7 we review the performance of each method for each of the three perfor- mance metrics in turn.
When interpreting the results shown in Figs. G 10 it may be useful to recognise the 'ideal/ desired' curve as one which takes on high values on the reft-hand-side representing yielding results in the top quantile(s)- In contrust, curves which take on high values on the right hand side represent those with high probability of yielding poor results in the lower quantile(s), (s), relative to the other methods.
Note:
high--roobbbllte of the method
6.3.1 Q-Methods - MSE
Beginning with Fig. 6 the results for the outcome model Q-methods on the LF dataset variants are shown in the first column. In Fig. Gh we see that our Q-MN achieves the highest probability of being in the best quantile for MSE when used as an outcome model Q for LF (v1) n = 500, followed closely by Q-LR and Q-SL, and Q-TVAE and Q-D in the second-best quantile. In contrast, Q-D without targeted regularization, Q-T, and Q-S all had higher probabilities of yielding results in the later quantiles (i.e., their performance was worse). Increasing the sample size to 5000, and considering Fig. G1, we see similar results, with MN again yielding the highest probability of the achieving the best results, with Q-D, Q-S, Q-T, and Q-D without targeted regularization performing the worst. Finally, for LF (v1) n = 10000, we see in Fig, Ge that Q-MN is superseded by Q-LR and Q-SL, followed by Q-TVAE. Q-T, Q-S, and Q-D perform poorly again.
These results suggest that Q-LR and Q-SL perform consistently well over different sample sizes, and that Q-MN can perform well in small sample sizes, but may start to overfit as the sample size increases. Recall that the task of causal inference is different from the typical supervised learning task, and more data does not necessarily imply that it is casier to estimate the difference between two response surfaces, particularly when this difference (which is the treatment effect) is of low-complexity relative to the response surfaces themselves.
Now consider Figs. i. m. p) for LF (v2), which introduces additional non-lincarity into the outcome model. We initially observe similar results for n = 500 in Gj with Q-MN, Q-L.R, and Q-ST. achieving the best results, and Q-S, Q-D without targeted regularization, and Q-T populating the later quantiles. Incrcasing the sample size to n = 5000, we see in Fig. 6n that Q-TVAE now becomes the most likely to yield the best results, followed by Q-SL and, interestingly, Q-D without targeted regularization. Q-D, Q-T, and Q-S, however, still perform poorly. Finally, for n = 10000, we see Q-TVAE maintain the lead, once again followed by Q-SL. The worst performers were, again, Q-D, Q-S, and Q-T. This suggests once again that Q-SL provides consistent performance across sample sizes, and that Q-MN is a good option for smaller sample sizes.
For the Gen dataset, we see from Figs. B(a.d.)) that the results are again quite mixed, although there is some indicate that Q-LR, the logistic regressor, performs reasonable well, and, for instance, that Q-TVAE performs poorly as the sample size increases.
15

---

Published in Transactions on Machine Learning Research (03/2023)
Q-methods
G-methods
U-methods
(a)
(b)
(c)
0.8
ONN
0.Dnotrep
GAN
GDnotreg
U-mijt meanver
U-ores
Q-AL
GBL
LF (vr) n=500
0.6
0.2
9,8
(d)
0.D
(
GD
i
U-Sup
pmotreg
Gnncteeg
O.S.
TYME
(v1] (v1)0=00000 n=5000
GLR
0.4
0.2
0.0
1.0
(g)
(h)
0.8
OINN
Dmooreg
G.Dntrreg
GAN
U-m Bn-nneeee
0.50
G-SL
0.6
O
GP
(v1) n=11
O-R
GiR
0.4
Probabilty of MSE quantile given methodology
0.2
0.0
Q.D
(k)
GD
U-su3
G-Dnotrip
GDnotren
LF (V2) n=500
0.6
OETMAE
GLR
0.4
0.0
(m)
(n)
(o)
0.8
0PHN
Dmotreg
GAN
O.S.
GSL
Bade
0.6
THAE
GP
(v2) n=5000
GLR
0.4
0.2
(p)
(q)
(r)
GU
Disus
0.8
DeDroer
G.ncttreg
GAN
U-m IE Bn-mmeerver
U-ones
(v2) n=10000
0.8
OUTWAE
OAR
0.5
0.2
0.0
2
Quantile
Figure 6: After recording the MSE for each Q (oulcome), G (propensity score), and U (update step) method combination, we rank order them (from lowest to highest MSE), and calculate p(O]M) where O the MSE quantile, and M the method. For quantiles, this enables us to find e.g., the probability of getting a MSE in the best quantile given a particular method p(O 0JM m). If method performs well, we expect have high probability of achieving MSE in the top two quantiles.
16

---

Published in Transactions on Machine Learning Research (03/2023)
For the IHDP dataset, we use fixed sample size of n 747, and the results shown in Fig. 9 Here it can be seen that Q-T and Q-TVAE achieve the best results, followed by Q-S and Q-MN. The worst performer was Q-LR. These results are consistent with previous work which highlighted state-of-the--rt performance of TVAE on IHDP (Vowels et al. 2021 Similarly, the fact that LR did so poorly possibly highlights the non-linearity of the data generating process for IHDP. The fact that Q-S and Q-T did so well is surprising given their relatively poor performance on the LF datasots described above. Such dataset dependence for the performance of causal estimators has also been previously noted by Curth et al. (2021b).
6.3.2 G-Methods MSE
The MSE results for the propensity score G-methods can be seen in the second column of Figs. and g and in plots (b, e and h) in Fig. Interestingly, there is very little dependence between the performance of the different methods. Arguably, there is some evidence that G-MN performs slightly worse than other methods in Fig. 1 and that G-D performs worse in Fig. D but the differences are not convincing. For instance, G-D performs quite well in Fig. g (b) on the 'Gen' dataset when 500, but not noticeably well as sample size increases. This suggests that, at least in our experiments, the MSE results are relatively robust to the
choice of propensity score model.
6.3.3 U-Methods MSE
The MSE results for the update U-methods are shown in the third column of Fig. 0 for the LF datasets. In Fig. D we see that the U-Base model and the U-multi update methods perform the best, with the U-ones model close behind. The submodel update is more likely to be the lower quantiles. As the sample size increases to n 5000 and n = 10000 in Figs. B and G we see the U-sub and, to a lesser extent, the U-ones performance shift. This behaviour has been observed before in work by [Neugebauer van der Laan (2005); who found that the performance of U-ones increased with sample sizes. Indeed, their proposition for multistep update process also performed more consistently in small samples, does our U-multi. Similar patterns of performance are seen in Figs. 6 6 and 6 for the LF (v2) dataset.
In Fig. we again see mixed results on the 'Gen' dataset. Whilst, in Fig. 9 there is some indication that the U-sub and U-ones performed approximately equally well and better than U-multi and U-Base on the IHDP dataset, the results are still far from indicating a clear pattern in performance.
6.3.4 Q-Methods s.e.
The standard error (s.c.) results are shown in Fig. and the bottom row of plots in Fig. 9 Starting with Fig. 7p, we find the methods yielding the tightest distributiin of estimates for the LF (v1) dataset n 500 are Q-MN, Q-LR, and Q-TVAE, followed by Q-D, Q-SL, and Q-T. At the lower end we find Q-D without targeted regularization, and Q-S. As the same size increases to 5000 Q-MN provides estimates which are even more likely to be the tightest, followed again by Q-LR, Q-TVAE, and Q-SL. Q-D is not far behind, with Q-S, Q-T, and Q-D without targeted regularization performing the worst. With n 10000, Q-MN is overtaken by Q-LR in terms of the tightness of the estimation, which is understandable given that Q-MN has a large mumber of hyperparameters (Q-LR has none), which contributes to variability in performance. Q-TVAF once again follows closesly behind, with the worst. performers being Q-D without targeted regularization, Q-S, and Q-T. Interestingly Q-MN exhibits a rise in the probability of being one of the worst performers, suggesting that there may exist better or worse combinations of G- and U-met with Q-MN. Once again, it is worth consulting the full set of rank-ordered results in the Appendix. With the results for LF (v2) in Figs. 7) 7 m, and 7p we see a similar pattern of results, in spite of the introduction of additional non-linearity in this dataset variant.
There some indication in Fig. ((jm that Q-SI and Q-LR perform well, whilst Q-D and Q-T have a higher chance of yielding results in the less competitive quantiles. However, the results are quite mixed, once again. Finally, for the IHDP results in Fig. D we see Q-LR and Q-SL provide the tightest estimates, followed by Q-MN, Q-D without targeted regularization, then Q-TVAE, Q-S, and Q-T.
17

---

Published in Transactions on Machine Learning Research (03/2023)
Q-methods
G-methods
U-methods
1.0
(a)
O-b
GD
(c)
O.Dnetreg
Gonatreg
Ulones
0.8
0.5L TWAE
U-Bado
0.6
LiL
Q
1500 LMA
0.4
0.2
0.0
1.0
(c)
(
O-Dnotreg
G-Dnotreg
U-ones
0.8
OPNN
GPNN
0.5L
G.S.
U-Bade
0.8
TTME
G-P
LF
0.4
n=5000
0.2
0.0
1.0
Q-D Diootrng
(9)
(h)
GD
Usub
G.nnttreg
Uinnes
0.8
Probability of SE quantile given methodalogy
O.S.
U-Bac2
QUI
GLR
(11)) m= 100000 10000
0.4
0.2
1.0
(k)
0-0
Dnotroo
G.nobrop
Ulnes
0.8
ON 0-5L
GNN
Uimith In meaewwr
L-Baae
0.6
22
7=500
0.2
0.0
(m)
O.D
GD
U-uub
O.Dnsoroog
G-Dncbrop
U-onés
U'muli
0.8
0.4
155000
0.0
1.0
(p)
0-D
(q)
GD
(1)
Douh
0.Drotrcg
G.Dnetreg
U-onés
0.8
U-Bace
O
200 1 100000000
QULI
0.4
0.2
Quantile
Figure After recording the standard error (s.c.) of the 100 ATE estimates for cach LF datasct and for each Q (outcome), G (propensity), and (update step) method combination, we rank order them (from low to high), and calculate p(O)M) where O is the quantile, and M is the method. This enables to find the probability of getting s.e. in the best quantile given particular method p(O OJM m). If method performs well, we expect to have high probability of achieving an s.e. in the top two quantiles.
18

---

Published in Transactions on Machine Learning Research (03/2023)
(a)
(b)
(c)
C-Onatre
GAN
G.LL
QHR
U-TWKE
Gen m0:000
Q-1
Probabiity of MSE quantée given methodology
0.0
(8)
(f)
Usi
Beokog
OWNN
G-SL
Gen n0000
0.8
0.2
9.8
(g)
(h)
(1
O.Dratreg
Usdrne
0.5L TIME
Uibase
Gen non1000000
0.4
2
Quantile
()
(K)
())
0.8
GAN
0.5L
G-L
O-TVAR
GER
Gen n=:00
Probability of quantile given methodology
0-S
1.0
(m)
0.0
(n)
GD
G-Orobep
(o)
Douh
0.8
MN
GAN
G-LL
U-Eate
TIME
0.6
Gen n=5000
0.4
9.8
(p)
q)
(r)
OUTONE OILR OCT
Gon 11:
DAS
Quantile
Figure 8: After recording the MSE and standard error (s.e.) for each Q (outcome), G (propensity score), and U (update step) method combination the 'Gen' dataset variations, to highest MSF),, and calculate P(O where O is the MSF quantile, and M is the method. For 5 quantiles, this enables us to find e.g., the probability of getting a MSE or s.e. in the best quantile given a particular method p(O M m). If a method performs well, we expect to have high probability of achieving an MSE in the top two quantiles.
rank order them (from lowest
19

---

Figure 9: After recording the MSE and standard error (s.e.) the 100 ATE estimates for the IHDP dataset and for each Q (outcome), G (propensity score), and (update step) method combination, we rank order them (from lowest highest), and caleulate M) where O is the MSE (top row) or s.e. (bottom row) quantile, and M is the method. For 5 quantiles, this enables to find the probability of getting MSE s.e. in the best quantile given particular method p(O O|M m). If method performs well, we expect to have high probability of achieving an MSE and/or s.c. in the top first or second quantiles, and a low Published in Transactions on Machine Learning Research (03/2023)
Q-methods
G-methods
U-methods
Probabiitt of MSE quantlee
1.0
(a)
(b)
(c)
DDnottteg
Gineatoog
given methodology
0.8
CHN
Uemahalis maanoar
O SL
GSL
Bade
10/-- HHHP
OUR
GLR
0.4
0.0
Probability of SE quantile
1.0
(d)
0n
GD
(f)
Dauhh
given methodology
0.8
O Dnotreg
OHN
G Gnettrrp
GAN
U-mul--Immmeeree inmmtstwvr
Unnes
OSL
G-L
U-Basc
1 dOH dOHI
D-TWAE
0.4
0.0
0
probability of achieving an MSE and/or s.e. in the last quantiles. Best viewed in colour.
6.3.5 G-Methods s.e.
The s.e. results for the choice of propensity score G-method can be found in the central columns of Fig. and Fig. (0,, and Figs. gk.n,q). As was found for the MSE results, the choice of G-method was not decisive, besides the poor performance of G-MN for IHDP dataset, and for the n 10000 LF datasets. It is interesting to again find that the choice of G-method does not have a strong impact on the tightness of the estimates.
6.3.6 U-Methods s.e.
The s.e. results for the choice of update U-method are presented in the right-hand column of Fig. and Fig. g.. In contrast the choice of G-method, the choice of U-method had a significant impact on the tightness of the associated estimates, and the pattern of performance is similar to the pattern for MSE. For low sample sizes, it can be from both Figs. 7 and that the tightest estimates are achieved using U-multi and U-Base, with U-sub yielding the least tight estimates. Inereasing the sample size shifts the performance of U-sub and U-ones, making them competitive with the other methods.
For the 'Gen' dataset, Fig. Lo,r) indicate that the base models perform well, whilst U-sub seems to be one of the least good performers. For the IHDP datasot, can be seen in Fig. D that the choice of U-method had little impact the tightness of the estimates, but the best performers were U-Base (f.e., no update), and U-multi.
6.3.7 Q- G-, U-Methods Normality
The results evaluating the normality of the estimates are provided in Fig. 10 for the LF and Gen dataset variants, and Fig. [ for IHDP. For the LF and Gen datasets, plots (a-i) provide the proportion results from the respective method which yielded normally distributed estimates (p > 0.01) for each of the different dataset sizes (500, 5000, 100003. In Fig. 10 can be seen that most Q-methods performed well across all sample sizes with LF (v1), with the exception of Q-D which was less likely to yield normally
20

---

Published in Transactions on Machine Learning Research (03/2023)
distributed estimates, and we observe a drop in performance for Q-MN as sample size increases. Once again, and as indicated by Fig. T0p, the choice of G-method was not found to impact the likelihood of normally distributed estimates. Figs. [0] indicates that the likelihood of U-Base and U-multi yielding normally distributed estimates dropped slightly with sample size, with U-ones yiclding consistently normally distributed estimates regardless of sample size.
For LF (v2), the results in Figs. LO indicate more variability, possibly as result of the additional
non-linearity in the outcome model. When n 500, the outcome Q-method most likely to yield normally
distributed estimates was Q-T, followed by Q-D and Q-MN. However, for 5000 and n 10000, the only methods not yiclding consistently normally distributed results were Q-D and Q-MN. For the propensity score G-methods, the method most likely to yield normally distributed results with n = 500 was G-LR, followed by G-SL. The other methods did not perform well until the sample size was increased to n = 5000 or n = 10000 for which all methods performed equally well. For the U-methods, the best performing result across all sample sizes was U-ones, followed by U-sub, U-multi, and finally U-base.
For the Gen dataset, the results in Figs. L0 10 indicate that Q-MN and Q-D were two of the methods least likely to yield normally distributed estimates, particularly as the sample size increased. On the other hand, the G-LR performed well as a propensity score method in this regard, and our MultiStep worked well as an update method.
Finally, the likelihood of achieving normally distributed estimates for the IHDP dataset are shown in Fig, II The sample size is fixed for this dataset, and the results for the Q- G-, and U- methods are presented together (hence the different graph format). Tt can be seen that Q-D provided the highest likelihood of normally distributed estimates, with the other methods yiclding comparable (and low) likelihood. Similarly, G-D yielded the highest likelihood of normally distributed estimates, with the other G-methods being relatively equal (and low). Finally, none of the U-methods provided a high likelihood of normally distributed estimates.
6.3.8 Summary of the Main Evaluation
Note that in some figures, certain methods may not have a monotonic probability which starts high and ends low, or vice versa. For example, in Fig. Dp. Q-LR has a u-shaped probability, suggesting that for some combinations of Q-LR with certain other G- and U-methods, its performance is good, and with others it is poor. In such cases it may be more informative consult the full results in the Appendix, to attempt to understand whether there is any particular combination dependence.
Overall, the results were quick mixed across datasets and metrics, highlighting that no single method out- performs all others across datasets of different functional form, structure, and sample size.
MSE Summary: Our Q-MN performed well on the LF datasots, particularly in smaller samples. We found that both Q-LR and Q-SL also performed consistently across the different sample sizes, even with the introduction of non-linearity with LF (v2). Indeed, with the introduction of this non-linearity, we found Q-TVAE to yield good performance, and this competitive edge held up with IHDP as well. We did not find that the choice of G-method had large impact on the results, although G-MN tended to do slightly worse. With smaller sample sizes = (500, 5000) and/or simpler datasets (LF v1), our U-multi performed the best as an update method. As sample size increased, we found that the onestep U-ones became the best performer, and similar behaviour has been found in other work NNugebauer & van der Laan 2005 For more complex datasets like IHDP, we found that U-ones and U-sub performed well.
Standard Error Summary: Once again, our Q-MN provided the tightest estimates, and did so consistently over all sample sizes and datasets except IHDP. The next best and most consistent estimator (including good performance on IHDP) in terms of the tightness of its estimates, was Q-SL. Once again, we did not find that the choice of G-method had a large impact on the results, but G-MN tended to do slightly worse than othrrs Our U-multi yielded consistently tight estimates across all datasets (including IHDP), although in general, the base models (without update steps) also performed well in this regard. As with the MSE results, U-ones and U-sub performed more competitively as the sample size increased.
Normality Summary: The choice of Q-method did not have big impact on the likelihood of normally distributed estimates for the LF datasets, although Q-D performed poorly, and the performance of Q-
21

---

Published in Transactions on Machine Learning Research (03/2023)
1.0
0.8
0.6
Q-D
O'Dnoctig
0.4
GD
Probability of normalty given methodology
0.2
Uemuli-nemmurrrr
(a)
G-R
(b)
Uffass
(c)
0.8
0.6
0.D
G-Onotreg
0.NN
GD
0-91
C.roorigg
GMN
Usuhh
0.2
Q-5
GT
(d)
(e)
(f)
1.0
0.8
100
0.4
OMN
GD
QSL
GAN
Useb
0.2 Q4R
GSL
Unewe
O-T
(g)
GF
(h)
(i)
0.0
US
GLR
500
5000
10000
500
5000
10000
500
5000
10000
Sample Size
Figure 10: Probability of p 0.01 for the Shapiro-Wilk test of normality for each Q (outcome). G (propen- sity score), and U (update step) method with the LF and general 'Gen' datasets (500, 5000, 10000). Because undertook all combinations of G and U, each point represents marginalization over the other dimension(s). For instance, for the Q methods, Q-D (DragonNet) is an average probability result when combining Q-D with all possible other G and U methods. Best viewed colour.
22

---

Published in Transactions on Machine Learning Research (03/2023)
1.0
Q-DD
G-D
0.8
Probabilty of normality given methodology
0.6
0.4
U-mutt-lin-mmanvar
Q-IR Dnotrag
U-sub
U-ones
0.2
Q-MIL
G-P G-MN
U-Base
Q-TVAE
G-B8
GDDotreg
0.0
Q-methods
G-methods
U-methods
Figure 11: Probability p 0.01 for the Chapiro-Wilk test of normality for each Q (outcome). G (propensity score), and U (update step) method with the IHDP dataset. Because we undertook all combinations G and U. cach point represents a marginalization over the other dimension(s)) For instance, for the Q mothods, Q-D (DragonNet) an average probability result when combining Q-D with all possible other G and U methods. Best viewed in color.
MN dropped as sample size increased. Surprisingly, these results reversed for the IHDP dataset, with Q-D providing the most frequently normally distributed estimates, with the other methods yielding generally poor performance. Both G-LR and G-SL worked well propensity score models for the LF-datasets, yielding a high likelihood of normally distributed estimates. However, on IHDP only the propensity score estimates from G-D were found to work well. U-ones and U-sub were found to yield consistently normally distributed errors across the LF datasets, with our U-multi unfortunately yiclding little advantage over the base model.
In some ways, the relatively disappointing results with respect to the normality of the estimates is not surprising. Benkeser et al. (0117) and van der Laan (2014) showed that the double-robustress property relating to normal limiting distribution which is afforded by estimators satisfying the efficient influence function does not apply when data-adaptive estimators are used (such as superlearners). In order for the double-robustness property to hold (with respect to the normal limiting distribution) with data-adaptive estimators, additional conditions must be satisfied. The failure to yield normally distributed estimates for many of the evaluated methods in this work thus may well be due to some degree of misspecification in the treatment or outcome models (or, indeed, both). One would expect that using the additional update steps proposed by Benkeser et al. (2017) and van der Laan (2014) would yield improved results and this presents promising direction for future evaluations and development.
6.4
Results of the Shapley Value Analysis
In addition to the presentation of the results given in Figs.
well as those given below in Figs. 20 26
we also explored whether a meta-analysis using the Shapley value approach (Shapley 1953 Lundberg et al. 2020 2017 Lundberg & Lee 2017) could provide additional insights into the performance variation across methods and datasets. Indeed, one of the limitations of the way the earlier results are presented is that they involve marginalization over one or more methodological components (e.g., obtain quantile probabilitics for the Q-methods, we have to marginalize over all G- and U-methods). This can make it difficult to between different combinations of
identify higher-order interactions performance components and datasets.
The results presented in Figs. 12 as well as those in the Appendix in Table and Figs. 17 19 represent the output from the SHapley Additive exPlanations (SHAP) machine learning explainability technique [und-] berg et al. 2020). The process as follows: we take the full set of factorial results from our experiments (all combinations from Table [) and use the choice dataset and choice/combinatte combination of Q-, G-, and U-methods predictors in regression with each of the MSE, ATE est imate standard error, and Shapiro-Wii Wilk test. p-vauue
23

---

Published in Transactions on Machine Learning Research (03/2023)
Ukoall Aevestnnn
scmple etde
GAMY
QdR
Ueeis
DMN
Q.0
ditanet
0.0
077
T
Get
GAM
OUR
Unaub
GAIN
0-1
GUR
Figure 12: Shapley predictor interaction heatmap for MSE outcome.
results as diffcrent outcomes. We use a random forest algorithm (Breiman 2001 as the regressor, with the default values the scikit-learn implementation which have been shown to yield stable and consistent per- formance (Pedregosa et al.) 2011 Probst et al. 2018). Then, the SHAP package conceives of the regression task as a game, with each predictor representing an agent which, in collaboration with other predictors, seek to maximise the performance of the regressor. The output of the Shapley analysis includes global predictor importances (which tell us, average, how useful each predictor is in explaining the regressor output), the individual impact of each datapoint in the dataset on the regressor's output, and a quantification of the degree of interaction between predictors. These results therefore help us identify whether there exist strong interactions between the choice of methods and/or datasets.
Figs. 12 and 14 show Shapley interaction heatmaps for the three outcomes MSE, ATE estimate standard error, and y- valuus respectively. Brighter values indicate the presence of interaction between predictors (where the combined information tells us something more about the likely value of the outcome than the individual predictors alone). A good method should arguably be one which does not exhibit strong interaction with the dataset and/or the sample-size. Any such interactions indicate that the method is otherwise sensitive to these factors, and therefore does not provide consistent performance.
For MSE Fig. 12 we see interactions between the sample size, G-MN, Q-LR, U-Multi and, to lesser extent, U-sub and Q-MN. Otherwise, practically no other methods show strong interactions with each other the datasets apart from the combination of G-MN with U-Sub. For the ATE estimate standard error
24

---

Published in Transactions on Machine Learning Research (03/2023)
scmple etdr
OrOnaing
Leeao
D4N
ouR
debbst
0.8
000
Q-L
077
C4NN
07144
USaee
DD
Unidd
GAMN
04R
Hoas
Unmol aoveena
G9L
0-1
0AR
6.0
Figure 13: Shapley predictor interaction heatmap for ATE estimate standard error outcome.
results in Fig. 13 very similar results, indicating strong interaction between the sample size and U-sub, G-MN and Q-LR, in particular. Otherwise, there are few notable interactions indicated by these results.
For the Shapiro-Wilk Wilk test for normality p-value results Fig. 14 size playing a role, indicating that the structural and/or /or functional form of the dataset impacts the ability of method yield normally distributed results, particularly when used in combination with Q-D. Remember here that the p-value is empirical test. for the normality of the distribution of the estimates that compute over the number of simulations, which is fixed across datasets/experiments.
we predominantly the dataset and sample
Finally, the overall variable importance summaries in Figs. 17 19 provided in Appendix I confirm the sensitivity of performance to the sample size and (for the p-value) the dataset.
7 Discussion
In this paper we have introduced some key aspects semiparametric theory and have undertaken a compre- hensive evaluation the potential semiparametric techniques to provide 'free' performance improvement for existing estimators without needing more data, and without needing to retrain them. We also proposed pseudo-ensemble NN method 'MultiNet' for simulating ensemble approach with single network,
25

---

Published in Transactions on Machine Learning Research (03/2023)
Ukoall levestann
desooet
0.7
040
Ugga
Uidvee
aLB
0g
0/786
UBBee
D4R
Umddi
GAN
GAN
Unees
OUR
DKR
Figure 14: Shapley predictor interaction heatmap for Shapiro- p-value outcome.
a new update stop variant 'MultiStep.. Our evaluation included a discussion of the choice of outcome 'Q' method, propensity score 'G' method, and the update 'U' method.
The summary of results indicates that the best results were subject to variation across datasets and sam- ple size. This particularly evident when comparing the results the LF datascts with those of the IHDP dataset, and when reviewing the Shapley interaction results. Indeed, the Shapley results quite telling, in that the most predictive variables for MSF, s.e. and normality were sample size and dataset struc- fure/functional form. This highlights a dependence of the performance on the method-dataset combination which is difficult to alleviate. That said, some methods were more stable across datasets than others. A further review of the Shapley interaction plots can be useful to researchers in selecting a group of Q- G-, and U-methods which are, least basod on these results, relatively inscnsitive to the dataset structure/ crructureffunctional form and sample size. Juggling the relativel performance metrics is difficult, however, because method might perform well in terms of MSE, and be insensitive sample size, for example, and yet fail to yield normally distributed errors. Such an estimator would be difficult to justify if one is interested in performing reliable statistical inference (such as that related to null-hypothesis significance testing).
More generally, dataset dependence was previously highlighted in the context causal effect estimation by Curth et al. (2021b), and it something which practitioners should be aware of, especially in the causal inference setting where we do not have access ground-truth, Researchers developing such methods should also, of course, be aware of this issue, because it significantly inform the evaluation design for testing and
26

---

Published in Transactions on Machine Learning Research (03/2023)
comparing different methods. These caveats notwithstanding, we found our MultiNet method to perform well as outcome method, yielding state of the art a number of evaluations, and performing particularly well on datasets with smaller sample sizes. Indeed, to an extent these results conflict with those of Farrell ct al. p021p, who showed that relatively basic neural networks were capable of excellent performance when combined with semiparamelric techniques and evaluated on their own simulated data as well as data for a direct mail marketing campaign. Arguably, their conclusions and results are less mixed than ours, although it is worth remembering that we restricted the set of methods used for the main evaluation to those with already competitive performance, and the remaining spread of our 'mixed' results may already be reassuringly tight. Unfortunately, it is difficult to say what is an acceptable level of performance, although recent large-scale work by Gordon et al. (2022) suggests that the primary challenge will be in satisfying identifiability that is, ensuring that our estimand can be expressed as a function of the observational data, that our model is structurally well-specified, and that no unobserved confounders exist which otherwise bias our estimates.
Many of the methods failed to yield normally distributed estimates. This somewhat expected given that the double robustness guarantees do not apply to the nature of the limiting distribution. [Benkeser et al. (0177) and van der Laan (2014) provide a means to augment the update step frameworks to include additional conditions which, when satisfied, extend the double robustness guarantees to the (normal) limiting distribution of the estimates. Furthermore, there are a number of conditions which must hold for the IF update techniques to work. Firstly, our estimator must be regular and asymptotically linear such that the second order remainder term (-(-) tends in probability to zero sufficiently quickly. These properties concern the sample size, the smoothness of the estimator, and the quality of the models we are using to approximate the relevant factors of the distribution. Clearly, if our initial model(s) is(are) poor/misspp moor/misppeiffied then a lincar path (or equivalently, first order VME) will not be sufficient to model the residual distance from the estimand, and the update steps may actually worsen our initial estimate. As long as our initial estimator is 'good enough' (insofar as it is regular and asymptotically linear), we can describe any residual bias using IFs. Unfortunately, we are not currently aware of a way to assess 'good onough'-ness, particularly in the causal-inference setting, where explicit supervision is not available. There may exist a way to use the magnitude of the IF to assess the validity of the assumption of asymptotic normality, and use this as a proxy for model performance, but we leave this to future work.
Many open questions remain: a similar set of experiments should be undertaken for other estimands (such the conditional ATE). Also, one may derive higher order IFs (Carone et al. 2014 van der Laan et al. 2021: van der Vaart 2014 Robins et al.) 2008) which introduce new challenges and opportunities. Additionally, it may be possible to use IFs to derive a proxy representing 'good enough'-nest í.e., whether the initial estimator close enough to the target estimand for the remaining bias to be modelled linearly. This, in turn, may also provide way to assoss the performance of causal inference methods, which would be highly advantageous given that explicit supervision will rarely be available in real-world causal inference settings. The extensions of Benkeser et al. (2017) and van der Laan| (2014) also represent an interesting avenue for further development, particularly in relation to the goal of undertaking valid statistical inference with nonparametric estimators.
7.1 Choasing Estimators
In practice, researchers are faced with a dilemma. Whilst the theory underpinning semiparametric and double-robust methods has not been called into question, we find that the real-world performance of the chosen estimators can vary across datasets, algorithmic choices (e.g., which candidate learners to included in an machine learning ensemble estimator), and sample sizes. This means that appliod rescarchers may be unable simply check the recent literature to find and use the most recent, stale-of-the-art methods, as determined by a finite set of experiments on benchmark datasets. Inevitably, there is no free lunch, and researchers are therefore forced to undertake their own evaluation process if they wish to maximise the chances that their methods will yield meaningful and robust results. The question then remains: What should be involved in this evaluation process to best expedite the selection of a method?
We discuss two scenarios: Firstly, the case when practitioners can feasibly simulate datascts with similar characteristics to the phenomenon under study. In this situation, practitioners can follow a similar process to the one documented in this work, and evaluate a set of candidate methods and select one or handful
27

---

Published in Transactions on Machine Learning Research (03/2023)
of methods which exhibit adequate exemplary performance. They can then apply these methods to their problem, and, following process of "triangulation', identify convergence in the results suggested by these chosen methods. The benefit of this approach is that (a) one can hedge over multiple (of the already better performing) methods rather than risking it all with one potentially poor performer, and (b) one can get a feeling for the consistency and sensitivity of the results to the choice of method. Results which are, for example, widely inconsistent across methods indicate that only a subset or none of the methods are working well for the specific application, and more simulation and evaluation are needed to establish a set of suitable methods.
Secondly, we consider the case when practitioners have little-to-no understanding of the underlying phe- nomenon (in particular, no understanding of the aar-linearity/complee /otlt//uunttiin form relating variables), and are therefore unable to perform meaningful simulations of the underlying phenomenon of interest. In this case, it is especially important that practitioners utilize multiple methods and, again, 'triangulate' order to understand the sensitivity of their results to the chosen methodology. Until one soes reasonable consistency across methodologies, it is necessary for researchers to consider again their choice of methods. We advise that in this second case, it is also important that researchers evaluate diverse range of meth- ods, with varying degrees of complexity. A review of the results presented in this paper does confirm the well-established notion that smaller datasets can benefit from the use of less complex methods (thus helping to prevent over-fitting). The use of cross-validation techniques can also be useful in mitigating challenges associated with fitting lighlyddata-ddaptiii learners to small datasets.
Finally, some of the methods explored in this paper did not compete well, and these can potentially be ruled out in order to reduce the choice. We say this tentatively because these methods may perform well with the right set of hyperparameters, and one cannot rule out the possibility that longer exploration of the hyperparameter space could change the order of results. The hyperparameter selection problem is yet another challenge which practitioners face. Altogether, the number of combinations of these hyperparameters, the choice of Q-, G-, and U-methods, and the choice of candidate learners included in the individual estimator ensembles reinforces the conclusion that practitioners may have to bite the bullet and undertake thorough evaluations of their methodologies, before deploying them on their real-world applications, wherever possible.
8 Broader Impact
It is always important to remember that the reliability of cansal inforence depends on strong, untestable assumptions (nol least because there is rarely any access to ground-truth in the domain of causal inference). Given the variability of the performance of the evaluated methods across datasets, in particular with regards to the normality of the estimates (and therefore also the validity of subsequent inference) any practical application of causal inference methods must be undertaken with cantion. Indeed, we recommend researchers establish the extent to which their inference depends on the methods used, by undertaking the same analysis with multiple approaches/ /estimatos...
References
A.M. Alaa and M. van der Schaar. Validating causal inference models via influence functions. ICLR, 2019.
A.M. Alaa and M. van der Schaar. Discriminative jackknife: Quantifying uncortainty in deep learning via higher-order influence functions. arxiv preprinl, arXiv:2007.13481v/ 2020.
N. S. Altman. An introduction to kernel and neighbo nonparametric regression. The American Statistician, 46(3):175-185, 1992. doi: 10. 10.1000/00001300 1999.109999992.775777
Philipp Bach, Victor Chernozhukov, Malte S. Kurz, and Martin Spindler. DoubleMI. - An object-oriented implementation of double machine learning in Python. Journal of Machine Learning Research, 2333): 1-6, 2022. URL http://jmlr. wtg/papers/vv///1-1a1.r 20wp:////2210022.mmmm.
D. Benkeser, M. Carone, M.J. van der Laan, and et al. Doubly robust nonparametric inference on the average treatment effect. Biomelrika, 104(4):863 880, 2017. doi: 10. 10.1093/biond (aamme//axx00 /asx053.
28

---

Published in Transactions on Machine Learning Research (03/2023)
R. Bhattacharya, R. Nabi, and I. Shpitser. Semiparametric inference for causal effects in graphical models with hidden variables. arXiv:2003. 12059v1, 2020.
I. Bica, A.M. Alaa, C. Lambert, and M. van der Schaar. From real-world patient data to individualized treatment effects using machine learning: Current and future methods to address underlying challenges. Clinical Pharmacology and Therupeutics, 109(1):87 100, 2020. doi: 10.1002/ept.199
P.J. Bickel, C.A.J. Klasson, Y. Ritov, and J.A. Wellner. Efficient and Adaptive Estimation for Semipara- metric Models. Spinger-Verlag, New York, 1998.
M.J. Blanca, R. Alarcon, and R. Bono. Current practices in data analysis procodures in psychology: what has changed? Pronliers in Psychology, 2018. doi: (3389//ppyg...... 2018.02555
V. Borisov, T. Leeman, K. Sebler, and J. Haug. Deep neural networks and tabular data: A survey. arXiv preprinl, ArXiv:2110.018880 2022.
L. Breiman. Random forests. Machine Learning, 45(1):5-32, 2001. doi: 10.1023/ L0023/A:10109334443222
M. Carone, I. Diaz, and M.J. van der Laan. Higher-order targeted minimum loss- based estimation. U.C. Berkeley Division of Biostatistics Working Paper Series, 2014.
H. Chen, T. Harinen, Lee J-L., M. Yung, and Z. Zhao. CausalML: Python package for causal machine learning. arXiv preprint, 2002. 1002.11634 16311, 2020.
V. Chernozhukov, D. verikor, M. Demirer, E. Dullo, C. Hansen, W. Newey, and J. Robins. Dou- ble/ /debiased machine learning for treatment and structural parameters. Econometrics Journal, 21(1): C1-C68, 2018. doi: 10.1111/cctj.120
A. Curth and M. der Schaar. Nonparametric estimation of heterogeneous treatment effects: From theory to lcarning algorithms. AISTATS, 130, 2021.
A. Curth, A.M. Alaa, and M. van der Schaar. Estimating structural target functions using machine learning and inffuence functions. arXiv preprint, 22Xiv:2008.0646111: 2021a.
Curth, D. Svensson, J. Weatherall, and M. van der Schaar. Really doing great at estimating CATE? a critical look ML benchmarking practices in treatment effect estimation. 35th Conference onf Neural Information Processing Systems (NeurIPS 2021), 2021b.
M. Dehghani, Y. Tay, A.A. Gritsenko, Z. Zhao, N. Houlsby, F. Diaz, D. Metzler, and O. Vinyals. The benchmark lottery. 2021, arXiv:2107.07002x1. arXiv preprint.
V. Dorie. Non-parametrics for causal inference. http:://gthub..oom/u/rrr////// com/udorie/npci, 2016.
B. Efron and G. Gong. A leisurely look at the boostrap, the jackkise, and cross-validation. The American Statistician, 37(1):36-45 1983. doi: Noi:10.23072665586
R.J. Evans and T.S. Richardson. Smooth, identifiable supermodels of discrete DAG models with latent variables. Bernoulli, 25(2):848-875 876, 2019. doi: 10.3150/17-BEJJ1000
M. Ezzati, A.D. Lopez, and C.J.L. Murray (eds.). Comporative Quantification of Health Risks: Global and Regional Burden of Disease Attributable to Selected Major Risk Factors, chapter Effects of multiple interventions. World Health Organization, Geneva, 2004.
M.H. Farrell, T. Liang, and S. Misra. Deep neural networks for estimation and inference. Econometrica, 89 (1):181 111,121, 2021.
A. Fisher and E.H. Kennedy. Visually communicating and tcaching intuition for inftuence functions. wwXiv:1810.03220
2019.
Y. Freund and R. Schapire. A decision-theoretic generalization of on-line learning and application to boosting. Journal of Computer and System Sciences, 19(1))119-139 1997. doi: 10.1006/jess. 1etiss.999... 1997.1504.
29

---

Published in Transactions on Machine Learning Research (03/2023)
J. Friedman. Greedy function approximation: A gradient boosting machine. The Annals of Statistics, 29(5). 2001.
Y. Gal and Z. Ghahramani. Dropout as a Bayesian approximation: representing model uncertainty in deep learning. ICML, 2016.
B.R. Gordon, R. Moakler, and F. Zettelmeyer. Close enough? large-scale exploration of non-experimental approaches to advertising measurement. arXiv preprint, 1rXiv:201.00055x2 2022.
C. Guo, Pleiss, Y. Sun, and K.Q. Weinberger. On calibration of modern neural networks. ICLR, 2017.
R. Guo, J. Li, and H. Liu. Learning individual causal effects from networked observational data. Association for Computing Machinery, 2020.
F. R. Hampel. The influence curve and its role in robust estimation. Journal of the Americon Statistical Association, (9(34)):383-333 393, 1974.
X. Han, B.C. Wallace, and Y. Tsvetkov. Explaining black box predictions and unveiling data artifacts through influence functions. arXiv preprint, arXiv:2005.06675v1, 2020.
L. Henckel, E. Perković, and M.H. Maathuis. Graphical critoria for efficient total effect estimation via adjustment in causal linear models. arXiv preprint, SEEi::1907.024355: 2020.
J. L. Hill. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 2011.
O. Hines, O. Dukes, K. Draz-Oraz, and S. Vansteelandt. Demystifying statistical learning based on efficient influence functions. arXiv preprint, arXiv:2107.00681 2021.
K. Hornik. Some results neural network approximation. Neural Networks, 6:1069-1072, 1993.
K. Hornik, M. Stinchcombe, and H. White. Multilayer feedforward networks are universal approximators. Newral Networks, 2:359-300 366, 1989. doi:
Y. Huang and M. Valtorta, Pearl's calculus of intervention is complete. Proceedings of the Secony-Secood Conference Uncerlainly in Artificial Intelligence, arXiv:1206.6831:211 224, 2006. doi: 10.5555 10.5555/3020411 0320549.. 3020446.
H. Ichimura and W. Newey. arXiv:1508.01378v2, 2021.
The influence function of semiparametric estimators.
arXio preprint,
G.W. Imbens and D.B. Rubin. Causal inference for statistics, social, and biomedical sciences. An Introduc- tion. Cambridge University Press, New York, 2015.
E. Jones, T. Oliphant, P. Petereson, and et al. SciPy: Open source scientific tools for Python. htt////wwscpppt.to org, 2001.
Y. Jung, J. Tian, and E. Bareinboim. Estimating causal effects using weighting-basee estimators. The 34th AAAI Conference on Artificial Intelligence, 2020.
A. Kadra, M. Lindaner, F. Hutter, and J. Grabocka. Regularization is all you need: simple neural nets can excel on tabular data. NewrIPS, 2021.
E.H. Kennedy. Semiparametric theory and empirical processes in causal inference. arXin:1510.0474000, 2016.
E.H. Kennedy. Optimal doubly robust. estimation of hetereogeneous causal effects. arXiv preprinl, a7X1v:2004.14497V,, 2020.
D. P. Kingma and J. L. Ba. Adam: method for stochastic optimization. arXiu:1412.69800 2017.
30

---

Published in Transactions on Machine Learning Research (03/2023)
P.W. Koh and P. Liang. Understanding black-box predictions via. influence curves. PMLR, 2017.
N. Kreif and K. DiazOrdaz. Machine learning in policy evaluation: now tools for causal inference. 1rXiv 070::1900.00/00011, 2019.
S. R. Kunzel, J.S. Sekhon, P.J. Bickel, and B. Yu. Meta-learners for estimating heterogeneous treatment effects using machine learning. aixi preprint, #XXv:1706.0046116 2019.
C.F. Kurz. Augmented inverse probability weighting and the double robustness property. Medical Decision Making, 2021. doi: 10.17770272980X21100777
J. Levy. Tutorial: Deriving the efficient influence curve for large models. 1rXii:1903.0170000 2019.
Li, S. Rosete, J. Coyle, R.V. Phillips, N.S. Hejazi, Malenica, B.F. Arnold, J. Benjamin-Chumn A. Mertens, J.M. Colford, M.J. van der Laan, and A.E. Hubbard. Evaluating the robustness of tar- geted maximum likelihood estimators via realistic simulations in nutrition intervention trials. Stalislies in Medicine, (2), 2022. doi: https://doi.org. 10. 10.1000/s/i.m344 .348.
C. Louizos, U. Shalit, J. Mooij, D. Sontag, R. Zemel, and M. Welling. Causal effect inference with deep latent-variable models. 31st Conference on Neural Information Processing Systems, 2017.
S.M. Lundberg and S-I. Lec. A unified approach to interpreting model prodictions. 31st Conference Newrul Information Processing Systems, 2017.
S.M. Lundberg, G.G. Erion, and S-I. Lee. Consistent individualized feature attribution for tree ensembles. Proceedings of the 34th International Conference on Machine Learning, Sydney, Australia, 2017.
S.M. Lundberg, G. Erion, H. Chen, A. DeGrave, J.M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal, and S-I. Lee. From local explanations to global understanding with explainable AI for trees. Nature Machine Intelligence, 2:56 67, 2020. doi: (00033//42256-09-00000
M.A. Fernadzer,,ndez, M. Schomaker, B. Rachot, and M.E. Schnitzer. Targeted maximum likelihood estimation for binary treatment: A lutorial. Statislics in Medicine, 37(16):2530 2546, 2018. doi: 10.1002/sim.7628.
R. Neugebauer and M.J. van der Laan. Why prefer double robust estimates? illustration with causal point treatment studies. Journal of Statistical Planning and Inference, 129(1):405 426, 2005.
W. Newey, Semi-parametric efficicency bounds. Journal of Applied Econometrics, 5:99 135, 1990.
W. Newey. The asymptotic variance semi-parametric estimators. Econometrika, 62:1349 82, 1994.
J. Pearl. Causality. Cambridge University Press, Cambridge, 2009.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, and B. et al. Thirion. Scikit-learn: Machine learning in Python. JMLR, 12:2825-2889 2011.
M. Petersen, Balzer, D. Kwarsima,, N. Sang, G. Chamie, J. Ayieko, J. Kabami, A. Owaraganise, T. Liegler, F. Mwangwa, and K. Kadede. Association of implementation of universal testing and treatment inter- vention with HIV diagnosis, receipt of antiretroviral therapy, and viral suppression in East Africa. Journal doi: 10.1001/jama.2017.377
of American Medical Association, 317(21): 2212122220000 2017.
K.E. Porter, S. Gruber, M.J. van der Laan, and J.S. Sekhon. The relative performance of targeted maximum likelihood estimators. International Journal of Biostatistics, 7:1034, 2011.
P. Probst, M.N. Wright, and A.-L. Boulesteix. Hyperparameters and tuning strategics for random forest. Wires Data Mining and Knowledge Discovery, 2018. doi: 20.1002/widm.133
T.S. Richardson and P. Spirtes. Causal inference via ancestral graph models. In P. Green, N. Hjort, and
S. Richardson (eds.), Highly Stractured Stochastic Systems. Oxford University Press, Oxford, 2003.
31

---

Published in Transactions on Machine Learning Research (03/2023)
T.S. Richardson, R.J. Evans, J.M. Robins, and I. Shpitser. Nested Markov properties for Acyelic Directed Mixed Graphs. arxiv preprini, arXiv:1701 7077000006662 2017.
J. Robins, A new approach causal inference in mortality studies with a sustained exposure period application to control of the healthy worker survivor effect. Mathematical Modelling, 7:1393-155 1986. doi: 00016/0000025588000
J.M. Robins, Li, E.J. Tchetgen, and A.W. van der Vaart. Higher order influence functions and minimax estimation of nonlinear functionals. Probability and Statistics: Essays in Honor of David A. Freedman, pp. 335-421, 2008.
A. Rotnitzky and E. Smucler. Efficiett adjustment sels for population average treatment effect estimation in non-paramttric causal graphical models. JMLR, 21(88)), 2020.
D. B. Rubin. Causal inference using potential outcomes: Design, modeling, decisions. Journal of the American Statistical Association, 100(469):322-33,, 2005. doi: 10.1198 041198/0162145040000001111
N. Sani, J. Lee, R. Nabi, and I. Shpitser. A semiparametzic approach to interpretable machine learning. arXiv preprint, arXiv:2006.04F5 04732 Search... arXiv:2006.0475 047722220 Search... arXiv:2006.04732, 2020.
C. Shalit, F. D. Johansson, and D. Sontag. Estimating individual treatment effect: generalization bounds and algorithms. ICML, 2017.
S.S. Shapiro and M.B. Wilk. An analysis of variance test for normality (complete samples). Biometrika, 52 (3-4):591-611, 1965. doi: 10.1093, 00.109//boome/55..302 b52.me1/52.3-4.591.
L.S. Shapley, A value for n-person games. Contributions to the Theory of Games, 2(28):307 317, 1953.
C. Shi, D. M. Blei, and V. Voitch. Adapting neural networks for the estimation of treatment effects. 33rd Conference on Neural Information Processing Systems, 2019.
C. Chi, T. Xu, and W. Bergsma. Double generative adversarial networks for conditional independence testing, arXiv:2006. 0tle:20006.0215v1, 2020.
Shpitser and J. Pearl. Identification joint interventional distributions in recursive semi- Mami-Maroovian causal models. Proceedings of the National Conference on Artificial Intelligence, 21:1219-1226, 2006.
R. Shwartz Ziv and A. Armon. Tabular data: Deep learning is not all you need. Information Fusion, 81: 84-90, 2021. doi: 10.1016/j.infful 1.000///.nffus2222111...1111 2021.11.011.
B. Siegerink, W. den Hollander, M. Zeegers, and R. Middelburg. Causal inference in law: an epidemiological perspective. European Journal of Risk Regulation, 7(1):175 186, 2016. doi: 101007/51187999X00007777
C. Szegedy, W. Liu, Y. Jía, P. Sermanet, S. Reed, D. Anguclov, D. Erhan, V. Vanhoucke, and A. Rabinovich. Going deeper with convolutions. CVPR, 2015.
J. Tian and Pearl. A general identification condition for causal effects. AAAI, 2002.
A. Tsiatis. Semiparametric Theory and Missing Data. Spinger, New York, 2006.
M. J. van der Laan and Gruber. Targeted minimum loss based estimation of causal effects of multiple time point interventions. Int. J. Biostat, Art 9(11,, 2012.
M. J. van der Laan and S. Rose. Targeted Learning Causal Inference for Observational and Experimental Data. Springer International, New York, 2011.
M. J. van der Laan and R. J. C. M. Starmans. Entering the era of data science: targeted learning and the integration of statistics and computational data analysis. Advances in Statistics, 2014.
M. J. van der Laan, Z. Wang, and van der Laan. Higher order targeted maximum likelihood estimation. arXiv:21010022000 2021.
32

---

Published in Transactions on Machine Learning Research (03/2023)
M.J. van der Laan. Targeted estimation of muisance parameters to obtain valid statistical inference. Inter- national Journal on Biostalisties, 10:29-57, 57, 2014.
M.J. van der Laan and D.B. Rubin. Targeted maximum likelihood learning. The International Journal of Biostatistics, 2(1); 2006. doi: 00.2222/155--469110453
M.J. van der Laan, F.C. Polley, and A.F. Hubbard. Super Learner. Statistical Applications of Genelios and Molecular Biology, 6(25), 2007. doi: 10.2202/1544-000 1544-6115.1100.
A.W. van der Vaart. Higher order tangent spaces and influence functions. Statistical Science, 29(4):679-66) 686 2014.
T. Verma and J. Pearl. Equivalence and synthesis of causal models. Proc. 6th Conf. on Uncertainty in Artificial Intelligence, 1990.
M. J. Vowels. Misspecification and unreliable interpretations in psychology and social science. Psychological Methods, 2021. doi: 10.1037 10.1037/met0000 0000424..
M. J. Vowels, N.C. Camgoz, and R. Bowden. Targeted VAE: Structured inference and targeted learning for causal parameter estimation. IEEE SMDS, 2021.
Y. Wen, D. Tran, and J. Ba. BatchEnsemble: An alternative approach to efficient ensemble and lifelong learning. ICLR, 2020.
D.H. Wolpert and W.G. Macready. No free lunch theorems for optimization. IFFE Transacions on Evolu- tionary Computation, 1(67), 1997. doi: 10.1109/4435.555888
P.A. Wu and K. Fukumizu, Causal mosaic: cause-effect inference via nonlinear ICA and ensemble method. AISTATS, 108, 2020.
P.A. Wu and K. Fukumizu. Intact-VALL Estimating treatment effects under unobserved confounding, ICLR, 2022.
L. Yao, S. Li, Y. Li, M. Huai, J. Gao, and A. Zhang. Representation learning for treatment effect estimation from observational data. 32nd Conference on Neural Information Processing Systems (NeurIPS), 2018.
L. Yao, Z. Chu, S. Li, Y. Li, J. Gao, and A. Zhang. A survey on causal inference. ACM Transactions on Knowledge Discovery from Data, 15(5):1 46, 2020. doi: 10.1145/3144944.
J. Yoon, J. Jordan, and M. van der Schaar. GANITE: Estimation of individualized treatment effects using generative adversarial nets. ICLR, 2018.
Zhong and J-L Wang. Neural networks for partially linear quantile regression. arXiv preprind, MrXIv:2106006225, 2021.
Zou and T. Hastie. Regularization and variable selection via the elastic net. R. Statist, Soc., 67()). 301-320, 2005.
33

---

Published in Transactions on Machine Learning Research (03/2023)
A Causal Assumptions
The causal quantity can be estimated in terms of observational (and therefore statistical) quantities if number of strong (but common: Yao et al. 2020: Guo et al. 2020 Rubin 2005 Imbens & Rubin 2015;
Vowels et al. 2021) assumptions hold: (1) Stable Unit Treatment Value Assumption (SUTVA):: the potential outcomes for each individual or data unit are independent the treatments assigned to all other individuals.
(2) Positivity: the assigument of treatment probabilities are non-zero and non-delerministic P(T 4 x,) 0, x, (3) Ignorabilityy Unconfoundedness/ Conditional Exchangeability: There are no unobserved confounders, such that the likelihoods of treatment for two individuals with the covariates equal, and the potential outcomes for two individuals with the same latent covariates are also equal s.t. T IL
B Statistical Inference with Influence Functions
Machine learning techniques, while parameterized per se, often non-parametric insofar the estimates they yield are not directly parameterizable (e.g.) a Gaussian with mean and variance. However, such parameterization is extremely helpful in facilitating statistical inference, including the ubiquitous mull hypothesis significance test, which is straightforward when the parameter being tested is normally distributed. Whilst other approaches to statistical inference with non-parametric methods exist, such as the bootstrap approach [Efron & Gong 1983), influence functions have been recognized as valid approach for some time
Bickel et al. 1998 Tsiatis 2006).
Following van der Laan & Rose 2011 7.75) can derive 95% confidence intervals from the influence
function to be (assuming normal distribution):
1 Σ
Var(o) (0) = Σ olz,)
12 {Zz,)
Var(ø)
50
(8)
V (P..)) + 1.96se.
V (Pn)
Peal = 2
1 - +
where Te (Px) is the estimated target quantity after bias correction has been applied, is the CDF of a normal distribution, is the standard error, and Poat is the p-value.
C IFs for General Graphical Models
In this paper, we focus on the estimation of average treatment. effect in the setting of Fig D.. However, the methods discussed in this paper can applied for more complex estimands with an arbitrary causal graph structure, as long as the estimand at hand is causally identifiable from the observed data. In this scction, we discuss the derivation of IF's for general form of an estimand in general graphical model.
C.1 Influence Function of an Interventional Distribution
The causal identification of interventional distributions is well-studied in the literature. In the case of full obsorvability, any interventional distribution is identifiable using (extended) g-formula (Ezzati al.] 2004 Robins 1986). some variables of the causal system are unobscrved, all interventional distributions are not necessarily identifiable. Tian & Pearl (2002) and Shhiitser & Pearl (2006) provided necessary and sufficient
34

---

Published in Transactions on Machine Learning Research (03/2023)
conditions of identifiability in such models. The causal identification problem in DAGs with unobserved (latent) variables can equivalently be defined on acyclic. directed mized graphs (ADMGs) (Richardson & Spirtes 2003 Richardson et al. 2017 Evans & Richardson 2019). ADMGs are acyclic mixed graphs with directed and bidirected edges, that result from a DAG through a latent projection operation onto a graph over the observable variables (Verma & Pearl 1990).
Pearl's do-calculus shown to be complete for the identification of interventional distributions Huang & Valtorta 2006). Let V denote the set of all observed variables. Starting with an identifiable interventional distribution P(y [T t')), an identification functional of the following form is derived using do-calculus:
P(y]do(T= = =t)) 1l, b
(9)
(
S
where a, b Cj, and dj are realizations of A₁, Bj. Cj, and Dj, respectively, and A,, B,, Cj. Dj. S are subsets of variables such that for each i and j. An and nj Ø. Note that the sets B; and Dj might empty. The I symbol in Eq. indicates a sunmation over the values of the set of variables $ in the discrele case, and an integration over these values the continuous setting. To derive the influence function of Eq. g we begin with conditional distribution of the form P(a[b). o Ø, we can write
Pv(v) = (1 - -e)p(v) + t,v()),
Pc(a(b) P.(a,,b b)
Pc(b)
dab(a
82,6(a,b) b) -P(a.b)
P(a, b)(,,b)) - P(b)]
(10)
de
P(b)
P²(b)
k==
(00.6(a,))
6g(b)
= P(a(b)
P(a,b) P(b)
where V is the point that we compute the influence function at, and a, b are the values of sets of variables A, B C V that are consistent with v. For an empty b, using similar arguments, we have:
(P)
(a(a) - P(a) 1).
= P(a)
(11)
de
e=0
With slight abuse of notation, for b 2, we define the IF of Eq. g
b.(bb P(b)
1. Using Eq. 10 and Eq. 11
we can now derive
o(v.P) d((1 c)P + cdo)
de
l==
(12)
(ay X
6 b, (an, bi) P(ai, bi)
65, (b;) P(bi)
Se (c,, (c) dj) P(cj,dj) dj) P(d,)
(d,(d)
Σ
n,P((ydd))
(.(cd) Pf(cld) Equation 12 is the foundation to the approach that shall be
Note that we used de P.(e discussed in the following section for deriving the IF of general class of estimands.
C.2 Influence Function of General Estimand
We have so far discussed the influence function of a causal effect of the form P(y|do(T t')). In this section, show how IF's can be derived for any general estimand of the form:
V Ep(x(P)),
(13)
35

---

Published in Transactions on Machine Learning Research (03/2023)
where x(.) is functional. Then we have:
P. = P. (1 ()-bPp
(P.))
dv (P.)
/
dPe
dk
dPe
((P.)
Pe
dv
de
dc
dP.
dc
(14)
1=0
/ (M(P) +
dk
dP:
p
dv
dp
de
| x(P) de
dPc
dk
dP.
dv +
Ep
dP
de
l==0
The value of de c= can be plugged into Eq. 14 which completes the derivation of the IF for the estimand in Eq. 13 As an example, if the queried estimand is the average density of variable Y, that is, K is the identity functional, then:
/ pd (y))dy..
((P)
(P.)
(P+1
dP.
dy
de
de
c=0
c=0
dP:
2P(y)
dy-
de
e==0
Algorithm summarises the steps of our proposed automated approach to derive the influence function of an estimand of the form presented in Eq. 13 given a general graphical model. Note that if the effect is identifiable, this algorithm outputs the analytic influence function, and otherwise, throws failure. A demonstrative example can be found in the associated code repository in the form of a notebook, and/or / in the attached supplementary code.
Algorithm IF of identifiable effect.
input: An estimand (TP)) of the form of Eq. 13
interventional distribution P, causal graph g
output: The analytic IF of (P()) if P is identifiable, fail o.w.
1: if P is identifiable then
2: P + the identification functional of (Eq. 9 using do-calculus
3:
Φ t the IF P as in Eq. 12
4:
(P()
+ the formulation as in Eq. 14
de
4 - P-ug @ into dV(P.)
5:
de
return
clsc
return FAIL
D Influence Function Update Methods
D.1 One-Step and Submodel Approach
Using the one-slep approach, the original estimator (P(P.) can be improved by a forward application of the Von Mises Expansion (VME) of Eq. 3 one takes the initial estimator and adds to it the estimate of the IF to yield an updated estimator which accounts for the 'plug-in bias'. In the case of the ATE, this
yields the augmented inverse propensity weighted (AIPW) estimator (Hines et al. 2021 Neugebauer & van
der Laan 2005
Kurz 2021
36

---

Published in Transactions on Machine Learning Research (03/2023)
The second submodel approach updates the initial estimate by solving Σi=1 done by estimating the degree to which the principal estimator Q is being biased, and correcting for it. This approach works by first constructing a parametric submodel in terms of the plug in estimator Q(t, X) and a function H of the propensity score G (which represents the biasing quantity), and with these can be used to derive an updated plug-in estimator (", (I. In the expressions which follow, we assumee binary treatment, T and have replaced the Dirac delta functions with indicator functions. First, the updated estimator Q*(1, x) can be expressed in terms of biasing quantity H as follows:
=0, 0, which can be
Q' = 1,x,) Q(T = 1,x,) + H T = 1), 1, (1)
where Hz, T 1) G(x)
The equivalent for when treatment is 0 can be expressed as: -T= Q(T = 0,,, + (z:, 0), where Hz,,T 0) 1 1, (0)
(15)
1 - G(x)
H(z, 4;) is known as the clever covariate. The parameter f is estimated as the coefficient in the associated intercept-free free 'maximum- -enximmm-iikelihoodd linear regression' reflected in the first line of Eq. 15 above. When the up- dated Q" is substituted into Eq one solves in one update what is known as the efficient influence function, and following the update, the mean of the IF, as well as the residual bias will be zero. In practice the two methods yield different results with finite samples (Porter et al.) 2011 Benkeser et al.] 2017). In particular, the one-step / AIPW estimator may yield estimates outside of the range of values allowed according to the parameter space, and be more sensitive to ncar-positivity violations (i.e., when the probability of treatment is close to zero) owing to the first term on the RHS of Eq. 4 "Luque-Fernande et al.] 2018). In contrast, the submodel approach will not, because it is constrained by the intercept- frree regression step for instance, if the outcome Y is binary, the logistic link function will prevent any 'out-of-bounds' convergence behavior which might otherwise result in unstable estimate for the coefficient 4. Thus, both the one-step and submodel approach achieve the same aim: solving for the efficient influence function. However, the one-step approach does so naively according to the VME, and the second does so according to an additional regression stage (but nonotheless still only requires a single update to the original outcome model Q)
D.2 Targeted Regularization
Targeted regularization has, to the best of our knowledge, only been used twice in the NN literature, once in DragonNet (Shi et al. 2019), and once in TVAE (Vowels et al. 2021,. both of which were applied to the task of causal inference. The idea is to solve the efficient influence curve during NN training, similarly to Eq. [5] on a per--aaae batch basis. The parameter & in Eq. 15 is treated as a learnable parameter, trained as part of the optimization of the NN. The submodel update in Eq. 15 is thereby recast as a regularizer which influences the weights and biases of the outcome model Q(t, x). In total, then, the training objective is to minimize the sum of the negative ligelioel (NLL) of the outcome model Q(t,x) which has parameters 0 (which comprises NN weights and biases), and the NLL of the updated outcome model (*(,, x)), which is parameterized by both I and 9. As the second NLL term involves the clever covariate H, which in turn involves the plug-in estimator for the propensity score G(Z), we also need a model for the treatment which may be trained via another NLL objective, or integrated into the same NN as the one for the outcome model. Due to the paucity of thcoretical analysis for NNs, is not clear whether targeted regularization provides similar guarantoes (debiasing, double-robustness asymptotic normality) to the one-step and submodel approaches, and this is something we explore empirically.
37

---

Published in Transactions on Machine Learning Research (03/2023)
E LF (v1 and v2) and Gen Dataset Details
The generating equations for [F(11) are:
Xi Be(0.5), X2 Be(0.65),
X3 in,[U(),4)], X4 ~ ~
(0,
Be(pr), where
(16)
Pr 0.05, +
Y1 0.5 -
+0.15x2x.1,, -0.15X2x4),
where [nt is an operator which rounds the sample to the nearest integer, Be Bernoulli distribution, U is a uniform distribution, is the sigmoid function, and Yi and Yo are the counterfactual outcomes when T and T 0, respectively. Covariate X1 represents biological sex, X2 represents age category, X3 represents cancer stage, and X4 represents comorbidities.
We create variant (v2) of this DGP by introducing non-linearity into the outcome, and then into the treatment assignment as follows:
Yi (erpp|111 0.1X,, 0.35X, + 0.25X, + 0.2X, + 0.15X2x.)).
(17)
Figs.
provide information on the propensity scores for the vl and v2 variants (the second version has the same propensity score generating model as v1).
15 and 16
4.0
LF (v1. v2)
3.0
Density
2.5
1.5
1.0
0.5
0.0
0.0
0.4
0.6
0.8
1.0
Probability
Figure 15: Marginal propensity scores for the LF (v1) and LF (v2) datasets. Note that the minimum probability of treatment in random draw from the DGP is 0.007. The datasets are intentionally designed such that certain subgroups are unlikely to reccive treatment, resulting in near-positivity violations.
38

---

Published in Transactions on Machine Learning Research (03/2023)
1=
Density
2
0.0
0.6
0.8
1.0
Probability
Figure 16: Propensity scores by treatment assignment for sample from the LF (v1) dataset.
Finally, the generating equations for the generalized version of the LF dataset, which we refer to as 'Gen', are:
Xi Be(0.5), X2 Be(0.55),
X3 int(U(O, (0, X4
Ux. N),, X; = 0.2%1 + Ux.:
h ~ 1)0,1),
1 ~ 1).
R1 ~ 1)0,1),
R2 ~ N(011)
UT NVO01), 1), T~ Be(pr). where
(18)
Pr of o--5 0.05X, 0.25X; 0.6X4 0.4X2X4 +0.1XX, 0.17 0.157 +.11
UM N)0,1), Mi -0.8 0.8+0.155
Mo m0.55UM
Y1 = (erpl-11 M1 -
0.25X, 0.2X, + 0.15,,,, R1 R2l),
Y, = o 11 Mo - 0.1X, - 0.35X, + 0.25X, + 0.2X, 0.15,,X, + Rı R2),
Y = Y + Yo(1-X)
~ 1), C= + 0.4X 0.4Ucc
where int[.] is operator which rounds the sample to the nearest integer, Be is Bernoulli distribution, U is a uniform distribution, N is a normal distribution, a is the sigmoid function, I and l2 are instrumental variables, R₁ and Rg are risk variables, Mj and My are the counterfactual mediators and Y1 and Yo are the counterfactual outcomes when T and O, respectively. C collider, and finally, X1-Xş are confounders.
F Hyperparameters
The LR and SL approaches implemented using the default algorithns in the scikit-learn package (Pe- dregosa et al. 2011), whilst the the DragonNet, S-learner, T-learner, and P-learner, the CausalML. package (Chen et al. 2020). For DragonNet the number of neurons per layer was set to 200, the learning rate set to x 10-1, number of epochs = 30, and batch size = 64. For TVAE the dimensionality of all latent variables was set to 5, the number of layers set to 2, batch size = 200, number of epochs = 100, learning rate = 5 x 10-4, and targeted regularization weight of 0.1.
implemented using
For CFR and MN, we undertake a Monte-Carlo train-test split hyperparameter search with 15 trials, for every one of the 100 samples from the DGP. The best performing set of hyperparameters is then used to
39

---

Published in Transactions on Machine Learning Research (03/2023)
Table 2: Hyperparameter search space for CFR and MN based methods.
Parameter
Min
Max
Batch size
10
64
1.2 Weight Penalty
1e-5
1e-3
No. Iterations
2000
10000
Learning Rate
1e-5
1e-2
No. Layers
2
14
Dropout Prob.
0.1
0.5
No. Neurons per Layer
5
200
train CFR and MN on the full dataset. For the hyperparameter search itself, we undertake 15 trials on a arain/test split for cach of the 100 samples from the DGP, and additional, separate hyperparamcter scarches are undertaken for methods using largeted regularization. The hyperparamelers explored for CFR and MN based methods are shown in Table 2 Note that the iteration count is not in terms of epochs it represents the number of batches sampled randomly from the dataset. The number of iterations can be multiplied by the batch size and divided by the dataset size to approximately determine the equivalent number of epochs this represents.
Note that, unlike in traditional supervised learning tasks, using the full data with causal inference is possible because the target estimand is not the same quantity the quantity used to fit the algorithms (Farrell et al., 2.21 Indeed, whilst cross-fitting is used for the hyperparameter search, subsequent use of the full data has been shown to be beneficial, especially in small samples (Curth & van der Schaar 2021). It is reassuring to note that overfitting is likely to worsen our recorded estimates, rather than misleadingly improve them. Similarly, even though the SL is trained and the corresponding weights derived using a hold-out set, the final algorithm is trained on the full dataset for estimation. Logistic regression is simply trained on the full dataset without any data splitting. This is what motivated us to ask the question as to whether or not it is possible to have 'free lunch' with IFs. Indeed, if no additional data is required, but improve our estimates and achieve valid statistical inference (for the purposes, for example, of null hypothesis significance tests), then this would represent a valuable gain. For all treatment models, we bound predictions to fall in the range [9755 (Li et al. 2022).
nonetheless
G Things that Did Not Work
G.1 Calibration
One of the initial possibilities that we considered which might explain why some methods (e.g., CFR) were not performing as well as others, was that the calibration of the output might be poor (Guo et al. 2017). However, we tried calibrating the trained outcome and treatment model networks using temperature scaling, We found to be unsuccessful, and we leave an exploration of why it failed to future work.
G.2 Restricted Hyperparameter Search
Additionally, we tried only performing hyperparameter search with held-out test set once at the beginning of the 100 subsequent simulations for each model and dataset variant, rather than performing it for every single simulation. This did not work, and we found that if the first network 'designed' through hyperparameter scarch happened to be degenerate with respect to its performance as plug-in estimator (notwithstanding its potentially adequate performance as an outcome model), then it will be degenerate for all simulations, and yield incredibly biased results. However, performing hyperparameter search for every simulation more accurately represents the use of these algorithms in practice.
This problem also highlights the importance of fitting multiple neural networks on the same data. As supervision is not available, the usual metrics for hyperparameter search (based on e.g., held out data loss scores) can be poor indicator for the efficacy of the network as a plug-in estimator. By re-performing hyperparameter search, even on the same data (put perhaps, with different splits), one can effectively bootstrap to average out the variability associated with the hyperparameter search itself. Indeed, as the
40

---

Published in Transactions on Machine Learning Research (03/2023)
results show, the average estimates for the ATE using CFR net are close to the true ATE, even if the variance of the estimation is relatively high. We leave a comparison of the contribution of variance from hyperparameter search to further work.
G.3 MultiStep Update Variants
Relating to our proposed MultiStep objective, we also tried a non-linear, generalized variant with an objective which still attempls to minimize the mean and variance of the influence function but such that the update step is parameterized as follows:
go(ui Q(1, x:), (2H)
(19)
It can be seen that instead of optimizing over the domain of g E in Eq. where 0 are the parameters of a shallow NN function g- Here, v1 E {0, 1] and u2 E {0. 1) are hyperparameters determining whether the NN function ge should be taken over just the clever covariate H, over both the clever covariate and the outcome model Q.
0 we instead optimize over A O
In practice however, this approach did not yield good estimates. Furthermore, we found that MultiStep update step given in Eq. 5 setting 01 = 0 (?.e., no mean-zero penalty) also did not work well. This result was surprising because similar approach in NNugobauer & van der Laan (2005), which did not include a mean-zero penalty, yielded an improvement. However, it is also intuitive that the two properties of the Efficient Influence Function are (1) mean-zero and (2) minimum variance, then it makes sense that optimization objective should benefit from the inclusion of both of these conditions.
H Initial Evaluation
Table 3: Initial results over a restricted set of model variations. All update steps use the same propensity score G- algorithm as their Q-model algorithm, unless indicated by 'w/ G-SL', which indicates the use of SuperLearner. Mean Squared Errors (MSE) and standard error (s.e.) (lower is better) and Shapiro-Wai Wilk test p-values for normality (higher is better) for 100 simulations. Best results are those competing across all three dimensions. Bold indicates the best result for each algorithm. Multiple methods may perform approximately equally well.
Dataset Q Model
U-Base MSC
S.c.
U-oras MSE
5.0.
Ulsub MSE
S.c.
Treg MBE S.D.
Treg U-sub MSE
S.D.
U-ones w/ MSE
G-SL
U-sab w/ G-SL
MSE
s.c.
IR
.01
.0031
2012
.270
.0007
.003
.248
.0008
.008
.378
.0006
.008
.591
.0008
.003
.0008
.003
.651
.0009
.008
CFR
003
T.F (+1)
MI--ne
.003
391
.001
.729
.012
.003
.81
.001
.0008
.003
001
OR3
MN-Loe+LMM
.355
8009
.015
.141
.0007
.003
.578
.009
.OUS
10
.0117
.001
.967
.0011
.003
.669
.0008
.008
776
1000
JOLS
MIN-Cisc
4818
.014
.251
J8014
.012
11018
.013
JOS3
J0086
.007
.702
.0645
.034
.831
.0007
.003
3330
5000
MIN-Caec+TAM TM
K65
0158
006
0018
.072
K05
231
.0037
.0381
003
2095
013
.747
.0007
.008
4755
001
013
0003
.497
.003
SL. CFR
349
0017
.938
0008
.0033
.02
.0000
.003
0155
.01
.006
.015
12
.0151
.012
.315
.01
.008
JO192
.032
.023
.0007
.003
.005
0011
J003
LP (x2)
MIN Inc
.119
.001
.015
204
0006
.003
211
1008
.015
0002
(0009
0002
.029
.0008
.003
.358
5000
010
.049
.000S
.003
MN- ILM
0011
.438
00000
.0033
813
10011
.006
00025
.0005
.0009
.003
MIN-Clse
Cusc
.349
JX32
.003
.031
.0U7
.0112
.0066
(0133
010
LK
022
1818
.013
1)
J156
205
J8611
2014
.1022
0115
.0
.0597
.08
SL CIR
0856
.002
.0
.0311
.033
J0
.1134
.004
77h8
/08
.2455
.774
.0
.0450
.052
25.5
.3
-0
.0604
.051
.3125
003
1.7
.114
IHDP
MI-Lue
A1524
.044
.0243
A4425
.042
L.M
.043
.3849
.038
-0
.0243
.022
.042
MN. Case
.1977
J045
.0737
.04
10
.004
.04
2.9
.115
J0
.102
.042
.0816
.042
.0383
.047
n- 747
MI-CCoc+LM
4.7
.108
1.4
.085
.2118
.013
23.9
.164
.0
.1824
.06
-0
1.1
.078
4.7
.002
H.1 Initial Evaluation
We share initial results in Table 3 These results were used to inform a subsequent set of experiments with a restricted set of variants. Specifically, we used these to select the most successful variant of MultiNet.
For LF (v1), we see that the base CFR performs significantly worse in all considered metrics than LR and SL. Base LR and base ST. achieved the best results in terms of MSE and s.e., although note that none of
41

---

Published in Transactions on Machine Learning Research (03/2023)
the base algorithms achieve asymptotic normality. Notice that LR's base MSE performance LF (v1) is actually better than its MSF performance using the one-step and submodel updates. Such behaviour has been noted before by Luque-Fernandez et al. (2018), and occurs when the base learner is already close and/or when both outcome and treatment models are misspecified. Unlike CFR, our MN-Inc and MN-Casc variants worked well as either outcome or mreatment models, yielding the best results with the one-slep update. The other two of our MN- variants also performed well with the one-step and submodel updates but required a SL treatment model to do so.
The potential improvements for LR in combination with update steps is more striking for LF (v2). Here, the LR base outcome model is misspecified (LF v2 has an exponential outcome model). Combining the LR with the SL one-step and submodel update processes enabled the LR method to perform well in spite of the non- non-linearity of the outcome. This is a demonstration of double-robustun is misspecified, the treatment model not (or at least, is sufficiently correctly specified), owing to the use of a SL, and the estimates are improved. As with the LF (v1) dataset, combining CFR with IFs resulted in substantial improvement, especially when using an SL treatment model, yielding competitive MSE, s.e., and normally distributed estimates (thus amenable to statistical inference). These results demonstrate the power of semiparametric methods for improving our estimation with NNs, and again illustrate the double- robustness property: the CFR outcome model was poorly specified, but was able to recover with an SL treatment model. Similar performance for our MN- variants on LF (v1) was observed with LF (v2).
- even though the outcome model
Unfortunately, no method variant yielded normally distributed estimates with the IHDP dataset. The worst performing estimator across any combination of semiparametric techniques was LR. This makes sense given the non-linearity in the THDP outcome process (Curth et al., (021b). The SL with the one-step submodel updates performed equally (poorly) as the best CFR and MN- CNsC variants, although the SL provided a smaller s.e.. Overall, the best methods were our MN-Inc and MN-Inc+LM variants in combination with either a one-step update, or a one-step update using SL treatment model.
The MultiNet variant which performed the best and most consistently across all datasets was our MN-Ine (or equally, MN-Ine+LMM with the one-step update. Whereas other methods benefited from the help of SL treatment model, MN-Ine worked well as both an outcome and a treatment model, making it the best all-rounder across datasets, as well as the least dependent on the SL for correction. For all NN based approaches, largeted regularization made little difference, and sometimes resulted in instability and high MSEs. Further work is required to investigate this, although it may relate to which treatment model is used, and the associated sensitivity to positivity violations. A prior application also described the potential for the regularization to be inconsistent (Shi et al. 2019
For all base learners, we observe the potential for improvement using the semiparameteric techniques, pri- marily for improving the associated MSE. It is also worth noting that in general, the base CFR method has consistently higher (i.e., worse) s.e. than the MN-variants, although combining CFR with an udpdate step (e.g., one-step w/ ST.) significantly tightened the s.e..
Additional Results AAaalssss
In Table we also provide a set of results (without any marginalization over any of the G-, U-, or Q- dimensions) for a subset of the methods considered in the full-factorial design. Note that whilst we have tried to highlight the best results in bold and underline, many of the results are close/ compeetioeee and illustrate (again) that the performance is datasct and combination dependent, as well as that there exist multiple possible 'best' options for a given situation.
In Table 5 we provide results for the absolute difference between the true average treatment effect and the estimated average treatment effect (acATE) over the r 100 simulations. We define this as Σ |7 - Al. Note that the results here were conducted over different simulations to the results in the other tables (hence why we report again the p- values and standard errors).
42

---

Published in Transactions on Machine Learning Research (03/2023)
Table 4: Unmarginalized results over a restricted set of model variations including the Double Machine Learning (DML) method (Chernozhukov et al. 2018). Mcan Squared Errors (MSE) and standard error (s.c.) (lower is better) and Shapiro- Wilk test p-values for normality (higher is better) are provided ad computed over 100 simulations. Best results those competing across all three dimensions. Bold indicates the best result for each algorithm. Multiple methods may perform approximately equally well.
Dataset
Model
U-Base
G-SL U-nmlti
G-SL+U-aub
G-MN+U-multe
L-MN++U-uub
MSE
s.b.
MSF
R.C.
P
MSR
MSR
8.6.
P
MSE
Q-MN
0.985
0.0024
0.054
0.959
0.0023
0.054
0.027
0.0104
0.109
0.945
0.0024
0.044
0.125
0.0380
0.093
LF (v])
Q-SL
0.030
0.0032
0.057
0.034
0.0036
0.060
0.023
0.0081
0.102
0.043
0.0035
0.090
0.053
0.0081
0.995
G-TVA
0.617
0.0033
0.850
0.3%6
0.0042
0.052
0.074
0.0104
0.108
0.716
0.0046
0.052
0.010
0.0383
0.003
(n 500)
Q-T
0.700
0.000
0.084
0.772
0.0056
0.065
0.485
0.0084
0.050
0.843
0.0057
0.064
0.050
0.0073
0
DML
0.113
0.0055
Q-MNN
0.000
0.000
0.025
0.0
0.0005
0.024
0.599
0.0009
0.031
0.0
0.0006
0.025
0.583
0.0917
0.030
LF (v])
Q-SL
0.024
0.0004
0.022
0.055
0.0005
0,023
0.536
0.0009
0.033
0.012
0.0004
0.224
0.008
0.029
Q-TVAR
0.444
0.0007
0.024
0.182
0.0006
0.025
0.655
0.0009
0.031
0.368
0.0007
0.024
0.414
0.008
0.830
(n 5000)
Q-T
0.976
0.0010
0.032
0.989
0.0010
0.002
0.985
0.0010
0.033
0.823
0.0010
0.032
0.854
0.0010
0.033
DMG
0.717
Q-NNN
0.0
0.000
0.027
0.0
0.0007
0.026
0.366
0.0044
0.020
0.0
0.0006
0.025
0.648
0.0006
0.022
LP (v])
Q-SL
Q-TVAE
0
1.0
0.0004
0,0002
0.019
0.015
0,993
0.976
0,0002
0.0004
0,099
0.016
0.625
0,292
0.0004
0,0004
0.020
0.019
0.997
0.782
0.0004
0.0002
0.015
0.018
0.861
0.0004
0.0004
0.019
0.020
(n 10000)
Q-T
0.341
0.0004
0.020
01331
0.0004
0.021
0.326
0.0004
0.021
0.141
0.0004
0.021
0.122
0.00011
0.021
DML
0-777
QUNN
0.439
0.0018
0.043
0.149
0.0019
6.444
0.014
0.0199
0.114
B.076
0.0019
0.044
0.0
0.0102
0.108
LP (v2)
Q-SL
0.0
0,0027
0,034
0,001
0.0029
0,077
0.025
0,0097
0.103
0.0
0.0029
0.057
0.0
0,0095
0.103
T-TVAE
0.087
0.0075
0,945
0,031
0.0000
0.046
0.0225
0.114
B.01
0.0067
0.045
0.0
0.0112
0.112
(n 500)
Q-T
0.817
0.0066
0.064
0.721
0.0002
0.065
0.047
0,0094
0.102
0.753
0.0064
0.065
0.003
0.0110
0.108
DML
0.017
0.050
QUNN
0.0
0.0011
0.029
0.0
0.000
6027
0.925
0.0010
0.001
0.0
0.0011
0.027
0.662
0.0010
0.030
LP (v2)
Q-SL
Q-TVAE
0.765
0.007
0.0008
0.0007
0.028
0.028
0.689
0.795
0.0088
0.0007
0,227
0,028
0.921
0.912
0,0009
0,0009
0,030
0.031
0.799
0,661
0.0006
0.0007
0.028
0.027
0.678
0.072
0.000
0,000
0,030
(n 5000)
Q-T
0.793
0.0010
0.032
0.684
0.0000
0,002
0.724
0.0011
0.033
0.755
0.0010
0.032
0.0
0.0010
0.037
DML
0.824
0.0012
QEEN
0.0
0.0000
0.027
0.0
0.0008
6.025
0.707
0.0005
0.020
6.0
0.0011
0.025
0.783
0.0009
0.026
IP (v2)
Q-SL
Q-TVAE
0.859
0.0
0.0004
0.0177
0,412
0.019
0.448
0.629
0.0004
0.0001
0,021
0,019
0.984
0.981
0.0005
0,9004
0.020
0.020
0,240
0.584
0.0044
0.0004
0.021
0.019
0:227
0.676
0,004
0,0001
0.021
(n- 10000)
Q-T
0.0
0.0124
0.430
0.719
0.0006
0.022
0.566
0.0006
0.022
0.674
0.0005
0.022
0.285
0.0055
0.022
DML
0,500
0.005
CEN
0.673
0.0011
0.0319
0.1815
6,0086
0.0622
0.6488
0.0025
0.0129
0.000
0.003
0.074
Gen
Q-TVAE
Q-SL
0.5014
0,9082
0.0015
0.0035
0.0356
0.0438
0,8349
0.8841
0,0013
00035
0.0873
0.0001
0.0062
0.0779
0.1013
0,0077
0.4463
0.0014
0.0035
0.0968
0.0453
0.0001
0,0055
0,005
0.0725
(n 500)
Q-T
0.5179
0.0025
0.0496
0.5921
0.0024
B0002
0.00R9
0.00G8
0.0796
0.5256
0.0025
0.0490
0.0044
0.0045
0.0896
DML
0.7614
0.0024
0.0468
Q-MIN
0.0019
0.0317
0
0.0017
0,0135
0.0007
0.027L
0.00L5
0.0309
0.7184
0.000L
0.0227
Gen
Q-SL
0.5309
0.0003
0.0186
0.2679
0.0003
0.I188
0.4612
0.0006
0.0249
0.2474
0,0003
0.0189
0.942
0.0055
0.0235
T-TVEE
0:2001
0.0005
0.0218
0.3555
00000
B,200
0,302
0,0006
0.0259
0.467
0.0005
0.0216
0.9806
0,0005
0.0241
(n - 5000)
Q-T
0.748
0.0011
0.0248
0.7897
0.001
0,026
0.8191
0.0011
0.0254
0.5491
0.001
0.0244
0.9
0.0011
0.0253
DML
0.7733
0.0005
0.0198
Q-MIN
0.0015
0.0213
0,0033
0.223
0,101
0.0004
0.0193
0.0011
0.0221
0.0522
0.000S
0.0168
Gea
Q-SL
0,7922
0.0001
0.0123
0.4991
00002
0.013
0.0003
0.0LG7
0.0275
0.0001
0.0127
0.0051
0.0033
0.0164
Q-TYAE
0.0163
0.8935
00002
0.151
0,0025
0,0003
0.017L
0.3575
0.0002
0.0156
0.0191
0.0172
(n 10000)
Q-T
0.0717
0.000S
0.0174
0.0615
0.0006
0.0177
0.0663
0.0005
0.0176
0.0636
0.0005
0.0176
0.04
0.0005
0.0175
DML
0,5017
0.0003
0.0137
Q-MIN
0.0
0.0647
0.809
0.0
0.0451
0.465
0.0
0,0390
0.400
0.0
0.0515
0.411
0.0
0.614
0.844
IHDP
Q-SL
Q-TVAE
0.0
0.0440
0.317
0.0
0.0
0.0122
0.6053
0,404
0.346
0.00
00134
0.03400
0.405
0.344
0.0
0.0090
0.346
0.808
0.0
0.0
0.191
0.625
0.896
(n - 747)
Q-T
0.0
0.0124
0.630
0.0
0.0121
0,422
0.0
0,0118
0.424
0.0
0.0125
0.426
0.0
0.0293
0.436
DML
0.0
0.0129
0.441
1.1 Shapley Results
In addition to the interaction plots given in the text, here we provide results for the regressor performance and the predictor impacts.
1.1.1 Regressor Performance
Table show the performance of the random forest regressor in predicting the three outcomes (MSE, ATE estimate standard error, and p-value). These results help us understand whether there is any information in the set of predictors which is useful in predicting the outcome (and therefore, in turn, whether there exist any potentially meaningful patterns). We provide results for the fraction of explained variance, R², and MSE. The results are averages over a 10-fold valls-vald evaluation. is useful to interrogate this table first to understand whether any further investigation is needed if the predictive performance is poor,
43

---

Published in Transactions on Machine Learning Research (03/2023)
Table 5: Unmarginalized results over a restricted set of model variations including the Double Machine Learning (DML) method (Chernozhukov et al.] 2018). Average absolute error on the ATE (acATE) and standard error (s.e.) (lower is better) and Shapiro- Wilk test p-values for normality (higher is better) are Waip
provided ad computed over 100 simulations. Best results those competing across all three dimensions, although aeATE is prioritised. Bold indicates the best result for each algorithm. Multiple methods may perform approximately equally well.
Datawet
Q Model
U-Base soATE
G-SL U-multi
J soATE
P ncXTE
G-SL+U-aub
G-MNNUUmull
scATE
G-MN++Usuub
P noATE
Q-MN
0.3357
0.0385
0.0286
0.251
0.0379
0.0275
0.0325
0.0700
0.0613
0.1368
0.0378
0.0289
0.9045
0.0768
0.0605
Q.SL
0.00118
(.9229
0.1877
0.0351
0.0442
0.0391
0,0204
0.0758
0.0578
0,0375
0.0443
0,0888
0,0222
0.9758
0,0577
LP (v1)
Q-TVAR
0.0036
0.06-2
0.0364
0.1008
0.0552
0.0349
0.0728
0.081
0.0632
0.1043
0.0584
0.0346
0.0029
0.0805
0.059
(m-50))
Q-T
0.2005
0.063
0.0422
0.7785
0.0613
0.0134
0.4219
0
0.0552
0.7178
0.0018
0.0429
0.1244
0.0734
0.0626
DML
0.1225
0.0635
0.0433
Q-MI
0,064
0.0217
0.0144
0.0037
0.0178
0.0144
0.6148
0,0248
0.0173
0.011
0.0182
0.0133
0.5032
0.0237
0.0179
Q-SL.
0.0216
0.0163
0.013
0.0519
0.0178
0.0155
0.5245
0.0237
0.0175
0.0657
0.0166
0.0129
0.2467
0.0217
0.0164
LF (v1)
Q-TVAE
0,2875
0.0228
0.0179
0.2168
0.0202
0.0168
0.0859
0.0243
0.018
0.16
0.0213
0.0173
0.455
0.0227
0.018
(n=5000)
Q-T
0.9761
0.024]
0.0183
0.9719
0.0245
0.0189
0.9311
0.0254
0.0198
0.9919
0.0246
0.0194
0.6963
0.0254
0.0190
DMI
0.7166
0.0278
0.0199
Q-MIN
0.0188
0.0175
0,0189
0.0165
0.0109
0.0181
0.0189
0.2521
0.0159
0,0135
Q-SL
0,9998
0.012
0.0095
0,9372
0.0129
0.0097
0.013
0.016
0.011
0.9926
0.0118
0.0092
0.3251
0.0149
0.0108
LP (v1)
Q-TVAR
0.2716
0.0153
0.0118
0.3701
0.0143
0.0085
0.5634
0.0161
0.0112
0.6847
0.01488
0.0104
0.1576
0.0154
0.0112
[m= (0=10000)
Q-T
0.3408
00113
0.3254
0.0179
0.0117
0.3335
0.0177
0.0116
0.3688
0.0177
0.0117
0.4915
0.0171
0.0116
DML
0.777
0.0174
000339
Q-MIN
0.4621
0.0926
0.0248
0,1313
0.0326
0,0241
0.016
0,0801
0,0899
0,2789
0.0344
0.0238
0,003
0,0724
0,0688
Q-SL.
0.001
0.0396
04328
0.0001
0,0384
0,0342
0.0086
0,0723
0.0572
0.0072
0.0385
0.0345
0,0001
0,070S
0.0564
LF (v2)
Q-TVAE
0.0806
0.0708
01354
0.0662
0.061
0.0312
0.0055
0,0813
0.0679
0.0203
0.0655
0.0055
0
0.076
0.0601
(m-50))
Q-T
0,333
0,0629
00002
0.735
0
0,0344
0.0103
0.0893
0.0087
0.8407
0.0208
0.0380
0.0079
0,0703
0,0685
DML
0.0167
0,06M4
6,6368
Q-MIN
0.02F7
0.02447
0
0.12477
0.0225
0.8417
0.0251
0.0181
0.0283
0.02222
0.6607
0.0222
0,0220
Q.5L
0.7856
6,0150
0.9039
0,0219
0,0158
0.9209
0,0249
0.0179
0,702
0:0227
0.0158
0.5785
0.0241
0,0188
LP (v2)
Q-TVAR
0.908
0,0225
0.0159
0.6827
0.0218
0.0168
0.7211
0.0242
0.0184
0.4674
0.0224
0.0147
0.2977
0.8233
0.018
(m=5000)
Q-T
0.7928
0.0249
0.0192
0.8975
0.0252
0.0193
0.7337
0.0250
0.0197
0.7734
0.0257
0.0198
0.0001
0.0271
0.0254
DML
0.824
0.0275
0,0206
Q-MN
0.0287
0.0206
0
0.0295
0.0224
0.8711
0.018
0.0151
0.0288
0.0192
0.1258
0.0246
0.0185
Q-SL.
0.8728
0.0162
0.0113
0,7803
0,0170
0.117
0.9735
0.0175
0,0125
0,5777
0.0156
0.0115
0.0173
0.0118
LP (v2)
Q-TVAB
0.5869
0.0163
66116
0.3810
0,0158
00113
0.0688
0.016
0.0115
0.3413
0.0155
0.0108
0.2405
0.0165
0.011
10000)
Q-T
0.617
0.0178
0.0124
0.7051
0,0184
0.0131
0.5027
0,018
0.0151
0.3855
0.0179
0.0127
0.2885
0.0188
0,0133
DMI
0,7006
Q-MIN
0,0012
0.0033
61899
0.0417
015851
0.0004
0.0712
0,0559
0.045
0.0309
0.1512
0.062T
Q-SI.
0.9273
0.0512
040350
0.892
0.05m
0.086
0
0.0619
0.0444
0.7328
0.0496
0.0362
0.0503
0.046
Gen
Q-TVAE
0.0083
0.0221
0.2201
0,0679
0.0303
0,0182
0.0003
0,0754
0.0619
0.103
0.0022
0.0196
0.0583
0,0513
(n-500)
Q-T
0.0305
0.0389
0.0263
0.0733
0,0999
0.0247
0.1224
0.056
0.0436
0.0G15
0.0375
0.0248
0.0002
0.0513
0.0405
DML
0.6547
0.0387
0,0248
Q-MN
0.037
0.0307
0
0,0349
0.0344
0.7922
0.0206
0.0155
0.0312
0.0291
0.2530
0.0174
0.0137
Q Q.SL
0.514
0.0131
0.0107
0.5157
0.0137
0.0111
0.0922
0.0187
0.0148
0,2609
0.0133
0.0111
0.8739
0.0175
0.0138
Gen
Q-TVAE
0.9298
0.0130
0.9347
0,0195
0.0344
0.4052
0.0185
0.015
0.5204
0.0176
0.0130
0.6550
0.0187
0.0145
(a=5000)
Q-T
0.748
0.0276
0,0178
0.8857
0.0267
0.0178
0.8492
0.0275
0.0183
0,8432
0.0209
0.0170
0.0702
0.0291
0.0208
DML
0.808
0.0188
0.0127
ANN
0.0329
0,0305
0,0329
0.0109
0.0125
0.014L
03027
0,002
0.3189
0.01
0.0011
Q-SL.
0.9635
0,0083
0.00G8
0.068G
0.0085
0.0003
0.2335
0.0111
0.0104
0.9994
0.0082
0.0062
0.6388
0.0102
0.0077
Gen
Q-TVAR
0.6851
0.0134
0.0092
0.5345
0.0116
0.0066
0.6109
0.0117
0.0111
0.6207
0.0113
0.0076
0.7070
0.0115
0.0083
(a-10000)
Q-T
0.3724
0.0196
B.000
0.3895
0.0184
0.0106
0.4392
0.0191
0.0103
0.5078
0.0188
0.0106
0.3961
0.019
0.0105
DML
0.9919
0.0145
0.0091
Q-MI
0.2142
0.3337
0.2066
0.2972
0.1673
0.2141
0:244
0:476
0.3005
0.7104
Q-SL
0.1474
0.1480
0.1292
0.1253
0
0.1564
0.1277
0.1543
0,2658
0.1852
0.2405
DDP
Q Q-TVAE
0.1125
0,0772
0
0.0199
0.0071
0
0.0912
0,0707
0.1289
0,1890
0.2088
0.2032
(u=767)
Q-T
0.0968
0,0889
0.0667
0.0677
0.0655
0.1118
0.1742
0.1328
0.1681
DML
0.0919
0.0667
there is no point explaining the regressor's behaviour; if it is good, then it is worth investigating what the regressor is using to achieve that performance.
From Table it can be seen that all outcomes were somewhat well predicted, with the arguable exception of the MSE, which had large standard errors for the R² and fraction of explained variance suggesting that test data in some of the folds in the 10-fold cross-validation scheme were poorly predicted, whilst others were predicted well. Even though R²2 and explained variance are not reliable metrics in non-linear regression lasks, they can nonetheless be seen that a relatively high R2 is achieved, indicating that information about the MSE is predictable from the set of predictors. This was especially true for the prediction of the Shapiro-Wiil Wilk test p-values and ATE estimate standard error, which both had average fraction of explained variance and average R² greater than 0.7.
44

---

Published in Transactions on Machine Learning Research (03/2023)
Table 6: Meta-analysis results for random forest regression performance for MSE, standard error of the estimates (ATE s.c.) and p-values as the outcomes. Results for R², explained variance, and MSE are given as the mean + the standard derivation across the 10-fold voss-valioatio procedure.
Outcomes:
MSE
ATE
p-value
Rd
0.47#0.667
0.97=0.022
0.4720.667
Explained Variance
0.48.10.65 0.88 0.000510.0006
0.97:00027
0.48100668
MSE
0.000310000004
0.0005L00.0066
semple_siee
sampl_5e0
QLR
QUR
GAN
Umm. mew. w
Umuli-En meneo
Usab
G0
G0
GAN
delare
U-0aee
Ulase
G.
9.0
0.5L
3.7
QT
G4Dnltng
0 Dnohng
GLL
G
O.TVAE
TWE
6 Onoreg
0 Dnotnng
088
GP
G.P
0000
noar (BHAP valua)) (impart on outpul magniaale)
0000
0.00
5038
3300
-0.25 0.000
80-AP unlun (impurt on modted (uidpe))
0.08
8.360
0076
0.100
3.125
0.700
Figure 17: Shapley predictor random forest regressor impact results for the MSE.
1.1.2 Shapley Predictor Impact
Given the regression performance results indicate some patterns exist in the data, let us now turn to the Shapley explainability results. For each of the three respective outcomes (MSE, ATE estimate standard error, and p-values), Figs. 17 18 and 19 depict the global predictor impacts on the regressor (left-hand plots) as well as the per-prediee per-predictor, per-datapoint impact on the regressor (right-hand plots).
The right-hand side plots in Figs. 17 19 are useful in visualising the spread of impact, and in which direction this impact is (i.e., does it push the prediction of the outcome up down in value). For example, in Fig. 17 one can see that the dataset is the most important predictor for predicting MSE, followed by Q-CFR, G-MN, and U-sub. We see that most of the impact (for all predictors) is clustered around O, which indicates that they are stable in terms of their relationship to MSE. However, the right-hand side plot shows that some outliers exist, with certain combinations yielding higher (or, to a lesser extent, lower) MSEs.
In Fig. 18 for ATE estimate standard error we again see dataset the most important predictor, followed by Q-CFR. and sample size. There is a larger variance in the impact these predictors have on the outcome than for MSE. Finally, looking at Fig. 19 for the p-values, we see heavy dependence on sample size. Note that the p-value outcome here is a function of a different sample size to the one used in the regression - it is computed based on the number of simulations which is fixed for all experiments. Thus it is still interesting to note that sample size was important as a determinant of the Shapiro-Win Wilk test for normality p-value.
Perhaps more interesting than these predictor impact results are the interaction plots, which are given in Section 6.4
45

---

Published in Transactions on Machine Learning Research (03/2023)
pide
Usub
datare
dabase
GMN
QUR
Ucees
Uona
Cnollr
G0ndrng
Um. ve
Umut-In-mennwr
0.8
0-8
0.0
040
0.5L
0-8L
3-T
0.1
G.5L
G8
0 Dnotne
0.TYAE
QMN
CNN
U-Ba60
USase
0.0
600
GLR
0.RR
G.P
G.P
0.22
(/HHP w (impuer nupput megsisate)
(empact
Figure 18: Shapley predictor random forest. regressor impact results for the ATE estimate standard error.
size deaset
dabewl
QD
GD
Ummub mommarr
Umubiih-raamar
G.D
GNN
QAN
0-T
G-NN
G-NN
U-cnaa
Ukna
0.8
Q5
Uni
QLR
QLR
Q.norreg
Drotreg
Q-SL
Q88L
T.TE
0 TWEE
0.SL
G.8L
G1R
C-7
G.P
Dase
G Drslig
GDrcbng
3.0
60
3.6
0.3
010
-2.
60
mas ([HHAP wilel Gorpac mndoll magelal)))
SHAP vndua (impart coalpaa))
Figure 19: Shapley predictor random forest regressor impact results for the Shapiro Wilk test p-values.
1.2 Alternative Perspectives
In the main text we provided summary results by estimating the probability that a particular Q (outcome). G (propensity), or U (update stop) method would result in a performance advantage. This was done because the number results was large, making it difficult to judge the efficacy a method in isolation. In Figs. 20 26 we provide the complete results for each of the seven dataset variants: LF (v1) with (500, 5000, 10000).
46

---

Published in Transactions on Machine Learning Research (03/2023)
LF (v2) with (500, 5000, 10000), and the IHDP dataset n 747. For each figure provide the comparison of each Q-method with each G- and U-method, and include a red dashed line to include the base method (just the Q-method without the IF update step) for comparison.
47

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 20: LF (v1) n = 500 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (left sub-colume column) and each of the possible update methods U. The base performance (no update step and therefore no G or U) given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a margina marginalization over the other dimension. For instance, for Q-D (DragonNet), the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Graph best viewed in color.
3 Melhodd
GIR
Malhods
Medhadd
64
Maitods
G-SL
G Medhadd
CN
08
Ginggn
66
Methods
CAN
Method
0 Melhods
GONtiu AM
GEL
140
G Malhoot
68L
U Mathods
Unjatt
Mod
Dauii
O Mehrods
40
GLR
0000
48

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 21: LF (v1) n = 5000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y- (y-axis) for each of the possible propensity models G (left sub-column) and each of the possible update methods U. The base performance (no update step and therefore no G or U' given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a margina marginalizatiin over the other dimension. For instance, for Q-D (DragonNet). the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Graph best viewed in color.
Methode
G
C Meimods
Methods
Unmo
Louu
GUROUR
Meshods
0.0
CLSL
Mehrods
GD
6.bgpnn
02
Methods
0-0ngres
05L
GLR
Melhod
Merhadd
Q-T
GG Methoode
03000 GD
G-NN GLR
Upids
G-TWAE
G Methoots
GRL
O-Orotreg
Methods
Vaul
Methods
0-0
Menodd
COLL
nnto
0.25
0.004
0.000
49

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 22: LF (v1) n = 10000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (left sub-column and each of the possible update methods U. The base performance (no update step and therefore no G or U) given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a marginalizatiin over the other dimension. For instance, for Q-D (DragonNet). the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Graph best viewed in color.
D Ueteeds
GMN
400
G
Uieis
Undi
OLR
GNN
018
CWN
CBL
U-bib
C.S
Q-T
O-TYE
BLR
UHUR
GLR
60
Q-0
6D
anto
0.08
0006
0.34
MSE
50

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 23: LF (v2) n = 500 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (left tub-column. and each of the possible update methods U. The base performance (no update step and therefore no G or U) given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a marginalization over the other dimension. For instance, for Q-D (DragonNet), the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Best viewed in color.
3 Melhodd
34N
G-theting
GLR
9
Malipa
Unuai
O
Medhadd
GAR
Maitods
G-SL
O Medhadd
G-bhetigg
Methods
Method
J-mi
Melhods
Dhoinggs P
068
Sqqun
G Malhoot
BUR
6
Voless
U Malhods
Mod
GL
60
Mab
O Mehrods
5-N
8.SL
40
0036
0000
51

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 24: LF (v2) n = 5000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y- (y-axis) for each of the possible propensity models G (left sub-column) and each of the possible update methods U. The base performance (no update step and therefore no G or U' given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a marginalizatiin over the other dimension. For instance, for Q-D (DragonNet). the "U-ones" point is an average result for the onestep ep update process, using all possible propensity models G. Graph best viewed in color.
3 Melhodd
0LR
GdL
Dogus
Malhods
Methadd
GAN
Goyying
Maitods
G Medhada
OUR
0,L
Methods
25
6.Dgeng
0 Melhods
Usgies
Uniuge
G Malhoot
U Malhodd
Mod
Hith
Dolles
Ookggoo 13
e Mehrods
0000
52

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 25: LF (v2) n = 10000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (left sub-column and each of the possible update methods U. The base performance (no update step and therefore no G or U) given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a marginalization over the other dimension. For instance, for Q-D (DragonNet) the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Best viewed in color.
3 Melhodd
GLR
0
Malhods
Methadd
Diiilgg
Maitods
G Medhada
GLR S
Uebb
Methods
GER
C-N
0 Melhods
0UR
COLL
GDWN
Oapus
G Malhoot
OUR 6o
Unbh
U Malhodd
Mod
e Mehrods
6
53

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 26: IHDP n = 747 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (left sub-column) and each of the possible update methods U. The base performance (no update step and therefore no G or U') is given as a horizontal dashed red line. Because we undertook all combinations of G and U each point represents a marginalization over the other dimension. For instance, for Q-D (DragonNet), the "U-ones' point is an average result for the onestep update process, using all possible propensity models G. Best viewed in colour.
Mellods
Jomun
Jo0bb
G Mejtoos
94h
P4NN
Dob
0.000g
G Meltoots
20
Jemani
Medhoch
G Mainccd
C-NN
0.30
6LR
Method
Jofbe
Medtodss
Meelodd
20
9LR
Melhoss
G Merhodd
Jmmi
Sebb
G Methoöi
PLR
R WN
Shethodd
G Mellooa
BUR
G
0.25
0220
0.55
0.10
0.5
54

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 27: General graph = 500 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each the possible propensity models G (left column) and each of the possible update methods U. The base performance (no update step and therefore no G or U') given as a horizontal dashed rod line. Becanse we indertook all combinations of G and U, each point represents a marginalization over the other dimension. For instance, for Q-D (DragonNet). the -U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Best viewed in color.
3 Vernode
Methods
Ueed
Meihodk
Maethade
Danaa
02 Greiing
G Medhode
CD
mil
U Methoos
GLR
G-LL
Uglb
GDvooggg
C Mettoote
CAN
GUR
Smuid
Meiloods
Methaos
GAN
Jmmb
U Melhodl
Ufborss
Medrodd
OUR
Co
GBN
Medlrod
O Menacde
50
68L
D-Dooooo
200
204
3.03
3,22
301
MBE
55

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 28: General graph n = 5000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axis) for each of the possible propensity models G (Jeft sub-column and each of the possible update methods U. The base performance (no update step and therefore no G or U') given as horizontal dashed rod line. Bocanse we undertook all combinations of G and U, each point represents marginalization over the other dimension. For instance, for Q-D (DragonNet). the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Graph best viewed in color.
3 Vernodd
GL
Uhab
Methods
Meihodk
GANN
Maethade
G Medhode
Jerufl Counaa Co
G-LL
U Methods
Usab
G-hhooog MNN
OLR
'900
0.
Done
C Mettoote
Meiloods
Uadb
Metroit
CANN
EK
6.
0.5L
Ueub
Umabi
U Melhodl
Medrodd
CLR
Medlodd
O Menacde
AMN
C4R
Dridtnn
0-38
2005
204
0.33
3322
301
MBE
56

---

Published in Transactions on Machine Learning Research (03/2023)
Figure 29: General graph = 10000 results. For each outcome model Q (x-axis) we plot the corresponding Mean Squared Error (y-axiss for each of the possible propensity models G (Jeft sub-coluun column) and each of the possible update methods U. The base performance (no update step and therefore no G or U') given as horizontal dashed rod line. Bocanse we undertook all combinations of G and U, each point represents marginalization over the other dimension. For instance, for Q-D (DragonNet), the "U-ones' point is an average result for the onestep ep update process, using all possible propensity models G. Best viewed in color.
Uone
3 Vernodd
OUR
Dones
Utab
Methods
Meihods
G-NN
GLR
Maethade
Umult
G Medhode
GAR
Unet
U Methods
0-S
OUR cAN
C Mettoote
G-Orodiing GAR Gaa
Canis
Meiloods
Unud
Umati
Metroit
CANN
GLR
Looes
U Melhodl
Medrodd
CLRS
GOL
Medhodd
O Menacde
AIN
EL
2005
204
0.33
3,22
301
57


Reviews:
[{"review_id": "sb2UvKFlgq", "reviewer": "Reviewer_b9xN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nThe paper presents three contributions to average causal effect estimation by using neural networks. The proposed method, MultiNet, is to mimic the performance of an ensemble with a single NN. Moreover, the paper provides a comprehensive evaluation of the recent semiparametric methods. The experiments demonstrated the performance of the proposed methods. The structure of the paper is good and easy to follow.\n\n**Weakness**\n\nThe first contribution of the MultiySet update method is close to robust double machine learning (DML). It would be better to compare MultiySet with the standard DML method in the main text and experiments.\n\nThe causal graph in Fig 4 is too simple, i.e. there is not a variable that is a cause of $Y$, but not a cause of $W$. It would be better to consider a more generalised case. In my experience, more variables are involved in the covariates, and large biases of estimation will be obtained.\n\n\n\nOverall, the paper is well presented and can be considered to accept after a minor version.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Causal inference plays a critical role in decision-making and data analysis but it is challenging to be evaluated due to there being no ground truth on causal effects or data generation process. The paper studies the important problem of causal inference, i.e. estimating the average treatment effect from observational data. The paper makes a solid contribution in the area of causal inference. One of the conclusions, i.e. ``practitioners need to check the consistency of their findings'' is a pertinent suggestion.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Minors**:\n\n-> Section 6.2.4: \n\n-> Section 6.3: Figs. 18-24", "broader_impact_concerns": "The paper provides a good estimator for causal effect estimation from observational data.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2xrmTGCRB1", "reviewer": "Reviewer_YQ4B", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, the paper is well organized. With various estimation methods, it does necessary to have a paper to compare them and provide a comprehensive analysis to show under which cases one method is preferred to another.\n\nHowever, as a paper focusing on empirical evaluation and providing suggestions for users on how to choose estimators, I do not think the current analysis is enough. For instance, the authors should consider different datasets that contain various graph structures and various function types, and consider the cases where latent confounders or selection bias exist. Also, more detailed suggestions should be provided to the users on how to choose estimators.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper focuses on an empirical evaluation of influence functions for the average treatment effect estimation. Moreover, it proposes a new estimator \"MultiNet\" and a variation on the semiparametric update step \"MultiStep\".\n\nSpecifically, the authors empirically evaluate different estimators of the outcome model Q, different estimators of the treatment model G, and different update processes U. The evaluations are performed on three datasets with varying sample sizes and with three evaluation criteria.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "This paper focuses on empirical evaluations, so more analyses are needed.\n\n1. The authors should consider different datasets that contain various graph structures and various function types, and consider the cases where latent confounders or selection bias exist. \n\n2. Also, more detailed suggestions should be provided to the users on how to choose estimators.", "broader_impact_concerns": "No.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "wdaAZXaVdc", "reviewer": "Reviewer_kwkh", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strength\n\n- The writing logic of Sec.s 1 -3 is pretty good. It makes the target and contributions of the work clear.\n- The comparison experiments are appreciated, which require lots of effort.\n\n# Main concern 1\n\nOne of the main contributions is the empirical analysis and comparison for a good understanding of the consistency of different methods. However, the current Sec. 6 and experiments are about the comparison of different methods and report the best-performed methods. It is not unlikely that I missed some key take-away messages and key parts in Sec. 6. But regarding answering the question about a good understanding of the method consistency, it cannot be answered by reporting the better performance methods, but requires more investigation and insights to answer the question. From the current study of work, maybe the authors could make it clear what is the take-away methods regarding a good understanding of the method consistency. How is it different from the current understanding in the domain? And why is it important and interesting for the domain? Please elaborate on them in the rebuttal and paper. \n\n# Main concern 2\n\nThe introduction of the proposed methods in Sec.s 4 and 5 needs to be further improved and polished. \n\n## Sec. 4\n\nit is not clear what is H in Eqn. (6). And it is not clear how the proposed method works. Importantly, it would be helpful for me to introduce the proposed algorithm in a complete manner for the readers, e.g., an algorithm box. \n\nMoreover, it claims the improved robustness of using the two criteria (A) and (B) of IFs. However, it is lack of support for the point from the paper, theoretically and empirically. If there are, please elaborate in the rebuttal and paper. \n\n## Sec. 5\n\nRegarding MultiNet, the idea is simple and clear that is taking the outputs from each layer and weighting them together with linear coefficients to train the model. However, it is not clear to me where and how IF contributes to the method.\n\nIn Sec. 5.1, it mentions that the IF can worsen the performance of NN. Even though this can be mentioned in other works, it is important to be self-contained on the important points. Thus, what is the reason for it? How and why the proposed method can deal with the problem?  Moreover, what motivates the proposed method to deal with the problem?\n\nIn Sec. 5.2.2, it claims the benefits for overfitting and convergence issues, which is lack of support. If there are, please elaborate in the rebuttal and paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The work focuses on the average treatment effect estimation (ATE). It highlights one of the challenges, that there is limited access to the ground-truth. Thus,\n\n(1) it is essential to have a good understanding of the consistency of typical methods.\n\n(2) it is necessary to have an unbiased estimate with statistically valid confidence intervals.\n\nRegarding (1), the work provides an empirical comparison. As for (2), the work proposed to use influence functions for ATE and propose an update step based on two constraints of IFs. Moreover, it proposes the method MultiNet.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "Besides the mentioned concerns in the previous part, please also elaborate on the following points. \n\n# Clarification\n\nIn Sec. 1, it mentions the evaluation from 3 aspects (a)-(c). It requires justification and motivation for them. Especially, why these three are important and sufficient for the evaluation? Are there any other aspects to consider? Since the paper claims an empirical contribution, this is an essential point.\n\nIn Sec. 3.1, please elaborate on “Furthermore, even with correct specification of our plug-in, …., joint distribution).”\n\nIn Sec. 3.2, from my understanding, IF is used for capturing the bias of the estimated concerned quantities. According to Eqn.(4), it is based on the propensity score model and the outcome models. However, based on the current introduction in Sec.s 3.1 and 3.2, the Q hat in Eqn. (2) and (3) is modelling the conditional probability expectation. I would like to suggest making it clear how it is related to Eqn.(1), which is the target quantity (e.g., maybe make it clear they are  the results with do calculus?)", "broader_impact_concerns": "Looks good. I didn't have any concern on the ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper aims to tackle the challenge that there is no straightforward evaluation of treatment effect estimation methods. To this end, the authors check the consistency of methods by evaluating different estimators, including a newly proposed one, on different datasets. The main finding is that method performance depend on dataset properties, such as the sample size and data generating process.

The paper makes a good contribution for practitioners to understand the consistency of causal effect estimation methods and how to choose among them. The empirical evaluation and presentation has been significantly improved according to reviewers' comments. 

The authors have promised reviewer YQ4B to add more empirical studies on more general graph and different level of nonlinearities, but this has not yet been done. I would suggest "acceptance with minor revisions". The authors should complete this part before the paper is finally accepted.
Decision: Accept (poster)