In many clinical trials or social science studies, comparison of treatment effect between two groups is often regarded as an important part of the research. This effect can be measured by using difference or ratio between each treatment group, and average causal effect (ACE) is often of primary interest since individual treatment effect is difficult to obtain especially when the study is not randomized. Under conventional randomized experiments, the treatment effect can be easily evaluated, since all factors except the treatment group (confounding factor) can be controlled by experimenters. In an observational study, however, many confounding factors are present and they may distort the causal effect severely, because the distribution of covariates among subjects between treatments can be different (Rubin, 1978).

In such cases, the inverse probability weighting (IPW) method can be used to balance the propensity score distributions between two treatment groups by dividing the observed outcome with the estimated propensity scores (Austin, 2011). The resulting estimator may approximate the effect as if it were obtained from a randomized experiment, but it requires a well-posited propensity score (PS) model for stable analysis (Rosenbaum and Rubin, 1983). However, problems due to model misspecification often arise in real applications, since the true propensity model is hard to know exactly in an observational study. Alternatively, the augmented inverse probability weighting estimator (AIPWE), a hybrid approach of regression and IPW method, has been extensively studied as a generalization of the IPW method (Robins et al., 1994; Kang and Schafer, 2007; Funk et al., 2011). AIPWE is well-known to be doubly-robust (DR), because it still produces a consistent result if either the regression or IPW model is correctly specified (Tsiatis, 2007).

However, AIPWE could be worse than a single IPW method when both models are incorrectly posited (Kang and Schafer, 2007). There have been many studies to overcome the limitation of this DR estimator. In this case, machine learning techniques are a natural choice to accommodate potential model misspecification for causal regression modeling (Lee et al., 2010; McCaffrey et al., 2013). Gruber and van der Laan (2010) proposed a targeted maximum likelihood estimation (TMLE) that is another semiparametric DR method, equipped with nonparametric super-learner (SL) algorithms. Using only SL models can infer false confidence intervals if severe bias occurs and one of the two models is incorrectly specified. Therefore, TMLE takes a semiparametric approach to the SL model to compensate for this limitations to have an overall unbiased estimator at all levels (Schuler and Rose, 2017). Causal forest (CF) (Athey et al., 2019) is a newly developed estimation method, based on ensemble models for causal inference. This method is developed from decision trees and random forest to estimate the treatment effect. It generates a large number of trees for each treatment group, then provides appropriate weights and variations for each group and finally generates an ACE estimator.

Despite various proposals for improving causal estimators, researchers may suffer from inadequate performance comparisons between these methods. In this article, we present an empirical study to compare several recently developed causal inference models for estimating ACE. We briefly overview recent developments in causal analysis, including regression-based estimator, IPWE, AIPWE, TMLE and CF and explore which methods perform better under various situations via simulation experiments. Our simulation study covers a wide range of scenarios that may change ACE and investigates existing and newly adapted machine learning models. We replace the regression with generalized smoothing model (GSM) and IPW models with several machine learning models with correct, mis-specified and omitted models, and show how to have lower bias and average error rate in a wide range of scenarios.

The rest of the article is organized as follows. Section 2 presents basic notation and underlying assumptions and summarizes various recent causal effect estimating methods. Section 3 provides extensive simulation results that compare performance improvements with different combinations of the resulting regression and propensity score models. In Section 4, several causal models are applied to the percutaneous coronary intervention (PCI) data for illustration. Section 5 provides a short explanation and discussion of our results.

2.1. Preliminaries: Data and assumptions

In causal inference, we are often interested in clarifying the causal relationship in observational data. However, to allow for a precise presentation of the causal treatment estimators, we need additional notation and assumptions for the observed data (Hernán and Robins, 2020). For i = 1, 2, . . . , n individuals, let Y_i and A_i be the observed outcome and the treatment or exposure group variable, respectively, and X_i be the p-vector of observed pretreatment covariates. Then, the observed data are represented by {O_i = (X_i, A_i, Y_i), i = 1, . . . , n} and the data distribution follows an unknown distribution P. For simplicity, we assume that A_i takes a binary value, A_i ∈ {0, 1}, but the method can be easily extended to the case with multiple treatments (McCaffrey et al., 2013). Under the counterfactual framework (Rubin, 1978), each individual will have two potential outcomes, say Y(a), which represents the outcome variable that corresponds to the exposure factor A = a. Suppose that we are mainly interested in inferring the average causal effect (ACE),

