summaries: This paper introduces MissDeepCausal method to address the problem of treatment effect estimation with incomplete covariates matrix (missing values at random -- MAR). It makes use of Variational AutoEncoders (VAE) to learn the latent confounders from incomplete covariates. This also helps encoding complex non-linear relationships in the data, a capability that is missing in the work of Kallus et al. (2018) -- the work which this paper extends. They employ the Missing data Importance Weight AutoEncoder (MIWAE) approach (Mattei & Frellsen, 2019) to approximate the posterior of their latent factors Z given the observed incomplete covariates X*. The main contributions of this work are presented in sections 3.2 and 3.3, where they use the approximated posterior derived from MIWAE to sample Z to be used for estimating outcomes and finally calculating the Average Treatment Effect (ATE). This is done according to the doubly robust estimator developed for data with incomplete covariates (Mayer et al., 2019b).         The paper considers average treatment effect estimation treatment T and an unobserved confounder Z causes the outcome Y with an added constraint that the observed X is a noisy measurement of the underlying Z and some of the entries of observed X are missing at random (MAR). Previous work (Kallus et al. 2018) on similar settings assumed a low rank model connecting Z and X along with some entries missing at random which we do not observe. Further, Y (outcome) is related to the treatment and Z with a linear model. They actually show that matrix factorization techniques with these assumptions form an unbiased estimator for the Average Treatment Effect. There is prior work on doubly robust estimators under ignorability assumptions. This contribution considers deep latent-factor models for causal inference -inferring the effect of a treatment- in the presence of missing values. The core challenge is that of confounders not directly observed, and accessible only via noisy proxys, in particular in the missingness. The contributed method relies on using the latent-factor model for multiple imputation in doubly-robust causal treatment effect estimators. As a consequence, it requires the missing at random assumption to control for the impact of imputation. Given that the confounders are not directly assumed, a first approach estimates their effect via an estimate of P(Z|X*)  (probability of confounder given observed data), which is then plugged in the doubly robust estimator in a multiple imputation strategy. A second approach uses heuristically the estimated latent confounders as regressors of non interest in a linear-regression model. The approaches are empirically compared to other imputation strategies used as plugins in the doubly-robust estimator. The contributed approach show marked benefits when the problem is highly non-linear.
abstractText: Inferring causal effects of a treatment, intervention or policy from observational data is central to many applications. However, state-of-the-art methods for causal inference seldom consider the possibility that covariates have missing values, which is ubiquitous in many real-world analyses. Missing data greatly complicate causal inference procedures as they require an adapted unconfoundedness hypothesis which can be difficult to justify in practice. We circumvent this issue by considering latent confounders whose distribution is learned through variational autoencoders adapted to missing values. They can be used either as a pre-processing step prior to causal inference but we also suggest to embed them in a multiple imputation strategy to take into account the variability due to missing values. Numerical experiments demonstrate the effectiveness of the proposed methodology especially for non-linear models compared to competitors.
references.authors: Susan Athey Julie Tibshirani Stefan Wager
references.title: Generalized random forests
references.venue: The Annals of Statistics,
references.authors: Helen A. Blake Clmence Leyrat Kate Mansfield Shaun Seaman Laurie Tomlinson James Carpenter Elizabeth Williamson
references.title: Propensity scores using missingness pattern information: a practical guide
references.authors: Victor Chernozhukov Denis Chetverikov Mert Demirer Esther Duflo Christian Hansen Whitney Newey James Robins
references.title: Double/debiased machine learning for treatment and structural parameters
references.venue: The Econometrics Journal,
references.authors: Ralph D’Agostino Jr. Wei Lang Michael Walkup Timothy Morgan Andrew Karter
references.title: Examining the impact of missing data on propensity score estimation in determining the effectiveness of selfmonitoring of blood glucose (smbg)
references.venue: Health Services and Outcomes Research Methodology,
references.authors: Ralph B D’Agostino Jr. Donald B Rubin
references.title: Estimating and using propensity scores with partially missing data
references.venue: Journal of the American Statistical Association,
references.authors: Ruocheng Guo Lu Cheng Jundong Li P. Richard Hahn Huan Liu
references.title: A survey of learning causality with data: Problems and methods
references.venue: arXiv preprint arXiv:1809.09337,
references.authors: Trevor Hastie Rahul Mazumder
references.title: softImpute: Matrix Completion via Iterative Soft-Thresholded SVD
references.authors: Trevor Hastie Rahul Mazumder Jason D. Lee Reza Zadeh
references.title: Matrix completion and low-rank svd via fast alternating least squares
references.venue: Journal of Machine Learning Research,
references.authors: Jennifer L Hill
references.title: Bayesian nonparametric modeling for causal inference
references.venue: Journal of Computational and Graphical Statistics,
references.authors: Daniel G Horvitz Donovan J Thompson
references.title: A generalization of sampling without replacement from a finite universe
references.venue: Journal of the American statistical Association,
references.authors: Stefano M Iacus Gary King Giuseppe Porro
references.title: Causal inference without balance checking: Coarsened exact matching
references.venue: Political Analysis,
references.authors: Guido W Imbens
references.title: Nonparametric estimation of average treatment effects under exogeneity: A review
references.venue: Review of Economics and statistics,
references.authors: Guido W Imbens Donald B Rubin
references.title: Causal inference in statistics, social, and biomedical sciences
references.authors: Fredrik Johansson Uri Shalit David Sontag
references.title: Learning representations for counterfactual inference
references.venue: In International conference on machine learning,
references.authors: Julie Josse Jerome P. Reiter
references.title: Introduction to the special section on missing data
references.venue: Statist. Sci.,
references.authors: Julie Josse Sylvain Sardy Stefan Wager
references.title: denoiser: A package for low rank matrix estimation
references.venue: Journal of Statistical Software,
references.authors: Julie Josse Nicolas Prost Erwan Scornet Gaël Varoquaux
references.title: On the consistency of supervised learning with missing values
references.venue: arXiv preprint,
references.authors: N. Kallus X. Mao M. Udell
references.title: Causal inference with noisy and missing covariates via matrix factorization
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Diederik P Kingma Max Welling
references.title: Auto-encoding variational bayes
references.venue: In International Conference on Learning Representations,
references.authors: Manabu Kuroki Judea Pearl
references.title: Measurement bias and effect restoration
references.venue: in causal inference. Biometrika,
references.authors: R.J.A. Little D.B. Rubin
references.title: Statistical Analysis with Missing Data
references.venue: ISBN 0471183865
references.authors: Christos Louizos Uri Shalit Joris M Mooij David Sontag Richard Zemel Max Welling
references.title: Causal effect inference with deep latent-variable models
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Jared K Lunceford Marie Davidian
references.title: Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study
references.venue: Statistics in medicine,
references.authors: Alessandra Mattei Fabrizia Mealli
references.title: Estimating and using propensity score in presence of missing background data: an application to assess the impact of childbearing on wellbeing
references.venue: Statistical Methods and Applications,
references.authors: Pierre-Alexandre Mattei Jes Frellsen
references.title: MIWAE: Deep generative modelling and imputation of incomplete data sets
references.venue: Proceedings of the 36th International Conference on Machine Learning,
references.authors: Imke Mayer Julie Josse Nicholas Tierney Nathalie
references.title: Vialaneix. R-miss-tastic: a unified platform for missing values methods and workflows
references.venue: arXiv preprint arXiv:1908.04822,
references.authors: Imke Mayer Stefan Wager Tobias Gauss Jean-Denis Moyer Julie Josse
references.title: Doubly robust treatment effect estimation with missing attributes
references.venue: arXiv preprint arXiv:1910.10624,
references.authors: Olli S. Miettinen
references.title: Theoretical epidemiology: principles of occurrence research
references.authors: Judea Pearl
references.title: Causal diagrams for empirical research
references.venue: Biometrika, 82(4):669–688,
references.authors: R R Core Team
references.title: A Language and Environment for Statistical Computing
references.venue: R Foundation for Statistical Computing, Vienna,
references.authors: Danilo Jimenez Rezende Shakir Mohamed Daan Wierstra
references.title: Stochastic backpropagation and approximate inference in deep generative models
references.authors: Thomas S Richardson James M Robins
references.title: Single world intervention graphs (swigs): A unification of the counterfactual and graphical approaches to causality
references.venue: Technical report, Center for Statistics and the Social Sciences, University of Washington,
references.authors: Geneviève Robin
references.title: mimi: Main Effects and Interactions in Mixed and Incomplete Data, 2019. R package version 0.2.0
references.authors: J.M. Robins A. Rotnitzky L.P. Zhao
references.title: Estimation of regression coefficients when some regressors are not always observed
references.venue: Journal of the American Statistical Association,
references.authors: Paul R Rosenbaum Donald B Rubin
references.title: The central role of the propensity score in observational studies for causal effects
references.authors: Paul R Rosenbaum Donald B Rubin
references.title: Reducing bias in observational studies using subclassification on the propensity score
references.venue: Journal of the American Statistical Association,
references.authors: D.B. Rubin
references.title: Estimating causal effects of treatments in randomized and nonrandomized studies
references.venue: Journal of Educational Psychology,
references.authors: D.B. Rubin
references.title: Multlipe Imputation for Nonresponse in Surveys
references.authors: Shaun Seaman Ian White
references.title: Inverse probability weighting with missing predictors of treatment assignment or missingness
references.venue: Communications in Statistics-Theory and Methods,
references.authors: Uri Shalit Fredrik D Johansson David Sontag
references.title: Estimating individual treatment effect: generalization bounds and algorithms
references.venue: In Proceedings of the 34th International Conference on Machine Learning-Volume
references.authors: Julie Tibshirani Susan Athey Stefan Wager Rina Friedberg Luke Miner Marvin Wright
references.title: grf: Generalized Random Forests (Beta), 2018
references.venue: R package version
references.authors: BETH Twala MC Jones David J Hand
references.title: Good methods for coping with missing data in decision trees
references.venue: Pattern Recognition Letters,
references.authors: S. van Buuren
references.title: Flexible Imputation of Missing Data
references.authors: Stef van Buuren Karin Groothuis-Oudshoorn
references.title: mice: Multivariate imputation by chained equations in r
references.venue: Journal of Statistical Software,
references.authors: Stefan Wager Susan Athey
references.title: Estimation and inference of heterogeneous treatment effects using random forests
references.venue: Journal of the American Statistical Association,
references.authors: Shu Yang Linbo Wang Peng Ding
references.title: Causal inference with confounders missing not at random
references.authors: Jinsung Yoon James Jordon Mihaela van der Schaar
references.title: Ganite: Estimation of individualized treatment effects using generative adversarial nets
references.authors: Corwin Matthew Zigler
references.title: The central role of bayes theorem for joint estimation of causal effects and propensity scores
references.venue: The American Statistician,
sections.heading: 1 INTRODUCTION
sections.text: 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 realworld 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 unconfoundedness with missing values, which can be difficult to assess in practice, the latter have been developed under strong para-
  metric 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 MissDeepCausal. 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 nonlinear 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.
