All models are wrong, but which are useful? Comparing parametric and nonparametric estimation of causal effects in finite samples

1 Introduction

In the past two decades, work emerging from the nonparametric statistics literature has allowed for the incorporation of machine learning algorithms in estimation while, in many cases, preserving theoretically valid statistical inference [e.g., 1–4]. A growing faction, though still a subset, of statisticians, epidemiologists, econometricians, and other applied scientists now take as a given that these nonparametric estimators, if appropriately applied, are typically superior to parametric estimators, particularly in terms of reducing or eliminating model misspecification bias, which can be substantial when using an overly simplistic parametric model in complex settings [e.g., 4–14]. (We note that we use the term “nonparametric estimator” to mean an estimator where each component of the causal model used in estimation [also called a “nuisance parameter,” e.g., an outcome regression] is estimated nonparametrically.) Minimizing all sources of bias is essential, especially because when bias remains on the same order as standard error, the probability of a corresponding hypothesis test rejecting the null will tend to one even when no effect is present.

Even so, an attitude of caution and skepticism remains prevalent, with many favoring simpler, parametric model-based approaches over approaches that are deemed by some to be unnecessarily complex [e.g., 15]. For example, Little (2013) observes that “Much modern nonparametric and semiparametric statistical theory lives in the ‘land of asymptotia’ [15].” This is accurate, and few have contributed evidence of the extent to which the assumption of nonparametric estimator superiority is borne out in real-world, finite-sample analyses, leaving the debate as to whether nonparametric methods confer any meaningful practical advantage in typical real-world applications largely unresolved [16].

Both nonparametric and parametric estimation methods generally involve using theory and subject-matter knowledge to inform the underlying causal model/graph [17] and the variables input into the model [18]. However, whereas parametric estimation methods would proceed by specifying parametric models (e.g., linear regression) for each nuisance parameter, specifying its functional form, interactions, higher-order variable forms, etc., nonparametric methods would typically involve using data-adaptive methods, such as machine learning algorithms, in flexibly fitting nuisance parameters.

Statistical theory informs when we would expect parametric versus nonparametric estimators to perform better in terms of bias. The bias of parametric estimators is a function of how close the parametric regression model is to the truth. For a simple parametric plug-in estimator (discussed further in Section 3 and Section S1 in the Supplement), bias is a function of how close the parametric outcome model is to the truth. For parametric estimators based on the propensity score, bias is a function of how close the treatment model is to the truth. For example, if the response and/or treatment assignment surfaces are nonlinear, then a parametric linear regression-based estimator that does not incorporate such nonlinearities would be expected to have appreciable bias [19]. The bias of nonparametric estimators also depends on how close the nonparametric regression models are to the truth, but the reliance on correct specification of the regression models is relaxed in at least one or two respects, depending on the estimator. First, the bias of nonparametric estimators that are also semi-parametric efficient estimators (e.g., the augmented inverse probability of treatment weighted [A-IPTW] estimator, the one-step estimator, double/de-biased machine learning estimator, targeted minimum loss-based estimator [TMLE]) is an integral of a product of two regression errors – e.g., one for the outcome regression and one for the treatment regression, a property called “double robustness” [20–22]. If the two errors are small, this product will generally be smaller than the single error involved in the definition of the parametric estimator bias. We discuss this further in Section S1 in the Supplement and TMLE in more detail in Section 3. Second, nonparametric estimators use machine learning algorithms in fitting nuisance parameters, which reduces reliance on correct parametric model specification, and therefore, can improve model fit and reduce bias even further. Consequently, because correctly specifying parametric models a-priori seems unlikely in the absence of clear mechanistic understanding, based on this asymptotic theory, nonparametric estimators are expected to have less bias, on average, than parametric estimators in very large samples when an accurate model is important to obtain an unbiased estimate, e.g., in the presence of significant, complex confounding or censoring [19]. However, asymptotic approximations may be poor in small or even moderately sized real-world samples that are afflicted by the curse of dimensionality, resulting in bias and under-coverage of confidence intervals [23]. The sample size at which asymptotic approximations reflect reality is generally problem dependent and unknown.