which measures average treatment difference between Y(1) and Y(0). However, the causal effect for an individual generally cannot be estimated because of the fundamental problem of causal inference: each individual will receive either A = 1 or 0, but not both at the same time.

To quantify the causal effect in observational studies, we shall make the following three common assumptions in causal analysis (Hernán and Robins, 2020),

• (i) consistency: Y = AY(1) + (1 – A)Y(0),

• (ii) exchangeability: Y(a) ? A | X,

• (iii) positivity: 0 < P(A = 1|X = x) < 1 for all .

The consistency assumption states that there are no multiple versions of treatment outputs, which means that the mechanism used to assign the treatments does not matter and assigning the treatments in a different way does not constitute a different treatment. The exchangeability assumption says the set of observed pretreatment covariates is sufficiently rich, such that there exist no unmeasured or unknown confounders, not included in X. Since there is no way to test this assumption with observed data, analysts should try to include a broad set of covariates to reduce the potential risk that a confounding variable is omitted inadvertently. Finally, the positivity assumption states the observed data for both treatments should be available for any given X, such that group comparison will make sense.

The commonly used randomized experimental designs allow researchers to estimate the average causal effect estimates closest to the theory, in which we may consider

Unlike randomized experimental data, however, the treatment group is unevenly distributed or associated with underlying confounders in observational data, making it difficult to use τ? for predicting the causal effect. There are several ways to solve this problem, such as matching, stratification of sub-classification, adjustment by instrumental variables, and treatment weighting (Stuart, 2010; Robins et al., 2007). In the following, we review a few recent developments for clarifying causal relationships between treatment and outcome in observational studies and provide some guidance on how we can implement machine learning techniques to improve the quality of causal estimates.

2.2. Classical approaches: Regression and inverse probability weighting (IPW) estimators

The estimation of treatment effects may proceed in two stages by modeling the conditional outcome Q_a(x) = E[Y|A = a, X = x] and the propensity score g(x) = P(A = 1|X = x). In principle, it is straightforward to use statistical regression models for Q_a(x) and g(x). First, we consider estimating the average causal effect through regression models. Since the average effect under A = a can be formulated by , a regression-based causal effect estimate can be obtained by

where Q?_a(x) is a regression estimate of the conditional mean outcome Q_a(x). We can use standard parametric regression models or more flexible semiparametric or nonparametric methods to predict the outcome Y_i from the treatment and covariates.

In our simulations, we primarily employ a generalized smoothing model (GSM) by Helwig (2020), which estimates parameters by estimating the unknown smooth function under a multiple and generalized nonparametric regression model from the sample of data. The GSM algorithm is available in the R package, npreg, which offers several smoother types and will select the smoothing parameter basically based on generalized cross validation (GCV) according to the nominated predictor’s class without specifying the model formula.

On the other hand, the inverse probability weight estimator (IPWE) can be described as a difference of weighted averages that assigns inverse probability weights to each treatment group. We can estimate the propensity score g(x) by ?(x) and obtain the treatment effect as

Conceptually, IPWE attempts to fully adjust for measured confounders by balancing the confounders across levels of treatment with treatment weight. When a function for the propensity score is unknown, logistic and multinomial regression models are commonly used for binary and multiple treatments, respectively.

The propensity scores, when properly estimated, ensure that covariates are assigned proportionally across two treatment groups, and may control very distant outliers. In this case, machine learning methods, such as boosted ensemble models, have been shown to yield causal effect estimates with many desirable properties (Lee et al., 2010; McCaffrey et al., 2013). We shall also explore generalized boosted model (GBM) for estimation of the necessary propensity score weights.

2.3. Augmented inverse probability weighting estimator (AIPWE)

An alternative and improved estimator is the augmented inverse probability weighted estimator (AIPWE) that combines the properties of both regression-based estimator and IPWE (Robins et al., 1994; Kang and Schafer, 2007; Funk et al., 2011). AIPWE requires two models: the outcome model and the propensity model. The outcome model typically uses a regression model and must specify an outcome model formula that states the relationship between covariates and outcome variables within each treatment group. Likewise, the propensity model specifies a model that relates covariate effects with treatment allocation. Specifically, the AIPWE takes the form

where