sections.heading: 2 SETTING, NOTATIONS AND RELATED WORKS
sections.text: 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 (Yi(0), Yi(1),Wi, Xi)i=1, ..., n with Wi ∈ {0, 1} a binary treatment, Xi = (Xi1, . . . , Xip) > ∈ Rp a vector of covariates, and (Yi(0), Yi(1)) ∈ R2 the outcomes we would have observed had we assigned control or treatment to the i-th sample, respectively. The observed outcome for unit i, Yi ∈ R is defined as Yi , WiYi(1) + (1 − Wi)Yi(0). The individual causal effect of the treatment is τi , Yi(1)− Yi(0) and the average treatment effect (ATE) is defined as
  τ , E[Yi(1)− Yi(0)] = E[τi].
  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:
  τ̂DR , 1
  n n∑ i=1 µ̂1(Xi)− µ̂0(Xi) +Wi Yi − µ̂1(Xi) ê(Xi) − (1−Wi) Yi − µ̂0(Xi) 1− ê(Xi) , (1)
  where µ̂w(x) are regression estimates of the conditional response surfaces µw(x) , E[Y (w) |X = x], w ∈ {0, 1}, and ê(x) is an estimate of the propensity score e(x) , P(Wi = 1 |Xi = x) (Rosenbaum & Rubin, 1983; Imbens & Rubin, 2015).
  Standard results state that if either (µ̂0, µ̂1) or ê is correctly specified, then τ̂DR 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:
  {Yi(1), Yi(0)} ⊥ Wi |Xi, for all i; (2)
  and the overlap assumption assuming the existence of some η > 0 such that η < e(x) < 1 − η, for all x ∈ X . 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 Mi ∈ {0, 1}p such that Mij = 0
  if Xij is observed and Mij = 1 otherwise. The matrix of observed covariates can be written as X? , X (1 −M) + NA M , with the elementwise multiplication and 1 the matrix filled with 1, so that X? takes its value in the half discrete space X ? , (R ∪ {NA})p. We model Mi as a random vector, and the possibility to infer causal effects with missing data now depends on additional assumptions on the joint law of (Yi(0), Yi(1),Wi, Xi,Mi)i=1, ..., n. Methods for causal inference with missing covariates can be classified into two categories.
  Unconfoundedness with missing values. Rosenbaum & Rubin (1984) extend the unconfoundedness hypothesis (2) to missing values as
  {Yi(1), Yi(0)} ⊥ Wi |X?i , for all i. (3)
  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:
  ∀x? ∈ X ?, e?(x?) , P(Wi = 1 |X?i = x?) , (4)
  which is a balancing score under (3). Consequently, an IPW estimator formed with estimators of e? 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., M ⊥ X , 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.
  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) identifiability of the ATE in the complete case, (ii) missing at random (MAR) values givenW and Y , (iii) correct specification 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
  Yi = Z T i α+ τWi + εi, (5)
  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 Ẑ and W leads to a consistent ATE estimator. Both techniques, multiple imputation and matrix factorization, rely on parametric (and linear) frameworks.
