Estimating Average Treatment Effects Utilizing Fractional Imputation when Confounders are Subject to Missingness

1.1 Missing Data

Despite the best intentions of researchers, missing data is nearly impossible to avoid in real-world settings. Fortunately, as prevalent as missingness is, so too are methods with which to address missingness; however, the type of missingness matters when selecting a method.

Missing data occurs by way of one of three mechanisms in observational data. The first and simplest type is Missingness Completely at Random (MCAR) [6]. In this setting, whether or not an observation is missing is independent of both the observed and missing data. The second is Missingness at Random (MAR). Here, whether or not an observation is missing depends only on the observable data but not the missing data. Lastly, missing data can be Missing Not at Random (MNAR). In this setting, even after conditioning on all observed data, missingness will still depend on the missing portion of the data.

Confirming which missingness type is present for a given data set can not be validated from only the observed data. Little [7] proposed a chi-square based test capable of checking if the MCAR assumption was violated, and more recently Mohan and Pearl have been able to use m-graphs to refute instances where the MAR and MNAR assumptions were proposed but shown not to hold [8, 9]. Despite not being able to validate if a particular missingness model is true, it is still common for arguments towards one of these assumptions to be made, relying on the knowledge of subject matter experts for the data at hand. If the observed data is sufficient to explain the missingness, the MAR assumption is plausible, and it is this assumption we carry forward for the majority of our discussion. In Section 6 we will discuss the extension to MNAR.

1.2 Approaches to Address Missing Data

Many methods have been proposed to address missingness when assuming a MAR pattern. Two of the most common approaches in practice are complete case estimation (CC) and multiple imputation (MI). In particular, MI was favored by the National Research Council in 2010 as one of its preferred means of addressing missing data in clinical trials [10]. Under CC, all records with missing data are excluded, and treatment effects are estimated only on fully observed cases. CC can be biased under MAR; more importantly, if multiple variables have missing values, there will only be a small portion of complete cases in the data. By throwing out a large portion of the data, the effective sample size shrinks, inflating variances. Therefore, CC suffers from loss of efficiency by utilizing less of the observed data in the final analysis [11].

MI, on the other hand, is traditionally recommended for MAR in part due to its improved efficiency over CC and due to it being applicable to a wider range of missingness mechanisms. With MI, the full joint distribution of the data is estimated (either empirically or modeled based on distributional assumptions), and from this, a series of M new imputed data sets are drawn. In each imputed data set, MI fills in each missing value with an imputed value by sampling from the posterior predictive distribution of the missing value given the observed values. Then, full sample analyses can be applied to each imputed data sets, and these multiple results are summarized by an easy-to-implement combining rule for inference [12].

MI has produced valid frequentist inference in a wide range of applications [13]. At the same time, Rubin’s variance estimator for MI has been shown to not always be consistent [14, 15, 16, 17, 18, 19, 20, 21, 22]. For inference using MI consistent, imputation must be proper (see [12] for a precise definition). Practically speaking, for any proper imputations, under sufficiently regular models, Rubin’s combining rule provides a consistent estimator of the parameter of interest and a weekly unbiased estimator of its asymptotic variance. When the imputation model is correctly specified and the MAR assumption holds, Meng [23] showed that a sufficient condition for imputation to be proper is that the imputation model and the analysis model are congenial. For example, the imputation model is correctly specified and the analysis is efficient under the same model. Congeniality is sometimes more elusive than it would appear. Even when the imputation model is correctly specified, Yang and Kim [22] showed that MI is not necessarily congenial for method of moments estimation. Therefore, some common statistical procedures may be incompatible with MI. From a causal inference perspective, this poses a problem as the validity of Rubin’s variance estimator has not been fully explored for many full sample estimation methods used widely in causal inference. Certain otherwise unbiased and consistent full sample causal inference methods (outcome regression, weighting, matching, etc.) may lose these properties when applied in conjunction with MI and MI-produced data sets. Many of the most common estimates for average treatment effects are based on method of moments estimators and are thusly susceptible to inaccurate variance estimates when using Rubin’s variance estimator for MI. For researchers desiring to make causal claims when utilizing MI, it is imperative for the variance properties of their estimators to, therefore, be either validated or an alternative method must be proposed.

1.3 Fractional Imputation as an Alternative