Notice that the first term in ${\widehat{\mu }}_a^{\text{AIPWE}}\hspace{0.17em}(a=0,1)$ is equivalent to ${\widehat{\mu }}_a^{\text{IPWE}}$ in (2.3), while the second term augments the IPWE to reduce variability and improve estimate efficiency by appending the regression estimator Q?_a(x) in a structured way. The AIPWE is known to be “doubly robust” in that it only requires either the propensity or outcome model to be correctly specified but not necessarily both (Tsiatis, 2007). The “augmented” term could complement the bias with zero even when either the propensity model or the outcome model goes wrong. AIPWE can also be aided by machine learning methods; for example, we may use GBM for the propensity score model instead of standard logistic model and GSM for the outcome model instead of linear regression model. In the last decade, AIPWE has gained much popularity and has been extensively studied due to its efficiency, reliability and simple implementation. Recently, Choi et al. (2021) proposed a pseudo-value-based AIPWE when the response is subject to censoring and showed that the double-robustness property is still valid under model misspecification.

2.4. Targeted maximum likelihood estimator (TMLE)

TMLE (van der Laan and Rose, 2011; Schuler and Rose, 2017) is another well-established doubly-robust maximum likelihood method for causal effect estimation. Like AIPWE, the TMLE procedure should also specify both the outcome and propensity models. It begins by estimating the conditional mean models given the treatment and covariates, Q₁(x) and Q₀(x), and the propensity score model g(x) with a prespecified method. Then one can create so-called “clever” covariate, which is very similar to the IPWE format and given by

Then, we generate “targeted” estimates of the set of potential outcomes, incorporating the clever covariate as an adjustment to reduce the bias. Specifically, by treating the initial estimate ?_a = &Qcirc;_a(x) as an offset, fit the following model

to obtain the estimate δ? of the fluctuation parameter δ. h(·) denotes a link function in generalize linear model. Note that H(x) is the term to compensate for the estimation of the outcome model and relies on the variation of the parameters of interest and δ? is the maximum likelihood estimator obtained by regressing H(x) on Y with the offset ?_a. Finally, we can estimate the ACE by

where ${\widehat{\mu }}_a^{\text{TMLE}}=1/n\hspace{0.17em}{\mathrm{\Sigma }}_{i=1}^n{\widehat{Q}}_a^{*}(x_i)$, (a = 0, 1). Here, ${\widehat{Q}}_a^{*}(x)$ can be obtained from $h({\widehat{Q}}_a^{*}(x))=h({\widehat{Y}}_a)+\widehat{\delta }\hspace{0.17em}H(x)$.

TMLE is a two-step semiparametric procedure that solves the efficient influence curve estimating equation and thereby yields an efficient or at least locally efficient solution to the parameter of interest. In addition, combined with the super-learner (SL) algorithm (van der Laan et al., 2007), TMLE has demonstrated its utility as a very powerful tool for efficient substitution estimator. By using machine learning methods, TMLE can have many advantages, such that it can reflect the underlying features of data smoothly needless of assuming restrictive conditions on a functional relationship of the outcome and exposure mechanisms. For example, several ensemble methods are considered to achieve better performance and then one of them can be selected according to data based on cross-validation and the ensemble algorithm can be applied to the data simultaneously and synergistically. Super learner is the final weighted results aggregation that combines several ensemble learnings by their performance, and is known to work very well with minimal cross-validated mean squared error. It can be implemented with the tmle package in R and statistical inference can be made by using the estimator’s influence curve or bootstrapping.

2.5. Nonparametric causal forest method

Causal forest (CF) is a causal inference learning method as an extension of Breiman (2001)’s random forest (Wager and Athey, 2018; Athey et al., 2019). In random forest, data are repeatedly split in order to minimize prediction error of an outcome variable. Causal forest is built similarly, except that instead of minimizing prediction error, data are split in order to maximize the difference across splits in the relationship between outcome variable and treatment variable. This is intended to uncover how treatment effects vary across a sample. Additionally, it can exploit the large feature space characteristic of big data and abstract inherent heterogeneity in treatment effects.

Fundamental idea of CF is constructing partitions of covariate act as sample came from randomized trial, which can be implemented as follows. Suppose there exist L partitions of covariate space (denoted by L-leaves) given an existing split rule. Furthermore, the leaf L(x) with test covariate value x is assumed to sufficiently small so that (Y_i, A_i) in i ∈ L(x) behave as they were gathered from a randomized trial. One can approximate the average treatment effect for any x ∈ L as