sections.heading: 3 MISSDEEPCAUSAL
sections.text: 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.
  Under this model, the unconfoundedness hypothesis (3) does not hold, so a standard treatment effect estimator using X∗ as covariates would be biased. On the other hand, we can express the treatment effect conditioned on X∗ as follows:
  E[Y (1)− Y (0) |X∗] = E[E[Y (1)− Y (0) |Z,X∗] |X∗] = E[E[Y (1)− Y (0) |Z] |X∗] .
  Consequently, if we have an unbiased estimator f̂(Z) of E[Y (1) − Y (0) |Z], the treatment effect conditioned on Z, and if we know P (Z |X∗), the conditional distribution of Z given X∗, then we can derive the treatment effect conditioned on X∗ by
  ĝ(X?) , E[f̂(Z)|X?] . (6)
  Furthermore, by expressing the ATE as
  τ = E[Y (1)− Y (0)] = E[E[Y (1)− Y (0) |X∗]] ,
  we can form an estimate of the ATE by E[ĝ(X?)]. 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 preprocessing step prior to causal inference by computing
  h(X?) , f(E[Z|X?]); (7)
  this can be seen as a non-linear extension of Kallus et al. (2018). Both estimators require sampling from the posterior distribution P (Z|X?). Consequently, we first describe in Section 3.1 how to learn the joint distribution of (Z,X) from X? using a variational autoencoder (VAE) with missing data, before turning to the details of each strategy.