Previous simulation studies comparing performance between parametric and nonparametric estimation have generally shown that nonparametric estimators incorporating data-adaptive, machine learning algorithms outperform parametric approaches, especially under the data structures enumerated earlier [19,24–26]. However, common critiques of simulation studies are that they (a) are typically designed to illustrate a difference in performance (or lack thereof), and thereby reflect hand-picked settings that are favorable to the method of interest; (b) are oversimplified “stylized models of reality”; and (c) may hold little value if the goal is to make a general statement about the degree to which the choice of nonparametric versus parametric estimation matters [27–31]. A similar criticism regarding the lack of generalizability can apply to the few previous real data analysis examples comparing nonparametric versus parametric estimators, sometimes finding meaningful differences and other times, not [14,32].

Responding to the critique that simulation-based evaluations of estimator performance have been conducted on oversimplified data too far removed from the complexity of real-world datasets, several groups recently proposed more comprehensive, yet tailored, evaluations of estimator performance by simulating datasets that closely mimic an observed dataset in all its complexity (similar to the idea of “plasmode simulation” [33]) [e.g., 27,28,34,35]. This work is premised on the belief that the optimal estimator will differ for different data scenarios; hence, the focus is on generating simulated data to match a particular observed dataset. However, by using the same data to both select an estimator and evaluate this estimator, these approaches are vulnerable to problems with post-selection inference. Given sufficiently large sample sizes, this problem can be circumvented via sample splitting [36]; however, practitioners can seldom afford to discard large amounts of data to perform estimator selection.

The current debate about how much the choice of a nonparametric estimator versus a parametric estimator may matter in the real world asks a more general question than the current literature can answer. We propose an approach to help answer such a question: what we call a large-scale Monte-Carlo simulation that instead of considering just a few data-generating mechanisms, summarizes estimator performance across thousands of data-generating mechanisms within a user-defined data setting. Unlike the plasmode simulation approaches referenced earlier, ours does not depend on the observed data. We apply this approach, comparing several parametric and nonparametric estimators to demonstrate how one can generate broader evidence of the degree to which choosing among these estimators can impact bias, confidence interval coverage, and mean squared error (MSE) in finite sample causal effect estimates.

Specifically, in our application of the proposed approach, we consider data-generating distributions drawn from the following setting: a binary treatment, a binary outcome, five binary covariates, and one multi-valued covariate. We compare performance in estimating the average treatment effect ( ATE , formally defined in Section 2). We consider: (1) a parametric “plug-in” estimator (also called a “g-computation estimator”), (2) a propensity score estimator that uses a parametric treatment model, (3) a nonparametric “plug-in” estimator, and (4) a nonparametric efficient estimator.

This article is organized as follows. In Section 2, we introduce notation and the causal estimand we consider. In Section 3, we describe our proposed large-scale Monte-Carlo simulation for comparing parametric and nonparametric estimators. Also in this section, we provide the specific parameters that define the space of data-generating mechanisms we consider in applying the proposed approach, the parametric and nonparametric estimators we compare, and how we quantify and summarize estimator performance in finite samples. In Section 4, we provide and discuss results. Section 5 concludes.

2 Notation and estimands

For simplicity, in this section, we focus on one component of the ATE contrast – the expected counterfactual outcome had treatment been set to some value, t , possibly contrary to fact, denoted E ( Y t ) , where T denotes a binary treatment variable, and Y denotes the outcome variable. The ATE would then be denoted E ( Y 1 ) − E ( Y 0 ) . We assume that the exchangeability assumption holds, Y t ⊥ ⊥ T ∣ X for t ∈ { 0 , 1 } , where X denotes the confounding variables. This assumption is necessary for the causal parameter, E ( Y t ) , to be identified from the observed data, Z = ( X , T , Y ) , by the statistical parameter θ ( P ) = E [ E ( Y ∣ T = t , X ) ] , where θ denotes the statistical parameter, which is a function of the data distribution, P . We let P be an element of the nonparametric statistical model defined as the set of all possible observed data distributions, noting that nonparametric statistical models can be correct. We also assume positivity, P ( T = t ∣ x ) > 0 for all x in the support of X , and for each treatment value in the causal parameter of interest. This is necessary for the statistical parameter to be well-defined.