Then, the corresponding ACE estimate with causal tree can be obtained by ${\widehat{\tau }}^{\text{Tree}}(x)={\widehat{\mu }}_1^{\text{Tree}}(x)-{\widehat{\mu }}_0^{\text{Tree}}(x)$. Likewise other random forest, one can use a bagging algorithm, aggregating the results of multiple causal trees from bootstrapped samples. Denote ${\widehat{\tau }}_b^{\text{Tree}}$ as an ACE estimate from the bth causal tree with bootstrapped samples for b = 1, . . . , B. Then, CF aggregates ${\widehat{\tau }}_b^{\text{Tree}}$ and yields the final causal estimator by averaging their scores, that is,

The CF algorithm is available in the R package, grf, where double-sample (validation data and training data) tree is the default base model; Athey et al. (2019) for further applications.

Unlike conventional tree algorithms, CF algorithm is not based on minimizing the mean squared error in a regression tree. Since the true τ is not observable, we cannot calculate the difference of τ? – τ. We instead calculate the conditional average treatment effect (CATE) for covariates in each final leaf L as (2.6) and (2.7). According to Athey and Imbens (2016), one can optimize CF algorithm by maximizing the variance of τ?(x) rather than minimizing the mean squared error in a regression tree. Therefore, if the variance of τ?(x) is maximized during the splitting step, the value of the splitting criteria can be selected as an optimal value (Gulen et al., 2020). Although the variance can increase due to smaller sample sizes, causal estimates can be free of overfitting and bias, as in random forest. Further, one can control the increased variance to some extent by more repeating the bagging procedure (2.7). According to Wager and Athey (2018), τ?(x) is a consistent, asymptotically normal estimate of CATE.

In this section, we provide empirical simulation results to evaluate the finite-sample performance of several causal inference methods under a wide range of scenarios. For each outcome model (m) and propensity score model (π), we consider three scenarios: (i) correct-model specification, (ii) incorrect-model specification and (iii) omitted-model specification. This setting yields 9 possible combinations of outcome and propensity models. Each model also involves some nonlinear and interaction terms, which are included to explore the performance of machine learning algorithms for causal inference in more realistic settings.

Our simulation involves three covariates, Z = (Z₁, Z₂, Z₃)′, which follows a multivariate normal distribution with mean zero, unit variance and Cor(Z_i, Z_j) = ±0.2, (i ≠ j). The correct-model case involves the correct outcome variable

where ε ~ N(0, 1). The treatment variable A ∈ {0, 1} is generated from a logistic regression model with P(A = 1|Z) = expit(Z₁ + 0.7Z₁I(Z₂> 0) + Z₃) as P₁ and P(A = 1|Z) = expit(Z₁ + 3.5Z₁I(Z₂> 0) + Z₃) as P₂, respectively, where expit(x) = 1/(1 + e^−x). Notice that P₂ has a greater magnitude of non-linearity than P₁. Under the incorrect-model case, we fit the outcome variable Y with A and covariates X = (X₁, X₂, X₃) in place of Z = (Z₁, Z₂, Z₃) in both outcome model (m) and IPW model (π), where

In the omitted model case, we assume that Z₂ is wrongly omitted in the outcome model (m) and IPW model (π), resulting in fitting the outcome variable Y with treatment variable A and covariates (Z₁, Z₃) instead of Z = (Z₁, Z₂, Z₃) in outcome model (m) and treatment variable A with covariates (Z₁, Z₃) in IPW model (π).

In this setting, we seek to find which of the proposed estimators provide higher robustness to model disturbances under nine scenarios, where both or only one of outcome model and IPW model are accurately or incorrectly specified or omitted. For ease of presentation, we consider two outcome models, linear regression (Reg) and GSM, with different specifications in outcome model (m). We cover the IPW model (π) with a logistic model (GLM) and three machine learning methods, including tree (Tree), gradient boosting machine (GBM), and super learner (SL). The SL model consists of predictions from five base learners, (i) ipredbagg, (ii) ksvm, (iii) polymars, (iv) rpartPrune, and (v) KernelKnn, to provide a significant opportunity for improvement in predictive power. The base learners are fundamental factors of the SL ensemble model, which are combined to enhance the predictive performance synergetically. Each of these base learners are available in the following R packages;