sections.heading: 3.1 DEEP LATENT VARIABLE MODELS WITH MISSING VALUES
sections.text: To estimate and sample from P (Z |X?), we use the missing data importance weight autoencoder bound (MIWAE) approach of Mattei & Frellsen (2019), which is summarized in the appendix A. 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 Qγ(Z|X?) (defined in Appendix A) plays a central role but is not necessarily a good surrogate for the posterior distribution Pθ(Z|X?). To
  sample from the true posterior distribution, we resort to importance sampling techniques using the variational distribution Qγ for proposal. More precisely, we can define, for any measurable function s,
  E[s(Z)|X?] = ∫ s(Z)pθ(Z|X?)dZ = 1
  p(X?)
  ∫ s(Z) pθ(X ?|Z)p(Z)
  qγ(Z|X?) qγ(Z|X?)dZ.
  This quantity can be estimated using self-normalized importance sampling with:
  E[s(Z)|X?] ≈ L∑ l=1 wls(Z (l)), where wl , rl r1 + ...+ rL , with rl , pθ(X ?|Z(l))p(Z(l)) qγ(Z(l)|X?) . (8)
  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 Z(1), . . . , Z(L) according to Qγ(Z|X?), compute the weights as in (8) and re-sample B << L with probability proportional to the weights.