As an alternative to CC and MI, there are likelihood-based methods that can be applied. When using these methods, the key insight is that under fully observed confounders, the full sample estimators are obtained by solving estimating equations. In the presence of partially observed confounders, the corresponding estimators can be obtained by solving conditional estimating equations which integrate out the missing confounders given the observed data. There are two difficulties in this approach. First, it requires consistent estimators in the conditional distribution of the missing confounders given the observed data, such as MLE. In the presence of missing values, an EM algorithm is typically used. Second, numerical integration is needed. Integration approximated by imputation was considered by many authors, such as Monte Carlo EM methods [24]. For Monte Carlo EM algorithms, in each E-step, the imputed values are regenerated, and thus the computation can be quite heavy. Also, the convergence of Monte Carlo sequence of the estimators is not guaranteed for fixed Monte Carlo sample size [25].

In practice, EM algorithms may not be feasible when the conditional expectation in the E-step is not available in a closed form. Instead, fractional imputation (FI) has been proposed to serve as a computational tool for implementing the expectation step (E-step) in the EM algorithm [24, 26, 27], which simplifies computation by drawing on importance sampling to obtain the fractional weights and reducing the iterative computation burden over other simulation methods such as Markov Chain Monte Carlo. See [28, 29, 30, 31] for applications of FI outside of the causal inference context. The main idea in FI is to produce a complete data set by imputation where each imputed value is associated with a fractional weight, by which the observed likelihood can be approximated by the weighted average of the imputed data likelihood. The resulting estimator approximates the maximum likelihood estimator.

For illustrative purposes, we can compare the format of the imputed data via FI to the more commonly encountered imputed data via MI. Suppose we have a data set where some of the records are subject to missingness. MI will attempt to address this missingness by creating M data sets where each imputed value is imputed exactly once in each of the M data sets as depicted in Figure 1. In this image X₁ is a fully observed covariate, and A and Y are fully observed treatment indicator and response variables respectively. X₂ is a covariate subject to missingness. R₁ and R₂ are missingness indicators for X₁ and X₂ where R_ij = 0 indicates X_j is missing for record i. Records with missingness are indicated in red on the left. The completed data sets (with X₂ imputed) are on the right with the records now including imputed values for X₂ highlighted in green.

Figure 1

Imputed Data via MI

FI instead seeks to generate a single imputed data set, but one where each imputed record also includes a fractional weight. That fractional weight is indicative of how likely the imputed data is to occur under the distribution of the completed data set. Figure 2 shares all the same features as Figure 1 except now includes an additional column in the imputed data for the fractional weight (ω) that will be incorporated into analyses.

Figure 2

Imputed Data via FI

At its most basic level, the difference in imputation approach is analogous to a debate between a wide vs a tall data set. It could even be argued that when recombining data sets imputed via MI an implicit "weight" of 1/M is assigned to every record. This implicit 1/M weight acts similarly to the fractional weights ω generated via FI. The important distinction between MI and FI lies in how the final "weights" assigned to each imputed record are produced. The weight generation process for FI will be discussed in more depth in Section 3. For the time being, we will assume these weights are generated appropriately and can be incorporated easily into further analyses.

FI is sufficiently versatile that, like MI, it can be deployed along-side both continuous and categorical variables. Under the same common regularity conditions as seen often with MI, we show that the FI estimators are consistent and asymptotically normal. The remainder of this paper focuses on expanding FI into the causal literature by developing an FI-based method for estimating causal treatment effects. Once developed we will validate the method relative to existing causal inference methods. Specifically, we investigate the comparative performance of FI vs MI and CC for estimating causal treatment effects when confounders are subject to missingness.

In summary, the proposed FI framework achieves desirable properties for causal inference:

1. the same fractionally imputed data set allows for applying general full-sample estimators (which solve certain estimating equations) of the average treatment effect, including regression estimators, inverse probability of weighting estimators, and augmented weighting estimators;

2. the FI estimators are asymptotically linear and therefore allow resampling methods for variance estimation and inference, which is simple to implement in practice;

3. the unified FI inference has theoretical guarantees and offers a solution to the uncongeniality issue of MI;

4. and lastly, FI is not only statistically efficient but also computationally efficient compared to MI, as demonstrated via simulation and real-world application.

The rest of the article is organized as follows. We begin in Section 2 with a description of notation and assumptions. In Section 3 we fully outline the FI process as well as derive the resulting variance estimator for the treatment effect estimator. We implement a simulation study in Section 4 comparing the accuracy and precision of treatment effect estimators when the missingness is addressed by CC, MI, and FI. Section 5 provides a demonstration of FI’s utilization with a real-world health data set. Finally, in Section 6, we end with a discussion of the results and of the implication they have on current and future causal work.