(i) improved predictive models by indirect classification and bagging for classification, regression and survival problems as well as resampling based estimators of prediction error (ipred), (ii) kernel-based machine learning methods for classification, regression, clustering, novelty detection, quantile regression and dimensionality reduction (kernlab), (iii) polynomial spline routines for the polynomial spline fitting routines hazard regression, hazard estimation (polspline), (iv) recursive partitioning for classification, regression and survival trees (rpart), and (v) the extended simple k-nearest neighbors algorithm by incorporating numerous kernel functions and a variety of distance metrics (KernelKnn).

Tables 1 and 2 show the simulation results of IPWE, AIPWE and TMLE estimators in the correct, incorrect and omitted model specification of the outcome and IPW models, respectively, to estimate the average causal effect estimator (ACE). Simulation results in Tables 1 and 2 are based on the previously described model specifications with P₁ and P₂, respectively. The data sample size is fixed at N = 500. The true parameter τ is approximated with pseudo-random estimates with sample size N = 10⁶, and given by τ? = 1.996 for Table 1 and τ? = 1.999 for Table 2. Overall, the robust estimators, such as AIPWE and TMLE, perform well on finite sample sizes. It can be seen that AIPWE and TMLE have usually small biases when at least one model is correctly specified. Compared to naïve GLM, ensemble methods can significantly reduce both measures of bias and variance, contributing to MSE. When both models are misspecified or informative covariates are wrongly omitted in both models, biases cannot be completely removed, but still machine learning methods help effectively reduce bias and MSE. In Table 2, we found that AIPW and TMLE work stably well even when A is severely unbalanced. Note that AIPW and TMLE aided with SL algorithm show the best performance in nearly all cases.

Table 3 summarizes the results of CF, TMLE-SL and AIPWE-SL estimators. Overall, CF provides a great performance without requiring any complex model specifications. However, when both outcome and propensity models are correctly specified, CF seems to slightly underperform TMLE-SL and AIPWE-SL with respect to both bias and MSE. With small samples, CF has a relatively large MSE, but as the sample size gradually increases, the difference gradually decreases and in some case, CF shows a better MSE performance. However, we note that the simulation scenarios in Table 3 is not comprehensive and care should be taken for generalizing the previous findings.

Finally, Figure 1 displays eleven ACE estimators when both the outcome and propensity models are misspecified. See the note of Figure 1 for model description: for example, “model 4: TMLE-Reg-SL” denotes the composite TMLE method by implementing Reg for the outcome model and SL for the propensity model. Notice that CF (model 1) performs the best in this case, while two IPWEs (model 2 and 3) clearly underperform the other competitors due to model misspecification. Although AIPWE and TMLE work well without much difference, their performance depends on propensity model specifications rather than outcome models. For example, AIPWE-Reg-SL (model 5) outperforms AIPWE-Reg-GLM (model 7), while the difference between AIPWE-Reg-SL (model 5) and AIPWE-GSM-SL (model 8) is not much noticeable. On the other hand, by comparing IPWE-SL (model 2) and AIPWE-Reg-SL (model 5), we can see that doubly-robust estimation helps reduce both bias and variance even though model misspecification is present. We can conclude that the use of flexible machine learning methods, in particular for propensity models, is beneficial in reducing potential biases and achieving more robust results to model misspecification.

As an illustration, we consider the data from a percutaneous coronary intervention (PCI) observational study with 996 patients at Ohio Heart Health, Christ Hospital, Cincinnati, in 1997 (Li and Shen, 2020). In the study, researchers continued to support Lindner Center for at least six months using two treatment methods, abciximab (expensive high molecular weight IIb/IIIa cascade blocker) for 70% of patients and general treatment alone with initial PCI for 29.9% of patients. The dataset has ten variables from 996 patients without omission. We consider mortality and cardbill as response variables; mortality represents life years, having 0 or 11.6 depending on whether died within 6 months or not and cardbill represents the total cost of treating a patient’s heart within the first PCI 12 months at the Lindner Center. The abcix variable indicates whether the patient is treated abciximab or normal treatment alone. In addition, there are seven confounding variables; stent, whether to use a stent (anti-collapse device) or not; height, ranging from 108 cm to 196 cm; sex; diabetic state, whether to be diabetic or not; acutemi, whether the patient suffered from acute myocardial infarction within the last 7 days; ejecfrac, the ratio of left ventricle ejection from 0% to 90%; and ves1proc, the number of vessels involved in the initial PCI of patients from 0 to 5. Our main objective is to approximate the causal effects of the treatment variable abcix on mortality or cardbill. We apply several causal inference estimation methods to PCI data to determine whether the treatment is cost-effective after calibrating confounders and ensuring that there are more survivors of abcix treatment than people with other treatment. We estimate the bootstrap standard error by generating 10,000 bootstrap resamples. Although nearly 70% of patients are treated with abcix treatment, it is not clear that abcix treatment are reasonably inexpensive and actually makes people live more than others. Given only the numerical value of PCI data, the cardbill for abcix treatment is much more expensive than general treatment alone and survivors of each group are similar.