We note that θ ( P ) constitutes the building block for the estimation of common marginal causal effects: the ATE , identified by E { m ( 1 , X ) − m ( 0 , X ) } ); the relative risk (RR), identified by E { m ( 1 , X ) } E { m ( 0 , X ) } ; and the odds ratio (OR), identified by E { m ( 1 , X ) } ∕ E { 1 − m ( 1 , X ) } E { m ( 0 , X ) } ∕ E { 1 − m ( 0 , X ) } , where m ( t , x ) represents the true conditional density function (i.e., regression) E ( Y ∣ T = t , X = x ) , and where we use subscripts with m ( t , x ) to denote the estimates of that conditional density based on a particular statistical model. For example, we use m β ˆ ( t , x ) to denote an estimate of that distribution that is a function of fitted parametric regression coefficients, β ˆ . Similarly, we use subscripts with θ to denote the estimates of the statistical parameter with the subscript denoting the particular estimator used. For example, we use θ T M L E to denote a TMLE estimate of the parameter θ ( P ) .

When the exposure takes values in a discrete set, like the binary exposure we consider here, the common marginal causal effects of the ATE , RR , and OR can be nonparametrically estimated, yet with the same asymptotic properties as if they were parametrically estimated, as long as the dimension of X is fixed [4]. Taking the ATE as an example, this means that it will be possible to find an estimator ATE ^ such that as the sample size n grows, n ( ATE ^ − ATE ) converges to a random variable that is normally distributed with mean zero and variance equal to the nonparametric efficiency bound. The efficiency bound is the smallest possible variance attainable by a regular estimator [1]. This means that the normal distribution can be used to approximate the sampling distribution of the estimator ATE ^ for finite sample sizes, which allows us to construct approximately correct confidence intervals and hypothesis tests.

In the next sections, we develop a simulation-based approach to evaluate if and when nonparametric estimators of θ ( P ) can be expected to outperform parametric estimators across a range of finite samples under minimal but reasonable assumptions on the nature of the true data-generating mechanism.

3 A large-scale Monte-Carlo simulation approach to systematically evaluate the finite sample performance of estimators

We propose the following general method to systematically evaluate and compare the performance of a set of candidate estimators in finite samples. Our method, which draws on the nonparametric Bayesian literature, compares the performance across a large number of data distributions within a setting defined based on the following: the number and type of variables, range of confounding bias, treatment effect heterogeneity, and degree of nonlinearity.

• For each j ∈ { 1 , … , J } , we draw a probability distribution P j using a minimally informative prior across the space of distributions (i.e., nonparametric model), ℳ . The minimally informative distribution on ℳ reflects the fact that investigators often have little knowledge of the data-generating mechanism before seeing the data. We describe this part of the procedure in more detail in Section 3.1.

• Then, for each of the J data-generating mechanisms, we generate S simulated datasets. So, for each j and s we have a dataset D j , s of some finite sample size n from P j .

• Then, for each of the L estimators being considered, we compute the parameter estimate θ ˆ l , j , s and the standard error estimate σ ˆ l , j , s ^[1]. Let θ ˆ 1 , j , s , … , θ ˆ L , j , s denote the parameter estimates and σ ˆ 1 , j , s , … , σ ˆ L , j , s denote the standard error estimates.

• We next compute the performance metrics based on sample averages across the S simulated draws. For distribution P j and sample size n , the Monte-Carlo bias, confidence interval coverage, and MSE of a given estimator θ ˆ l are

  B l ( n ) ( P j ) = 1 S ∑ s = 1 S θ ˆ l , j , s − θ ( P j ) , CICov l ( n ) ( P j ) = 1 S ∑ s = 1 S 1 θ ˆ l , j , s − c σ ˆ l , j , s n < θ ( P j ) < θ ˆ l , j , s + c σ ˆ l , j , s n , MSE l ( n ) ( P j ) = 1 S ∑ s = 1 S { θ ˆ l , j , s − θ ( P j ) } 2 ,

  where c denotes the 97.5th percentile of the standard normal distribution.