2.1 Treatment Effect Estimation Notation

Following Neyman [32] and Rubin [33] we use the potential outcomes framework. Treatment is denoted by A_i ∈ {0, 1}, where 0 and 1 are labels for control and treatment respectively. For each subject i, define a pair of potential outcomes {Y_i(1), Y_i(0)} which represent the outcomes if the subject was treated, Y_i(1), and if he or she was not, Y_i(0). Implicit in this notation, we make the stable unit treatment value assumption [34]. The observed outcome for subject i is then Y_i = Y_i(A_i). Let X_i be the vector of confounders for subject i. We assume that {X_i , A_i, Y_i(1), Yi(0)}i=1nare independent draws from the distribution {X, A, Y(1), Y(0)}, and therefore, {(Xi,Ai,Yi)}i=1nare independent and identically distributed. The conditional treatment effect is τ(X) = E{Y(1) − Y(0) | X}, and the average treatment effect is τ₀ = E{τ(X)}. The average treatment effect cannot be estimated without further assumptions, because for each subject only one potential outcome is observed. The common assumptions for identifying the average treatment effect [1] are as follows:

Assumption 1 (Ignorability). Y(a)⊥A | X for a = 0,1, where ⊥ means "is conditionally independent of".

Assumption 2 (Sufficient overlap). With probability 1, 0 < c₁ ≤ e(X) ≤ c₂ < 1, where e(X) = pr(A = 1 | X) is the propensity score.

Under Assumption 1, adjusting for covariates X creates a randomization-like scenario and removes confounding biases brought on by treatment selection (i.e., for each level of X, the treatment assignment is as good as randomization). However, in practice, there are often many variables in X, some of which are continuous; therefore, directly conditioning on each level of X is difficult. Alternatively, the propensity score has been proposed as a one-dimensional summary of X [1]. The central role of the propensity score lies in the fact that Assumption 1 implies Y(a)⊥A | e(X) for a = 0,1. Therefore, adjusting for the propensity score alone can remove confounding biases.

2.2 Treatment Effect Estimation Under Fully Observed Data

Reliance on the propensity score comes about naturally when we decompose the joint density of X, A, and Y into three particular components. Specifically

Based on this decomposition, a number of propensity score-based estimators have been proposed for estimating the treatment effect including propensity score matching, subclassification, or weighting. See [35] for a textbook discussion. If we limit the class of propensity estimators to only parametric estimates (of the form e(X | θ) where θ is the vector of parameters used to estimate the propensity score), the two most common estimates of τ from this class are that of Inverse Propensity Weighting (IPW) and Augmented IPW (AIPW). Both are illustrated here as examples. In both examples, let e(X | θ^)be the estimated propensity score where θ has been estimated by some consistent estimator θ^.In practice, e(X | θ^)is typically a logistic regression model.

2.3 Missing Data Notation

Whereas all of the above discussion holds under fully observed responses and confounders, in this article we consider the case where X contains missing values. To that end, let R be a collection of indicator variables R = (R₁, ..., R_p) corresponding to (X₁, ..., X_p) where R_ij = 1 indicates that X_j is observed for subject i and R_ij = 0 indicates X_j is not observed for subject i. Let R denote the collection of all possible missingness patterns. Let X_obs,i, the observed part of covariates for individual i, consist of X_ij with R_ij = 1. Similarly let X_mis,i, the missing part of covariates for individual i, consist of X_ij with R_ij = 0. Note that the dimension of X_obs,i and X_mis,i will vary across individuals. For notational simplicity, we suppress the subscript i for subject.

If we return to equation (1) and incorporate the parameterizations laid out in examples 1 and 2, the decomposition of the joint distribution can be naturally extended to account for missingness. To do so, we rewrite the joint distribution as:

where α is the collection of parameters used to estimate the missing portion of X given the observed portion of X, and ρ is the collection of parameters used in describing the missingness mechanism.

Finally, in this article we are only interested in the case where X_mis follows a MAR pattern, which leads us to our final assumption which completes the basis for Theorem 1 (the proof for which is made available in the appendix).

Assumption 3 (Missing at random). R⊥Xmis   ∣Xobs   ,A,Y

