Many methods have been developed to estimate the causal effect of an intervention, such as the administration of a treatment, on an outcome such as survival, from observational data, i.e., data that is potentially confounded by selection bias due to the absence of randomization. Classical ones include matching (Iacus et al., 2012), inverse propensity weighting (IPW, Horvitz & Thompson, 1952; Rosenbaum & Rubin, 1983) and doubly robust methods (Robins et al., 1994; Chernozhukov et al., 2018; Wager & Athey, 2018; Athey et al., 2019). More recent proposals use deep learning methods that ensure balance of the population at the level of representation (Johansson et al., 2016; Shalit et al., 2017), infer the joint distribution of latent and observed confounders, the treatment and the outcome (Louizos et al., 2017) or predict the counterfactuals with GANs (Yoon et al., 2018). For a detailed review of existing literature on treatment effect estimation we refer to Imbens (2004), Lunceford & Davidian (2004) and Guo et al. (2019).

However, state-of-the-art methods still suffer from important shortcomings. In particular, they seldom consider the possibility that covariates have missing values, which is ubiquitous in many real-world situations (Josse & Reiter, 2018) and has been widely discussed in different contexts (Mayer et al., 2019a; van Buuren, 2018; Little & Rubin, 2002). Although this question of missing attributes in the context of treatment effect estimation has been raised early in the development of causal inference (Rosenbaum & Rubin, 1984), there is still a lack of effective and consistent solutions addressing this problem, with a few notable exceptions such as Mattei & Mealli (2009); Seaman & White (2014); Yang et al. (2019); Kallus et al. (2018) which mainly focus on inverse propensity weighting (IPW) methods and Kuroki & Pearl (2014) who discuss identifiability of causal effects under measurement error or unobserved confounders. Recently, Mayer et al. (2019b), in addition to suggesting doubly robust estimators with missing data, classified the existing approaches into two families: the ones that adapt the causal inference assumptions to the missing values setting (D’Agostino Jr & Rubin, 2000; Blake et al., 2019) and the ones (Mattei & Mealli, 2009; Seaman & White, 2014; Kallus et al., 2018) that consider the classical machinery and missingness mechanisms assumptions (Little & Rubin, 2002). While the former are based on the assumption of unconfoundednesswith missing values, which can be difficult to assess in practice, the latter have been developed under strong parametric assumptions about the outcome, treatment and covariates models, in addition to relying on missing values hypotheses that can also be difficult to meet in practice (Yang et al., 2019).

To avoid relying on the hypothesis of unconfoundedness with missing values or being in the very parametric (and linear) framework of multiple imputation (Mattei & Mealli, 2009; Seaman & White, 2014) and matrix factorization (Kallus et al., 2018), we propose a new method for causal inference with missing data, which we call MissDeep- Causal. MissDeepCausal is inspired by the work of

Kallus et al. (2018) in the sense that we consider a model with latent confounders, and assume that we only have access to covariates with missing values that are noisy proxies of the true latent confounders. However, our approach generalizes and extends the work of Kallus et al. (2018) in different aspects: (i) instead of linear factor analysis models with missing values, we consider non-linear versions using deep latent variable models (Kingma & Welling, 2014; Rezende et al., 2014); (ii) we rely on the missing at random (MAR) (Rubin, 1976) assumption for the missing data mechanisms, and not on the stronger missing completely at random (MCAR) one; (iii) we take into account the posterior distribution of the latent variables given observed data and not only their conditional expectation. This latter point allows us to define a multiple imputation strategy adapted to the latent confounders model, and to couple it with doubly robust treatment effect estimation (Chernozhukov et al., 2018).

In the remainder of this article we first introduce the problem framework and recall existing work for handling missing values in causal inference in Section 2. We then introduce two variants of our MissDeepCausal approach in Section 3. Finally we compare MissDeepCausal empirically with several state-of-the-art methods on simulated data in Section 4.