• Finally, we summarize the performance across the J data distributions with a probability distribution of each performance metric (e.g., the bias). Specifically, the bias for estimator l at sample size n is summarized by:

  P ( ∣ B l ( n ) ∣ > b ) = 1 J ∑ j = 1 J 1 { ∣ B l ( n ) ( P j ) ∣ > b } .

  This function can be seen as an approximation of the so-called survival or reliability function, where the former name is common in the statistics literature and the latter in the engineering literature. We use the name reliability function, because it is closer to the intended interpretation for each estimation procedure. Because we consider the sample space of probability distributions as all possible data-generating mechanisms (i.e., phenomena under study), the probability distribution can be heuristically interpreted as probabilities across all possible studies.

In contrast to standard simulation studies, which often focus on a few data-generating mechanisms, this large-scale Monte-Carlo simulation will provide evidence of the behavior of the estimators across many data-generating mechanisms and, consequently, will be informative for assessing how the estimators perform across a broader range of problems that may be encountered within the defined data setting, ℳ .

4 Results

We stratify results by each estimator, the three sample sizes, and by complexity of the data-generating mechanism. To vary complexity, we examine (a) no interactions between variables and no treatment effect heterogeneity, (b) treatment effect heterogeneity but no interactions between variables, (c) up to third-order interactions between variables and no treatment effect heterogeneity, and (d) up to third-order interactions and treatment effect heterogeneity. We also stratify results by degree of practical violations of the positivity assumption.

We see in Figures 1 and 2 that at least one, and often, both of the nonparametric estimators, BART and CV-TMLE, perform better than the parametric estimators in terms of both bias and MSE. The one exception is for the smallest sample size of N = 100 and the simplest two data-generating mechanisms that do not have any interactions between variables, the parametric plug-in estimators perform similar to TMLE in terms of bias (Figure 1(a) and (b)), and BART performs relatively worse. In terms of MSE, however, at least one of the nonparametric estimators performs best across all settings (Figure 2), with BART dominating in all cases when N = 100 and performing at least close to the best in all other scenarios. We also include a figure summarizing variance in this simulation study in Figure S3, and note that in general, BART has the smallest variance and CV-TMLE has the largest.

Figure 1

Reliability function for the absolute bias of each of the four estimators (parametric IPTW, parametric plug-in [g-computation: (i) main terms only and (ii) including treatment–covariate interactions (TCI)], nonparametric plug-in [BART], and nonparametric efficient estimator [CV-TMLE]) when considering data-generating mechanisms characterized by one binary treatment, one binary outcome, five binary covariates, one numeric covariate, and various levels of model complexity. (a) No treatment effect heterogeneity, no interactions, (b) treatment effect heterogeneity, no interactions, (c) no treatment effect heterogeneity, up to three-way interactions, and (d) treatment effect heterogeneity, up to three-way interactions.

Figure 2

Reliability function for the absolute MSE of each of the four estimators (parametric IPTW, parametric plug-in [g-computation: (i) main terms only and (ii) including treatment–covariate interactions (TCI)], nonparametric plug-in [BART], and nonparametric efficient estimator [CV-TMLE]) when considering data-generating mechanisms characterized by one binary treatment, one binary outcome, five binary covariates, one numeric covariate, and various levels of model complexity: (a) no treatment effect heterogeneity, no interactions, (b) treatment effect heterogeneity, no interactions, (c) no treatment effect heterogeneity, up to three-way interactions, (d) treatment effect heterogeneity, up to three-way interactions.

The largest separation between the parametric versus nonparametric reliability curves in Figures 1 and 2 occurs for sample sizes N ∈ { 500 , 1,000 } and in the more complex settings of variable interactions and possibly treatment effect heterogeneity. For example, in the most complex setting with N = 1,000 (Figure 1(d)), the nonparametric estimators result in absolute bias > 0.1 in 0% of data distributions, but the parametric estimators result in absolute bias > 0.1 in 11.1 and 13.3% of data distributions for the parametric main-terms plug-in estimator and IPTW estimator, respectively. This result was anticipated, because in larger sample sizes, the asymptotic advantages of nonparametric estimators to model complexity and nonlinearities may take hold.

We also stratify by degree of practical positivity violations in Figures S4 and S5. We categorize practical positivity violations as follows: “minimal violations” are those in which the minimum (across all iterations) predicted conditional probabilities of treatment ≥ 0.02 ; “moderate violations” are those in which the minimum predicted conditional probabilities of treatment < 0.02 but ≥ 0.01 ; and “severe violations” are those in which the minimum predicted conditional probabilities of treatment < 0.01 . Although estimator performance degrades slightly in the presence of severe positivity violations, particularly for the IPTW estimator, the relative performance of the estimators remains the same.