sections.heading: 3.2 MISSDEEPCAUSAL WITH MULTIPLE IMPUTATION (MDC-MI)
sections.text: 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 (Z(j))1≤j≤B ∈ Rn×d from the posterior distribution P (Z|X?). On each sample, we evaluate the function f and aggregate the results: ĝ(B)(X?) = 1B ∑B j=1 f(Z
  (j)). 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 Z(j) of the latent variables and apply the doubly robust estimator from (1) on each table Z(j):
  τ̂ (j) = 1
  n n∑ i=1 µ̂ (j) 1 (Z (j) i )− µ̂ (j) 0 (Z (j) i ) +Wi Yi − µ̂(j)1 (Z (j) i ) ê(j)(Z (j) i ) − (1−Wi) Yi − µ̂(j)0 (Z (j) i ) 1− ê(j)(Z(j)i ) , (9)
  and get the final estimate for the causal effect by computing the mean of the estimators i.e. τ̂ = 1 B ∑B j=1 τ̂
  (j). The doubly 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 τ̂ .
sections.heading: 3.3 MISSDEEPCAUSAL WITH LATENT VARIABLES ESTIMATION AS A PRE-PROCESSING STEP (MDC-PROCESS)
sections.text: We also propose MDC-process as a non-linear extension of Kallus et al. (2018), where we estimate h(X?) defined in (7). For that purpose, we first approximate the expectation of the posterior distribution Ẑ(x?) , E[Z|X? = x?] (10) 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 Ẑ(x?) 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
  ẽ(z) = P(Wi = 1 | Ẑi(x?) = z),
  and similarly for the conditional response surfaces.
  However, note that this latter strategy would require Ẑ(x?) from (10) to be a confounder instead of Z as it is assumed (see Figure 2).
sections.heading: 4 SIMULATION STUDY
sections.heading: 4.1 METHODS
sections.text: We compare the following methods to handle missing values (the following acronyms are identical to the method labels used in Figures 3–5):
  • MissDeepCausal: – MDC.process: using the doubly robust estimator MDC-process from Section 3.3; – MDC.mi: using the doubly robust estimator MDC-mi Section 3.2.
  We extended the publicly available code of Mattei & Frellsen (2019) to implement both methods, using the default parameters for the MIWAE part. Using notations of Appendix A, we use L = 10, 000 for the importance sampling weights. In addition, for MDC-mi, we sample B = 200 observations from the estimated posterior distribution of Z|X?. • MICE: the multiple imputation approach as suggested in Mattei & Mealli (2009) and Sea-
  man & White (2014). We generate 10 imputations, using the implementation in the R (R Core Team, 2018) package mice (van Buuren & Groothuis-Oudshoorn, 2011). • MF: the matrix factorization approach of Kallus et al. (2018) using R packages softImpute (Hastie & Mazumder, 2015) (for continuous data) and mimi (Robin, 2019) (for mixed data) for the matrix completion based on nuclear norm penalty. We choose the dimension of the latent space via cross-validation.
  • MIA: the method proposed by Mayer et al. (2019b) which 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 case1.