In this section we start by quickly reviewing the problem of causal inference from observational data without missing data. We consider the potential outcomes framework (Rubin, 1974; Imbens & Rubin, 2015) where we have a sample of n independent and identically distributed (i.i.d.) observations with {0, 1} a binary treatment, a vector of covariates, and (the outcomes we would have observed had we assigned control or treatment to the i-th sample, respectively. The observed outcome for unit is defined as . The individual causal effect of the treatment is and the average treatment effect (ATE) is defined as

The ATE , i.e., the link between W and Y , can be estimated from observational data by taking into account the confounding factors X, i.e., the common causes of W and Y . A popular estimator of from observational data is the

so-called doubly robust estimator:

where are regression estimates of the conditional response surfaces {0, 1}, and is an estimate of the propensity score (Rosenbaum & Rubin, 1983; Imbens & Rubin, 2015).

Standard results state that if either or is correctly specified, then is an unbiased estimator of (Robins et al., 1994; Chernozhukov et al., 2018; Wager & Athey, 2018) under the following assumptions (Rosenbaum & Rubin, 1983): the ignorability or unconfoundedness assumption that states that all confounding factors are measured, i.e., conditionally on X, the treatment assignment is independent of the potential outcomes:

and the overlap assumption assuming the existence of some such that for all .

We now consider an extension to account for possible missing entries in the covariates. For that purpose, we denote the missingness pattern of the i-th sample as such that if is observed and otherwise. The matrix of observed covariates can be written as , with the elementwise multiplication and 1 the matrix filled with 1, so that takes its value in the half discrete space . We model as a random vector, and the possibility to infer causal effects with missing data now depends on additional assumptions on the joint law of . Methods for causal inference with missing covariates can be classified into two categories.

Rosenbaum & Rubin (1984) extend the unconfoundedness hypothesis (2) to missing values as

This implies the assumption, illustrated in Figure 1, that if a covariate is not observed, it is not a confounder. In particular, observations can have different confounders depending on their pattern of missing data. They define the generalized propensity score as:

which is a balancing score under (3). Consequently, an IPW estimator formed with estimators of can be an unbiased estimator of the ATE with missing values. Nevertheless, this method relies both on the fact that the covariates X are the appropriate set of confounders, which can be questioned without missing data (Kallus et al., 2018), and requires certain expert input and reasoning to verify that for each observation, treatment assignment and/or outcome values depend only on observed values of the confounders (Blake et al., 2019; Mayer et al., 2019b). Note in particular, that it is not because the missing data in the covariates are completely at random (MCAR), i.e., , that (3) is met. In practice, in addition, a difficulty with this approach is that estimating (4) requires fitting one model per pattern of missing values, which is unrealistic with classical tools (Miettinen, 1985; D’Agostino Jr & Rubin, 2000; D’Agostino Jr et al., 2001; Blake et al., 2019); Mayer et al. (2019b) address this problem using random forests adapted to covariates with missing values.

Figure 1. Unconfoundedness with missing values. X represents a the complete covariates, and M a missing data mechanism, represents the observed incomplete covariates, confounding the treatment assignment. The formalism of Pearl (1995) and

Missingness mechanisms assumptions. Multiple imputation is one of the most powerful approaches to estimate parameters and their variance from an incomplete data (Little & Rubin, 2002; van Buuren, 2018). Seaman & White (2014) show that when assuming (i) iden-tifiability of the ATE in the complete case, (ii) missing at random (MAR) values given W and Y , (iii) correct specifi-cation of the propensity score with logistic regression and of the Gaussian distribution of covariates, then multiple imputation gives a consistent estimate for the ATE estimated with IPW. An extension to doubly robust estimation has been proposed by Mayer et al. (2019b).

Instead of assuming that confounders are observed directly, Kallus et al. (2018) consider a more general model where observed covariates X are noisy and/or incomplete proxies of the true latent confounders Z. More specifically, they assume a low-rank model for the covariates and estimate the latent variables from the incomplete confounders using matrix completion methods (Hastie et al., 2015; Josse et al., 2016). Then, under the linear regression model

with random latent variables Z, missing values completely at random (MCAR) in X, unconfoundedness given Z, and some additional assumptions, they prove that regressing Y on and W leads to a consistent ATE estimator. Both techniques, multiple imputation and matrix factorization, rely on parametric (and linear) frameworks.

To avoid relying on the hypothesis of unconfoundedness with missing values (3) or being in the very parametric (and linear) framework of multiple imputation and matrix factorization, we propose MissDeepCausal, an approach based on deep latent variable models where the latent variables are assumed to be the confounders as represented in Figure 2.

Figure 2. Graphical representation of the model underlying MissDeepCausal. Z represents the unobserved latent confounders of the treatment W and the effect Y . X represents a proxy for the confounders, and M a missing data mechanism; the observed incomplete covariates.

Under this model, the unconfoundedness hypothesis (3) does not hold, so a standard treatment effect estimator using as covariates would be biased. On the other hand, we can express the treatment effect conditioned on as follows:

Consequently, if we have an unbiased estimator of , the treatment effect conditioned on Z, and if we know , the conditional distribution of Z given , then we can derive the treatment effect conditioned on by

Furthermore, by expressing the ATE as

we can form an estimate of the ATE by . We describe such an estimator in Section 3.2 below, which is reminiscent of multiple imputation techniques in the field of missing value imputation (Rubin, 1987).

Another strategy, described in Section 3.3, is to consider latent variables estimation as a pre-processing step prior to causal inference by computing

this can be seen as a non-linear extension of Kallus et al. (2018). Both estimators require sampling from the posterior distribution . Consequently, we first describe in Section 3.1 how to learn the joint distribution of (Z, X) from using a variational autoencoder (VAE) with missing data, before turning to the details of each strategy.

3.1. Deep latent variable models with missing values

Variational autoencoding Deep latent variable models can be defined as follows. Let be n i.i.d. random variables such that

The prior distribution of the latent variables or codes is often isotropic Gaussian . The function is a (deep) neural network called the decoder and is a parametric observation model, which we take to be multivariate Gaussian. The inference of deep latent variable models can be achieved by maximizing evidence lower bounds of the likelihood, such as the variational autoencoder bounds.

With missing values, the appropriate quantity to target for inference on , when the missing values mechanism can be ignored (Rubin, 1976; Little & Rubin, 2002), is the observed log-likelihood. Using Rubin (1976)’s notations, we define the partition of the data in realized observed and missing values given a specific realization of the pattern, it can be written as:

The corresponding evidence lower bound (ELBO) is:

with KL for the Kullback-Leibler divergence and the variational distribution

with the (parametric) variational distribution over . The function , called the encoder, is parametrized by a (deep) neural network whose weights are stored in .

To take into account missing values in deep latent variable models, Mattei & Frellsen (2019) suggest the missing data importance weight autoencoder bound (MIWAE) approach. They use a simple variational family where they impute the missing entries with a constant and show that using this class of distributions, it maximizes a lower bound of the observed log-likelihood. Specifically, they replace with

where is an imputation function chosen beforehand that transforms into a complete input vector .

Self-normalized importance sampling To estimate and sample from , we use the missing data importance weight autoencoder bound (MIWAE) approach of Mattei & Frellsen (2019), which is summarized above. They use a simple variational family where they impute the missing entries with a constant and show that using this class of distributions, it maximizes a lower bound of the observed log-likelihood. Note that their approach requires the classical missing at random (MAR) (Rubin, 1976) assumption to ignore the missing values mechanism when maximizing the observed likelihood for the VAE inference.

In the MIWAE approach, the variational distribution plays a central role but is not necessarily a good surrogate for the posterior distribution . To sample from the true posterior distribution, we resort to importance sampling techniques using the variational distribution for proposal. More precisely, we can define, for any measurable function s,

This quantity can be estimated using self-normalized importance sampling with:

Equation (8) is used in our second strategy described in Section 3.3, while for our first strategy (described in Section 3.2) we sample L samples according to , compute the weights as in (8) and re-sample B << L with probability proportional to the weights.

3.2. MissDeepCausal with multiple imputation (MDC-MI)

MDC-MI uses the importance sampling strategy presented in Section 3.1, to compute an approximation of (6) by Monte-Carlo as follows. First, we draw B i.i.d. samples from the posterior distribution . On each sample, we evaluate the function f and aggregate the results: . This approach can be viewed as a multiple imputation method, which consists in generating different imputed data sets by drawing the missing values from their posterior distribution given observed values, then estimating the parameters of interest on each imputed data set and aggregating the results according to Rubin’s rules (Rubin, 1987) to obtain a final estimate for the quantity of interest. Here we consider the samples of the latent variables and apply the doubly robust estimator from (1) on each table :

and get the final estimate for the causal effect by computing the mean of the estimators i.e. . The dou- bly robust estimator from (1) is asymptotically normal (under some mild assumptions) (Wager & Athey, 2018) which is required for the aggregation in multiple imputation procedures (Rubin, 1987). Note that this multiple imputation strategy additionally allows to reflect the variability due to the missing values in the variance estimation of the estimator .

3.3. MissDeepCausal with latent variables estimation as a pre-processing step (MDC-process)

We also propose MDC-process as a non-linear extension of Kallus et al. (2018), where we estimate defined in (7). For that purpose, we first approximate the expectation of the posterior distribution

to get estimates for the latent confounders. In a second step, we use them under the regression model (5) and accordingly regress the observed outcome Y on the estimated latent factors and the treatment assignment W to obtain an estimation of the treatment effect. This strategy is a heuristic extension of Kallus et al. (2018) to a non-linear case in the sense that the latent variables encode non-linear relationship between covariates.

An alternative, still heuristic, approach is to use the estimated latent confounders from (10) as inputs for standard techniques to estimate the average treatment effect. More precisely, for the doubly robust estimator (1), we replace the estimates for the propensity score with estimates for

and similarly for the conditional response surfaces.

However, note that this latter strategy would require from (10) to be a confounder instead of Z as it is assumed (see Figure 2).

4.1. Methods

We compare the following methods to handle missing values (the following acronyms are identical to the method labels used in Figures 3–7):

– MDC.process: using the estimations of the latent variables either in a regression adjustement estimator or in a double robust estimator as presented in Section 3.3;

– MDC.mi: using the doubly robust estimator MDC-mi Section 3.2.

Mattei & Frellsen (2019) to implement both methods. Throughout all experiments, we fix L = 10, 000 for the importance sampling weights. We choose hyperparameters, (variance of the prior on Z) and (dimension of estimated latent space), by cross-validation. We vary the number of draws B from the posterior for the MDC.mi approach from 50 to 500 (results only reported for B = 500).

Mattei & Mealli (2009) and Seaman & White (2014). We generate 20 imputations, using the python implementation in the scikit-learn (Pedregosa et al., 2011).

2018) based on nuclear norm penalty (python implementation inspired by the R package softImpute (Hastie & Mazumder, 2015)). The dimension of the latent space is chosen via cross-validation on the nuclear norm penalty parameter.