In terms of 95% confidence interval (CI) coverage (Table 2), we see that all estimators perform well for sample size N = 1,000 in the simplest scenarios without interactions, with slight over-coverage by IPTW and BART. However, with more complexity in terms of interactions, coverage suffers for the parametric estimators, particularly with increasing sample size; in contrast, the nonparametric BART and CV-TMLE estimators continue to attain at least the nominal coverage rate. We show in Table S1 the importance of using cross-fitting in attaining nominal coverage for TMLE when using data-adaptive regression estimators for asymptotic normality [4,57,58]. Interestingly, BART tends to over-estimate uncertainty leading to coverage above the 95% level. This is not generally expected to hold for BART, though potentially this could be addressed using improvements to the BART prior, such as the SoftBART prior [59], which is better suited to many realistic data-generating processes and has better asymptotic properties.

5 Discussion

We proposed a large-scale Monte-Carlo simulation method to bring more comprehensive evidence to the debate [16] as to whether nonparametric estimators that use data-adaptive machine learning algorithms in fitting nuisance parameters confer any meaningful advantage over simpler parametric methods in real-world finite sample analyses. Previously, others contributed finite sample evidence in favor of one estimator or class of estimators over another by conducting simulation studies across a smaller, more bespoke set of data-generating processes [19,24,25,27,28], which may not be representative of performance in general. Consequently, our proposed large-scale approach greatly expands the number of data-generating mechanisms considered from a small few to thousands within a data setting that a user can define based on some minimal knowledge about observed variables in their data (but without the need to first look at their data), thereby likely resulting in more generalizable – and thus, informative – evidence of estimator performance in finite samples where asymptotic properties learned from theoretical results may not provide good approximations.

We applied our proposed approach to compare the performance of two nonparametric estimators, BART and CV-TMLE, to two parametric estimators, an outcome regression-based plug-in estimator and IPTW, in finite samples of sizes N = 100 , 500, and 1,000 and across different degrees of model complexity. Corroborating findings from previous, smaller-scale simulation studies [19,24–26], in our application of the proposed approach, we found that even in small samples, nonparametric estimators nearly always outperform the parametric estimators in terms of bias and MSE. However, the advantage of nonparametric estimation attenuated with decreasing sample size and decreasing complexity, and in the simplest data-generating mechanisms and samples of N = 100 , the parametric plug-in estimator performed similar to the nonparametric efficient estimator in terms of bias and between the nonparametric estimators in terms of MSE.

Even though our results were learned from 1,773,000 data sets across 3,000 data-generating mechanisms, the space of data-generating mechanisms we considered was nonetheless limited. In particular, we only considered settings with a binary treatment, a binary outcome, and six confounding variables, only one of which was multi-valued. It is certainly possible that our conclusions would differ for more complex settings. In future work, we will develop a software tool for running simulations with user-specified outcome and covariate types and covariate dimensions. This way, users would be able to compare estimator performances over a space of data-generating mechanisms that are likely to contain the probability distribution corresponding to their real-world data sampling setting or a close approximation thereof. Such future work could be further extended to accommodate longitudinal data as well as common data complexities such as missing data, measurement error, and selection bias.

Finally, although we have shown that the choice of nonparametric vs parametric estimator may matter – in some cases more than others – all estimators are limited by the data input. Unmeasured variables and variables measured with error are significant and near-ubiquitous limitations that can thwart accurate causal effect estimation and inference. A relatively recent high-profile and high-stakes example involved data errors leading to inaccurate algorithmic predictions in the criminal justice system that resulted in unintended parole denials [60,61]. We can work to improve the estimation step of answering research questions, but the accuracy of our answers will be limited by the weakest link, highlighting the importance of theory, subject matter knowledge, identification, and data quality, in addition to estimation.

In this article, we have focused on strengthening the estimation link. Our results show that in the large space of settings we have considered, this can be accomplished by employing nonparametric estimators grounded in asymptotic theory to substantially reduce bias in large-sample settings with interactions and nonlinearities while compromising very little in terms of performance even in simple, small-sample settings.