sections.heading: 4.2 SETTINGS
sections.text: 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). Zi ∼ Nd(0, 1), covariates Xi are sampled from Np(µ(Z),Σ(Z)), where (µ(Z),Σ(Z)) = (V tanh(UZ + a) + b,diag{exp(ηT tanh(UZ + a) + δ)}) with U, V, a, b, δ, η drawn from Standard Gaussian distributions and Uniform distributions.
  Missing values are generated completely at random (MCAR), i.e., P(Mij = 1) = ρ, ∀ i, ∀ j, with ρ ∈ {0, 0.1, 0.3} and we consider the following problem dimensions: n = 10, 000, p = 10, and d = 3. Results are reported using 20 simulations for each setting. Throughout all experiments the true ATE τ is fixed at 1. Additional experiments with other choices of parameters are reported in the Appendix B.
sections.heading: 4.3 RESULTS
sections.heading: 4.3.1 REGRESSION ADJUSTMENT
sections.text: First, we assess the quality of our heuristic described in Section 3.3 concerning the non-linear extension of Kallus et al. (2018). For this we define treatment and outcome models with a logistic-linear model as follows: logit(e(Zi·)) = αTZi· and Yi ∼ N ((βTZi + τWi, σ2). An estimation of τ is obtained by regressing the observed outcomes Yi on the estimations of the latent factors Z (for MDC.process, MF) and on the imputed data Ximp (for MICE).
  1https://github.com/imkemayer/causal-inference-missing
  Figure 3 shows that as expected our proposed method, lin.MDC.process actually outperforms all other methods when the covariates are generated according to a DLVM model. The bias observed is small with respect to the one exhibited for completely observed covariates X . Additionally we observe that if the data is generated under the LRMF model, then our method performs as well as the initial proposal of Kallus et al. (2018) (results for this are not reported here).
sections.heading: 4.3.2 DOUBLY ROBUST ESTIMATION
sections.text: Now we turn to the more flexible framework which does not assume linear relationships (equation 5) between the outcome and the confounders. We generate the treatment and outcome using logisticlinear models and non-linear models based on non-linear transformations of the latent variables Z. We consider the doubly robust estimator (1) with the (imputed) covariates X for MICE and with the estimation of the latent variables Z for MF and MissDeepCausal.
  To estimate the regression surfaces (µ1, µ0) and the propensity score e required for the doubly robust estimator (equation 1), we use either a logistic-linear model or (generalized) random forests (Athey et al., 2019), indicated respectively by the prefixes loglin and grf in all figures.2 For the latter, we use the implementation of the R package grf (Tibshirani et al., 2018).
  Figure 4 illustrates that even when the latent variables are generated from matrix factorization, our approaches based on the VAE with missing values lead to estimates that are almost unbiased, given that the estimation of the propensity score and response surfaces is adapted to the considered models. The small bias observed for the matrix factorization pre-processing approach from Kallus et al. (2018) is not in contradiction with their theory since we use the doubly robust estimators and not the regression model. In addition, they require a (much) larger number of proxy covariates w.r.t. the number of latent confounders.
  Figure 5 shows that as expected, due to the flexibility of MissDeepCausal, the suggested approaches better handle highly non-linear relationships between the latent confounders and the observed (incomplete) covariates. It turns out that the multiple imputation strategy is particularly appropriate when the relationships between the outcome, the treatment and the confounders are highly nonlinear.
sections.heading: 4.4 IHDP DATA
sections.text: 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 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
  2More specifically, for estimating µw, we use all observations having Wi = w. 2In Figure 4b, all loglin estimations yield values around 6 and are therefore omitted for better readability.
  data are obtained by generating new potential outcomes. More precisely, we follow the scenario “B” of Hill (2011) , i.e., Y (0) ∼ N (µ0, 1) and Y (1) ∼ N (µ1, 1), with (µ0, µ1) = (exp(X + W )β, Xβ−ω) where ω is chosen to get an average treatment effect τ equal to 4. 3 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 σprior and d by 5-fold cross-validation.
  In addition to comparing the estimators considered in this paper that handle missing data, we also add the CEVAE estimator detailed in Louizos et al. (2017) as a baseline. 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. For comparability with their method and 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 µ(1)(Xi) and µ(0)(Xi) for all i): ∆ =
  ∣∣τ̂ − 1n∑µ(1)(Xi)− µ(0)(Xi)∣∣. Table 1 shows that the doubly robust estimators (either in the parametric regression, DRlog−lin, or the random forest form DRrf ) 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 nonlinear function of the (complete) covariates X , whereas the treatment assignment is taken from the (derandomized) experiment and can therefore well depend on latent variables. The results of MissDeepCausal are competitive with other
  3We use and adapt the corresponding code from V. Dorie: https://github.com/vdorie/npci/.
  approaches and greatly improve on CEVAE and MICE. Its performances when used with the double robust estimators are stable with respect to the percentage of missing values.