We begin by estimating propensity scores using logistic, Tree, GBM, and SL algorithms, including seven clinical predictors. Figure 2 displays the propensity score distributions for abcix treatment and normal treatment under four approaches, which shows the two distributions for each treatment group appear to well overlap in all cases. As in Li and Shen (2020), we can use Cochran’s 0.25 rule to measure how balanced the variables are for each group. The measurements for stent, acutemi and ves1proc are greater than 0.25 in standardized differences, but the balance in these groups can be achieved through the propensity score adjustment. Several causal effect estimates are considered to compare the causal inference effects between abcix and normal treatment. Table 4 summarizes the causal effects of treatment on 996 patients for two response variables, mortality and cardbill. In AIPWE-SL and TMLE-SL, we used a SL model for estimating the propensity scores by optimally mixing five machine base learners, (i) ipredbagg, (ii) ksvm, (iii) polymars, (iv) rpartPrune, and (v) KernelKnn.

The results in Table 4 imply similar conclusions among all causal estimators (except for the regression estimator): patients with abcix treatment generally survive more than those of other group and the cost of abcix is reasonably determined. Notice that Reg shows somewhat different results for cardbill, for which we suspect model misspecification. The results of IPWE are also less reliable due to large standard errors. The ACE estimates for mortality, when the general treatment (0) and abcix treatment (1) are compared after six months, is τ? = −0.051 with AIPWE-Reg-SL and τ? = −0.049 with TMLE-Reg-SL, and both are significantly lower than 0. The mortality rate of abcix treatment is also lower than that of general treatment alone, demonstrating that treatment of abcix is much more effective. In respect of the cost, the ACE estimate of AIPWE-Reg-SL is τ? = $983.534 and that of TMLE-Reg-SL is τ? = $968.201. The cost of abcix treatment would be a little higher than that of general treatment alone, but their difference is not statistically significant after considering confounding effects. In summary, our results suggest that patients treated with abcix lived longer while requiring a little higher treatment costs.

Doubly robustness (DR) is an important component in causal inference to correctly infer a target estimand after adjusting for potential confounders (Tsiatis, 2007). In observational studies, DR methods play an important role in estimating causal estimates if one of two working models (Q and π) is subject to model misspecification, but they may not work well if both components were not correctly specified. In this paper, we explore many DR inferential methods by replacing naïve linear regressions with flexible ensemble algorithms that are commonly used in causal analysis. For the outcome model, we mainly focused on GSM (Helwig, 2020), which is much more flexible than a simple regression model or generalized additive model (GAM). To handle the propensity model, we applied several ensemble algorithms, such as decision tree, GBM, and SL algorithms, and demonstrate their usefulness under model misspecification, but of course, other machine learning algorithms could be explored to effectively reduce potential bias. CF is a newly introduced causal method that is based on random forest. Unlike conventional causal inference methods, CF does not require a complex model selection procedure and is free of model misspecification. However, we found that when working models are close to true models, DR estimators, such as TMLE and AIPWE, often outperform CF.

In this paper, we tried to compare several methods as much as possible under different simulation configurations, but our results are limited in that we cannot cover all possible scenarios. Therefore, care should be taken not to overgeneralize our findings. Main implication from this paper is that DR estimators, such as AIPWE and TMLE, with support from machine learning algorithms work well even when both models are not perfect. Their performance seem quite similar and it is hard to say that one is superior to the other in our setting. CF works very reliably and shows relatively more consistenet results when there is a highly nonlinear pattern in the data. However, these causal estimators cannot completely remove biases especially when informative covariates are omitted in the model. Hence, investigators should more focus on study design to get more information about the target variable as much as possible. Finally, in our experiments, it is assumed that the outcome is completely observed. However, the outcome is often masked due to censoring in many observational studies. It is interesting to explore causal machine learning algorithms and measure their performance in survival analysis (Choi et al., 2021).