4.2. Settings

Under the latent confounding assumption (corresponding to the graphical model in Figure 2), we generate covariates according to two models:

• LRMF: The covariates are generated from a low-rank matrix factorization model as in Kallus et al. (2018).

• DLVM: The covariates are generated from a deep latent variable model as in as in Kingma & Welling (2014). , covariates are sampled from , where with drawn from standard Gaussian distributions and uniform distributions.

We define treatment and outcome models with a logistic-linear model as follows: and . We add an additive noise term in the outcome model such that the signal-to-noise ratio (SNR) is either 5 or 10 (results are only reported for the case SNR = 10).

Missing values are generated completely at random (MCAR), i.e., , with {0, 0.3, 0.5 0.9} and we consider the following problem dimensions: , and . Results are averaged over 30 replications for each setting. We only report results for n = 10, 000, experiments with other choices of parameters are reported in the Supplementary Material. Throughout all experiments the true ATE is fixed at 1.1

4.3. Results

4.3.1. REGRESSION ADJUSTMENT

First, we assess the quality of our heuristic described in Section 3.3 concerning the non-linear extension of Kallus et al. (2018). An estimation of is obtained by regressing the observed outcomes Y on the estimations of the latent factors Z (for MDC.process, MF) and on the imputed data (for MI).

Figures 3 and 4 show that our proposed method, MDC.process tends to slightly outperform all other methods when the covariates are generated according to a DLVM model. As expected the performances of all the methods decrease when the percentage of missing values increase, and both MF and MDC process better recover the latent structure when p is larger. Additionally we find that when the data is generated under the LRMF model, then our method performs as well as the initial proposal of Kallus et al. (2018) (results are reported in the Supplementary Material).