sections.heading: 5 CONCLUSION
sections.text: 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 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.
sections.heading: A APPENDIX
sections.text: A.1 DEEP LATENT VARIABLE MODELS WITH MISSING VALUES
  Deep latent variable models can be defined as follows. Let (Xi, Zi)i≤n be n i.i.d. random variables such that {
  Zi ∼ P (Zi) Xi ∼ Pθ(Xi|Zi) = Φ (Xi|fθ(Zi)).
  The prior distribution of the latent variables or codes Zi ∈ Rd is often isotropic Gaussian Zi ∼ N (0d, Id). The function fθ : Rd → H is a (deep) neural network called the decoder and Φ(·|η)η∈H 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 defineXi = (Xi,obs, Xi,mis) the partition of the data in realized observed and missing values given a specific realization of the pattern, it can be written as:
  `(θ) , n∑ i=1 log pθ (Xi,obs) = n∑ i=1 log ∫ pθ (Xi,obs|Zi) p(Zi)dZi.
  The corresponding evidence lower bound (ELBO) is:
  L(θ, γ) , n∑ i=1 EQγ [lnPθ (Xi,obs|Zi)]−KL (Qγ (Zi|Xi,obs) ‖Pθ (Zi)) ,
  with KL for the Kullback-Leibler divergence and the variational distribution
  Qγ (Z|Xobs) , Ψ (Z|gγ(Xobs)) ,
  with Ψ(·) the (parametric) variational distribution over Rd. The function gγ : X → K, 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 Qγ with Qγ (Z|Xobs) = Ψ (Z|gγ (ι (Xobs)) , where ι is an imputation function chosen beforehand that transforms Xobs into a complete input vector ι (Xobs) ∈ X .
sections.heading: B SUPPLEMENTARY RESULTS
sections.text: We only report the cases with logistic-linear specification of (e, µ0, µ1), we therefore use the logistic-linear doubly robust estimator and omit the prefix loglin in Figures 6 and 7. They show that when the number of covariates is large the MDC methods accurately recover the true ATE. In the DLVM case (bottom row), the MDC methods outperform all other compared methods.
  Additionally, we report results on the effect of changing the sampling parameter B for the MDC.mi strategy (all other parameters of the VAE are set to the default values). Again, we only report the cases with logistic-linear specification of (e, µ0, µ1). These results, in Figure 8, show that for small number of covariates (in this case, p = 5), the choice of B does not have a large influence on the final estimate, however for larger number of covariates (e.g., p = 100), increasing B reduces the bias and the variance.