Here we adopt the notion of Rubin’s MAR in the sense that MAR may hold for different variables depending on the missingness pattern. The scientific justification of this assumption may be difficult; however, theoretically it is the weakest condition under which the missingness process can be ignored [6]. Alternatively, one can consider the variable-based taxonomy of MAR [9], where X_mis represents variables that are subject to missingness and X_obs represents variables that are fully observed. Our method development is similar under this notion of MAR but using the new definition of X_mis and X_obs in equation (4).

2.4 Treatment Effect Estimation Under MAR

In the presence of missingness, we can let U¯(τ;Xobs,A,Y∣η^)=E{U(τ;X,A,Y∣η^)∣Xobs,A,Y}.We term U¯(τ;Xobs,A,Y∣η^)the "mean estimating function" of τ₀ given the observed data. Note that in defining U¯(τ;Xobs,A,Y∣η^),η^must now also expand to include any new parameters α^utilized in estimating E(X_mis | X_obs, A, Y; α). From Theorem 1, under Assumptions 1–3, U¯(τ;Xobs,A,Y∣η^)is an unbiased estimating function of τ₀. Therefore, a consistent estimator of τ₀ can be obtained by solving PnU¯(τ;Xobs,A,Y∣η^)=0.

3.1 Implementing Fractional Imputation

Earlier in this paper, Figure 2 depicted what a completed data set might look like after using FI to impute missing data during the design stage. Besides the shape of the data, the addition of an explicit weight to each record of the imputed data set stands out from more commonly encountered MI data sets. To generate fractional weights, FI deploys a three-step process. First, missing values are imputed from some proposed distribution, then fractional weights are updated, then model parameters for the full joint distribution are re-estimated (now under updated weights). The re-weighting and model updating steps cycle until convergence.

In the first step, sometimes referred to as the imputation step or the I-step [27, 40], every missing value is imputed M times by means of some proposal function h(X_mis | X_obs). This generates M new values Xmis,i⋆(j)for each partially observed record. The choice of h() is arbitrary and left to the imputer, but a convenient choice is f (X_mis | X_obs; α) to align with the decomposition from equation (5). This choice would necessitate being able to provide or estimate a value for α (for instance α = α^0where α^0is the MLE for α calculated only using complete cases). At this first step, every record where a value was imputed will have a fractional weight of ωij⋆=1/M.Note, at the conclusion of every step in FI, fractional weights for each individual, observation, etc. i are held to the condition∑j=1Mωij⋆=1.

In the second step, the weighting step or W-step, the fractional weights are updated proportional to the likelihood of the imputed value under the full joint distribution, divided by the likelihood of the imputed value under the proposal distribution h() from the I-step. If choosing the decomposition from equation (5) that would appear as

given the current parameter estimates for η^(t)=(α^(t),β^(t),θ^(t)).

In the third step, the maximization step or M-step, the parameter values used to estimate the full joint distribution are updated given the new values of ωij⋆ from η^(t) to η^(t+1).The W-step and M-step iterate, setting t=t+1 after each M-step, until the model parameters converge. The resulting data set with the final fractional weights included can then be passed on to the analysis stage as if it were a complete data set (similar to MI passing on M data sets). An example generating a complete data set via FI is included in the appendix.

3.2 Characterizing Fractional Imputation for Estimating Treatment Effects

For illustration, consider the case where X contains only two variables, X = (X₁, X₂),where X₁ is fully observed and X₂ is subject to missingness. Let R₂ be the response indicator of X₂. From Examples 3 and 4, we can obtain an estimator of τ₀ by solving the mean estimating equation

where U(τ;Xi,Ai,Yi∣η^)denotes either UIPW(τ;Xi,Ai,Yi∣η^) or UAIPW(τ;Xi,Ai,Yi∣η^).

In (6), the conditional expectation E{U(τ;Xi,Ai,Yi∣η^)∣Xobs,i,Ai,Yi}is often difficult to obtain. The basic idea of FI is to overcome this difficulty by creating a weighted set {(ωij⋆,Xi⋆(j),Ai,Yi):j=1,…,M}such that E{U(τ;Xi,Ai,Yi∣η^)∣Xobs,i,Ai,Yi}can be approximated by∑i=1n∑j=1Mωij⋆U(τ;Xi⋆(j),Ai,Yi∣η^).

3.3 Asymptotic Results

Because τ^is obtained through the method of estimating equations, we establish the asymptotic properties of τ^,in a manner similar to [20].