https://github.com/imkemayer/MissDeepCausal.

Figure 3. Estimated ATE via regression adjustment for varying amount of missing values; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z results are obtained using the true confounders Z. (n, p, d) = (10 000, 100, 2).

10 100 p

Figure 4. MSE of estimated ATE via regression adjustment for varying p; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z results are obtained using the true confounders Z. (n, d) = (10 000, 2).

Now we turn to the more flexible framework which does not assume linear relationships (5) between the outcome and the confounders. We consider the doubly robust estimator (1) with the (imputed) covariates X for MI and with the estimation of the latent variables Z for MF and MDC.

To estimate the regression surfaces and the propensity score e required for the doubly robust estimator (1), we use a logistic-linear model, either with or without additional regularization.

Figure 52 illustrates that even when the latent variables are generated from matrix factorization, our approaches based on the VAE with missing values lead to unbiased estimates. We note as well that all methods perform similarly, independently of the number of observed covariates p (results for the other values of p are in the Supplementary Material).

Figures 6 and 7 3 show that as expected, due to the flex-ibility of MissDeepCausal, the suggested approaches better handle highly non-linear relationships between the latent confounders and the observed (incomplete) covariates. MDC methods are the only ones achieving no biais or small bias under this non-linear model. This is all the more true as the number of variables p is large compared to the dimension of the latent space d. The matrix factorization approach fails in this setting to recover the confounders Z.

Figure 5. Estimated ATE via parametric doubly robust estimation for varying amount of missing values; covariates generated from a LRMF, (logistic-)linear model specification for (); results with Z are obtained using the true confounders Z. (n, p, d) = (10 000, 1 000, 2).

4.4. IHDP data

We assess our methodology on the Infant Health and Development Program (IHDP) benchmark data (Hill, 2011). The original data comes from a randomized control trial where the aim was to assess the impact of visits by specialists on children’s test scores. There are six quantitative and 19 binary variables, recorded for 985 individuals. Hill (2011) transformed the original experimental data into observational data by selecting a nonrandom subset among the treated, stratified along an ethnicity variable, which

Figure 6. Estimated ATE via parametric doubly robust estimation for varying amount of missing values; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z are obtained using the true confounders Z. (n, p, d) = (10 000, 1 000, 2).

Figure 7. MSE of estimated ATE via parametric doubly robust estimation for varying p; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z results are obtained using the true confounders Z. (n, d) = (10 000, 2).

leads to two unbalanced treatment groups. In total there are 139 treated and 608 control observations in the new data set. Then, keeping fixed the treatment variable, simulated data are obtained by generating new potential outcomes. More precisely, we follow the scenario “B” of (Hill, 2011) , i.e., and , with where is chosen to get an average treatment effect equal to 4. 4 After simulating the outcomes, we add missing values to the 25 covariates, assuming an MCAR mechanism. For the MIWAE part of our MDC methods, we select the parameters and by 5-fold cross-validation.

In addition to comparing the estimators considered in this paper that handle missing data, we also add two other approaches: the CEVAE estimator detailed in Louizos et al. (2017) as a baseline and the MIA.GRF estimator proposed in Mayer et al. (2019b). Note that CEVAE does not deal with missing values so that we replace the missing values by the mean of the variables. The CEVAE estimator is based on the difference between the two conditional expectations. The MIA.GRF estimator targets (4) and the generalized response surface analogue. It is based on estimation using random forests where missing values are encoded with missing incorporated in attributes such that the splitting rules in the random forests exploit the missingness pattern (Twala et al., 2008; Josse et al., 2019). We use the R package grf (Tibshirani et al., 2018) for the complete case and the implementation provided by Mayer et al. (2019b) for the incomplete case5.

Finally, we additionally apply a nonparametric doubly robust estimator, denoted by DR, on the approximated confounders (resp. imputed covariates) based on (generalized) random forests (Athey et al., 2019). For this part we use the implementation of the R package grf (Tibshirani et al., 2018).

For comparability with previous experiments on these data, we report the in-sample mean absolute error, i.e. the mean absolute difference between the estimated ATE and the sample ATE (by construction of the data we know the exact values of and for all .

Table 1 shows that the doubly robust estimators (either in the parametric regression, , or the random forest form ) systematically outperform the corresponding OLS estimator which highlights that the linear model is not appropriate, at least that it is not linear in the covariates X. Indeed, we know that the outcome is simulated as a non-linear function of the (complete) covariates X, whereas the treatment assignment is taken from the (de-randomized) experiment and can therefore well depend on latent variables. The results of MissDeepCausal are competitive with other approaches and greatly improve on CEVAE and MI. Its performances when used with the double robust estimators are stable with respect to the percentage of missing values.

Table 1. Methods on the IHDP benchmark data. Mean absolute error (with standard error) across simulations on all the data points (in-sample error). OLS corresponds to the estimator obtained by regression and DR to the doubly robust estimator(s).

In this work we have investigated the problem of treatment effect estimation with incomplete covariates. This problem of missing values is highly relevant for modern causal inference as it is exacerbated with high dimensional data. Yet most causal inference techniques do not address this issue; and complete case analysis, in addition to leading to potentially inconsistent causal effects estimators, is not an option anymore. We have proposed MissDeepCausal which borrows the strength of deep latent variable models to retrieve the latent confounders from incomplete covariates encoding complex non-linear relationships. We use a modular approach in the style of Bayesian propensity based methods for treatment effect estimation (Zigler, 2016), where the latent variables are used as inputs for doubly robust estimators. We suggest a multiple imputation strategy that allows to fully exploit the posterior distribution of the latent variables. Numerical results are very encouraging insofar as we obtain best relative performance in terms of bias and MSE whether the underlying model is well or badly specified compared to current state of the art. Open challenges include heterogeneous treatment effect estimation with missing values as well as the ambitious task of handling missing not at random type data.

Athey, S., Tibshirani, J., Wager, S., et al. Generalized random forests. The Annals of Statistics, 47(2):1148–1178, 2019.

Blake, H. A., Leyrat, C., Mansfield, K., Seaman, S., Tomlinson, L., Carpenter, J., and Williamson, E. Propensity scores using missingness pattern information: a practical guide. arXiv preprint, 2019.

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1–C68, 2018.

D’Agostino Jr, R., Lang, W., Walkup, M., Morgan, T., and Karter, A. Examining the impact of missing data on propensity score estimation in determining the effectiveness of self-monitoring of blood glucose (smbg). Health Services and Outcomes Research Methodology, 2(3-4): 291–315, 2001. doi: 10.1023/A:102037541.

D’Agostino Jr, R. B. and Rubin, D. B. Estimating and using propensity scores with partially missing data. Journal of the American Statistical Association, 95(451):749–759, 2000. doi: 10.2307/2669455.

Guo, R., Cheng, L., Li, J., Hahn, P. R., and Liu, H. A survey of learning causality with data: Problems and methods. arXiv preprint arXiv:1809.09337, 2019.

Hastie, T. and Mazumder, R. softImpute: Matrix Completion via Iterative Soft-Thresholded SVD, 2015. R package version 1.4.

Hastie, T., Mazumder, R., Lee, J. D., and Zadeh, R. Matrix completion and low-rank svd via fast alternating least squares. Journal of Machine Learning Research, 16: 3367–3402, 2015.

Hill, J. L. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1):217–240, 2011.

Horvitz, D. G. and Thompson, D. J. A generalization of sampling without replacement from a finite universe. Journal of the American statistical Association, 47(260):663–685, 1952. doi: 10.1080/01621459.1952. 10483446.

Iacus, S. M., King, G., and Porro, G. Causal inference without balance checking: Coarsened exact matching. Political Analysis, 20(1):1–24, 2012. doi: 10.1093/pan/ mpr013.

Imbens, G. W. Nonparametric estimation of average treatment effects under exogeneity: A review. Review of Economics and statistics, 86(1):4–29, 2004.

Imbens, G. W. and Rubin, D. B. Causal inference in statistics, social, and biomedical sciences. Cambridge Univer- sity Press, 2015.

Johansson, F., Shalit, U., and Sontag, D. Learning representations for counterfactual inference. In International conference on machine learning, pp. 3020–3029, 2016.

Josse, J. and Reiter, J. P. Introduction to the special section on missing data. Statist. Sci., 33(2):139–141, 05 2018. doi: 10.1214/18-STS332IN.

Josse, J., Sardy, S., and Wager, S. denoiser: A package for low rank matrix estimation. Journal of Statistical Software, 2016.

Josse, J., Prost, N., Scornet, E., and Varoquaux, G. On the consistency of supervised learning with missing values. arXiv preprint, 2019.

Kallus, N., Mao, X., and Udell, M. Causal inference with noisy and missing covariates via matrix factorization. In Advances in Neural Information Processing Systems, pp. 6921–6932, 2018.

Kingma, D. P. and Welling, M. Auto-encoding variational bayes. In International Conference on Learning Representations, 2014.

Kuroki, M. and Pearl, J. Measurement bias and effect restoration in causal inference. Biometrika, 101(2):423– 437, 2014.

Little, R. J. A. and Rubin, D. B. Statistical Analysis with Missing Data. Wiley, 2002. ISBN 0471183865. doi: 10.2307/1533221.

Louizos, C., Shalit, U., Mooij, J. M., Sontag, D., Zemel, R., and Welling, M. Causal effect inference with deep latent-variable models. In Advances in Neural Information Processing Systems, pp. 6446–6456, 2017.

Lunceford, J. K. and Davidian, M. Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study. Statistics in medicine, 23(19):2937–2960, 2004. doi: 10.1002/sim. 1903.

Mattei, A. and Mealli, F. Estimating and using propensity score in presence of missing background data: an application to assess the impact of childbearing on wellbeing. Statistical Methods and Applications, 18(2):257– 273, 2009. doi: 10.1007/s10260-007-0086-0.

Mattei, P.-A. and Frellsen, J. MIWAE: Deep generative modelling and imputation of incomplete data sets. In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 4413–4423, Long Beach, California, USA, 09–15 Jun 2019. PMLR.

Mayer, I., Josse, J., Tierney, N., and Vialaneix, N. R-miss-tastic: a unified platform for missing values methods and workflows. arXiv preprint arXiv:1908.04822, 2019a.

Mayer, I., Wager, S., Gauss, T., Moyer, J.-D., and Josse, J. Doubly robust treatment effect estimation with missing attributes. arXiv preprint arXiv:1910.10624, 2019b.

Miettinen, O. S. Theoretical epidemiology: principles of occurrence research. John Wiley & Sons, New York, 1985. ISBN 0827343132.

Pearl, J. Causal diagrams for empirical research. Biometrika, 82(4):669–688, 1995.

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011.

Rezende, D. J., Mohamed, S., and Wierstra, D. Stochastic backpropagation and approximate inference in deep generative models. pp. 1278–1286, 2014.

Richardson, T. S. and Robins, J. M. Single world intervention graphs (swigs): A unification of the counterfactual and graphical approaches to causality. Technical report, Center for Statistics and the Social Sciences, University of Washington, 2013.

Robins, J. M., Rotnitzky, A., and Zhao, L. P. Estimation of regression coefficients when some regressors are not always observed. Journal of the American Statistical Association, 89(427):846–866, 1994. doi: 10.1080/ 01621459.1994.10476818.

Rosenbaum, P. R. and Rubin, D. B. The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1):41–55, 1983. doi: 10.1093/ biomet/70.1.41.

Rosenbaum, P. R. and Rubin, D. B. Reducing bias in observational studies using subclassification on the propensity score. Journal of the American Statistical Association, 79(387):516–524, 1984. doi: 10.2307/2288398.

Rubin, D. B. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5):688–701, 1974. doi: 10.1037/ h0037350.

Rubin, D. B. Inference and missing data. Biometrika, 63 (3):581–592, 1976.

Rubin, D. B. Multlipe Imputation for Nonresponse in Surveys. Wiley, Hoboken, NJ, USA, 1987. ISBN 9780471655740.

Seaman, S. and White, I. Inverse probability weighting with missing predictors of treatment assignment or missingness. Communications in Statistics-Theory and Methods, 43(16):3499–3515, 2014. doi: 10.1080/03610926. 2012.700371.

Shalit, U., Johansson, F. D., and Sontag, D. Estimating individual treatment effect: generalization bounds and algorithms. In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pp. 3076–3085. JMLR. org, 2017.

Tibshirani, J., Athey, S., Wager, S., Friedberg, R., Miner, L., and Wright, M. grf: Generalized Random Forests (Beta), 2018. R package version 0.10.2.

Twala, B., Jones, M., and Hand, D. J. Good methods for coping with missing data in decision trees. Pattern Recognition Letters, 29(7):950–956, 2008.

van Buuren, S. Flexible Imputation of Missing Data. Chapman and Hall/CRC, Boca Raton, FL, 2018.

Wager, S. and Athey, S. Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523):1228– 1242, 2018.

Yang, S., Wang, L., and Ding, P. Causal inference with confounders missing not at random. Biometrika, 2019.

Yoon, J., Jordon, J., and van der Schaar, M. Ganite: Estimation of individualized treatment effects using generative adversarial nets. 2018.

Zigler, C. M. The central role of bayes theorem for joint estimation of causal effects and propensity scores. The American Statistician, 70(1):47–54, 2016.

In this supplementary material we provide additional results for our simulation study. Namely, Figures S.8 and S.9 are complementary to Figures 4 and 7 respectively where instead of varying the dimension of the ambient space p, we vary the number of observations n. Figure S.10 illustrates the result of similar performance of all methods in the case of data generated under a LRMF model mentioned in Section 4.

Figure S.8. MSE of estimated ATE via regression adjustment for varying n; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z results are obtained using the true confounders Z. (p, d) = (100, 2).

Figure S.9. MSE of estimated ATE via parametric doubly robust estimation for varying n; covariates generated from a DLVM, (logistic-)linear model specification for (); results with Z are obtained using the true confounders Z. (p, d) = (1 000, 2).

Figure S.10. Estimated ATE via regression adjustment for varying amount of missing values; covariates generated from a LRMF, (logistic-)linear model specification for (); results with Z results are obtained using the true confounders Z. (n, p, d) = (10 000, 100, 2).