Title: MISSDEEPCAUSAL: CAUSAL INFERENCE FROM IN-COMPLETE DATA USING DEEP LATENT VARIABLE MOD-ELS

Abstract: 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.

1. INTRODUCTION:

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.

2. SETTING, NOTATIONS AND RELATED WORKS:

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 (Y i (0), Y i (1), W i , X i ) i=1, ..., n with W i ∈ {0, 1} a binary treatment, X i = (X i1 , . . . , X ip ) ∈ R p a vector of covariates, and (Y i (0), Y i (1)) ∈ R 2 the outcomes we would have observed had we assigned control or treatment to the i-th sample, respectively. The observed outcome for unit i, Y i ∈ R is defined as

Y i W i Y i (1) + (1 -W i )Y i (0).

The individual causal effect of the treatment is τ i Y i (1) -Y i (0) and the average treatment effect (ATE) is defined as

τ E[Y i (1) -Y i (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:

EQUATION

where μw (x) are regression estimates of the conditional response surfaces

µ w (x) E[Y (w) | X =

x], w ∈ {0, 1}, and ê(x) is an estimate of the propensity score e(x) Rubin, 1983; Imbens & Rubin, 2015) .

P(W i = 1 | X i = x) (Rosenbaum &

Standard results state that if either (μ 0 , μ1 ) or ê 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:

{Y i (1), Y i (0)} ⊥ ⊥ W i | X i ,

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 M i ∈ {0, 1} p such that M ij = 0 if X ij is observed and M ij = 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 M i as a random vector, and the possibility to infer causal effects with missing data now depends on additional assumptions on the joint law of (Y i (0), Y i (1), W i , X i , M i ) 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

{Y i (1), Y i (0)} ⊥ ⊥ W i | 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:

EQUATION

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 given W 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

EQUATION

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.

3. MISSDEEPCAUSAL:

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:

X X M Z W Y {Y (0), Y (1)}

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

EQUATION

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[ĝ(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

EQUATION

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.

3.1. DEEP LATENT VARIABLE MODELS WITH MISSING VALUES:

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:

EQUATION

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.

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 (Z (j) ) 1≤j≤B ∈ R n×d from the posterior distribution P (Z|X ). On each sample, we evaluate the function f and aggregate the results: ĝ 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) :

(B) (X ) = 1 B B j=1 f (Z (

EQUATION

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 τ .

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 h(X ) defined in (7). For that purpose, we first approximate the expectation of the posterior distribution

EQUATION

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 Ẑ(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

ẽ(z) = P(W i = 1 | Ẑi (x ) = z),

and similarly for the conditional response surfaces.

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

4. SIMULATION STUDY:



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 4 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 Seaman & 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 .

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) . Missing values are generated completely at random (MCAR), i.e., P(M ij = 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.

Z i ∼ N d (0, 1), covariates X i are sampled from N p (µ (Z) , Σ (Z) ), where (µ (Z) , Σ (Z) ) = (V tanh(U Z + a) + b, diag{exp(η T tanh(U Z + a) + δ)}) with U, V,

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) . For this we define treatment and outcome models with a logistic-linear model as follows: logit(e

(Z i• )) = α T Z i• and Y i ∼ N ((β T Z i + τ W i , σ 2

). An estimation of τ is obtained by regressing the observed outcomes Y i on the estimations of the latent factors Z (for MDC.process, MF) and on the imputed data X imp (for MICE).

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). 

4.3.2. DOUBLY ROBUST ESTIMATION:

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.

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 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 2 More specifically, for estimating µw, we use all observations having Wi = w.

2 In 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) (X i ) and µ (0) (X i ) for all i):

∆ = τ -1 n µ (1) (X i ) -µ (0) (X i ) .

Table 1 shows that the doubly robust estimators (either in the parametric regression, DR log-lin , or the random forest form DR rf ) 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. 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.

5. CONCLUSION:

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.

A APPENDIX A.1 DEEP LATENT VARIABLE MODELS WITH MISSING VALUES:

Deep latent variable models can be defined as follows. Let (X i , Z i ) i≤n be n i.i.d. random variables such that

Z i ∼ P (Z i ) X i ∼ P θ (X i |Z i ) = Φ (X i |f θ (Z i )).

The prior distribution of the latent variables or codes

Z i ∈ R d is often isotropic Gaussian Z i ∼ N (0 d , I d ).

The function f θ : R d → 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 define X i = (X i,obs , X i,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 θ (X i,obs ) = n i=1 log p θ (X i,obs |Z i ) p(Z i )dZ i .

The corresponding evidence lower bound (ELBO) is:

L(θ, γ) n i=1 E Qγ [ln P θ (X i,obs |Z i )] -KL (Q γ (Z i |X i,obs ) P θ (Z i )) ,

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

Q γ (Z|X obs ) Ψ (Z|g γ (X obs )) ,

with Ψ(•) the (parametric) variational distribution over R d . 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|X obs ) = Ψ (Z|g γ (ι (X obs )) ,

where ι is an imputation function chosen beforehand that transforms X obs into a complete input vector ι (X obs ) ∈ X .

B SUPPLEMENTARY RESULTS:

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. 

:

https://github.com/imkemayer/causal-inference-missing

We use and adapt the corresponding code from V. Dorie: https://github.com/vdorie/npci/.

Figures and Tables:

Figure1: Unconfoundedness with missing values. X represents a the complete covariates, and M a missing data mechanism, X * represents the observed incomplete covariates, confounding the treatment assignment. The formalism ofPearl (1995) andRichardson & Robins (2013) is used.

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; X * represents the observed incomplete covariates.

a, b, δ, η drawn from Standard Gaussian distributions and Uniform distributions.

Figure3: Estimated ATE via regression adjustment; covariates generated from a DLVM, (logistic-)linear model specification for (e, µ 0 , µ 1 ); results with Z true results are obtained using the true confounders Z.

Figure 4: Estimated ATE for 10% of missing values; covariates are generated according to LRMF; Z true results are obtained using the true confounders Z.

Figure 6: Estimated ATE for different amounts of missing values with regression estimation; p = 100 covariates, n = 1, 000 observations, d = 3 latent confounders, logistic-linear model specification for (e, µ 0 , µ 1 ); top: LRMF; bottom: DLVM.

Figure 7: Estimated ATE for different amounts of missing values with doubly robust estimation; p = 100 covariates, n = 1, 000 observations, d = 3 latent confounders, logistic-linear model specification for (e, µ 0 , µ 1 ); top: LRMF; bottom: DLVM.

The results of MissDeepCausal are competitive with other 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).

Table content: <table><tr><td>% NA</td><td>Method</td><td>OLS</td><td>∆ DR log-lin</td><td>DR rf</td></tr><tr><td/><td>X (complete data)</td><td colspan="3">0.72 ± 0.02 0.13 ± 0.00 0.20 ± 0.01</td></tr><tr><td/><td>M F</td><td colspan="3">0.56 ± 0.03 0.14 ± 0.01 0.16 ± 0.01</td></tr><tr><td>0</td><td>M DC.process</td><td colspan="3">0.51 ± 0.03 0.15 ± 0.01 0.19 ± 0.03</td></tr><tr><td/><td>M DC.mi</td><td colspan="3">0.47 ± 0.03 0.16 ± 0.01 0.14 ± 0.02</td></tr><tr><td/><td>CEV AE(X)</td><td/><td>0.34 ± 0.02</td><td/></tr><tr><td/><td>M ICE</td><td colspan="3">0.85 ± 0.02 0.16 ± 0.00 0.24 ± 0.01</td></tr><tr><td/><td>M IA.GRF</td><td>-</td><td>-</td><td>0.23 ± 0.01</td></tr><tr><td>10</td><td>M F M DC.process</td><td colspan="3">0.50 ± 0.03 0.15 ± 0.01 0.15 ± 0.01 0.42 ± 0.02 0.15 ± 0.01 0.16 ± 0.02</td></tr><tr><td/><td>M DC.mi</td><td colspan="3">0.35 ± 0.02 0.17 ± 0.01 0.13 ± 0.02</td></tr><tr><td/><td>CEV AE(Xmean imp)</td><td/><td>0.31 ± 0.01</td><td/></tr><tr><td/><td>M ICE</td><td colspan="3">1.20 ± 0.02 0.30 ± 0.00 0.32 ± 0.01</td></tr><tr><td/><td>M IA.GRF</td><td>-</td><td>-</td><td>0.17 ± 0.01</td></tr><tr><td>30</td><td>M F M DC.process</td><td colspan="3">0.39 ± 0.02 0.16 ± 0.01 0.17 ± 0.01 0.37 ± 0.02 0.16 ± 0.01 0.15 ± 0.02</td></tr><tr><td/><td>M DC.mi</td><td colspan="3">0.30 ± 0.02 0.18 ± 0.01 0.13 ± 0.01</td></tr><tr><td/><td>CEV AE(Xmean imp)</td><td/><td>0.38 ± 0.02</td><td/></tr><tr><td/><td>M ICE</td><td colspan="3">1.54 ± 0.03 0.46 ± 0.01 0.42 ± 0.01</td></tr><tr><td/><td>M IA.GRF</td><td>-</td><td>-</td><td>0.19 ± 0.01</td></tr><tr><td>50</td><td>M F M DC.process</td><td colspan="3">0.28 ± 0.01 0.20 ± 0.01 0.21 ± 0.02 0.24 ± 0.01 0.20 ± 0.01 0.21 ± 0.02</td></tr><tr><td/><td>M DC.mi</td><td colspan="3">0.18 ± 0.01 0.22 ± 0.01 0.22 ± 0.03</td></tr><tr><td/><td>CEV AE(Xmean imp)</td><td/><td>0.38 ± 0.02</td><td/></tr></table>

Peer review (iclr2020):
Summary: 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. In this paper, the authors want to consider general non-linear relationships between Z and X with the same MAR (missing at random assumption) for missing entries. So they fit a latent variable model in the form of VAE to find the P(Z| observed part of X) using a slightly modified version of the ELBO Lower Bound. For missing entries, they just replace those entries by a constant and do the usual VAE fit. After the VAE fit, multiple Z*s are samples from the optimized P(Z| observed X) and then used in the doubly robust formula on each Z for estimating the average treatment effect and then finally the estimates averaged over the different Z*s. There is an alternative where the conditional mean of the latent variable is estimated from the VAE and used in the doubly robust computation. In many synthetic examples, authors compare this with existing baselines and show that their method improves. Pros: - Baselines compares are comprehensive enough from my perspective. - The paper is well written with clear pointer to existing work on doubly robust estimators with standard ignorability assumptions and the work for the linear, low rank model case by Kallus et al. 2018. Cons: Major Issues - There is no reason to believe that even for the synthetic experiments, that the VAE posteriors would asymptotically yield unbiased ATE*s which is provably the case in (Kallus et al. 2018) (of course for their restricted linear model/low rank assumptions). There is no reason to suppose Z*s used in eq (9) from the VAE satisfy ignorability even in the asymptotic limit. In this light, the paper in essence just estimates Z*s from some latent model that is fit and then use those latents to regress Y and then computes ATE. This seems a natural heuristic to try given such a problem. So I don*t find a big methodological novelty. I would be willing to increase my scores if the authors could convince me on this point. - For the LRMF model (Fig 4) the MF approach seems to do as well as the authors proposal and we have a guarantee for the MF case under those linear/low rank modeling assumptions. So the only demonstrated benefit is for the synthetic experiments for Fig 5 and Fig 3 (I agree that it is considerable particularly with large fraction of missing values in Fig 3) in whose settings we dont know about how unbiased it is in the limit. More synthetic experiments with different kind of generative models could be more convincing. - Some real world data set would have been more convincing - although I agree ground truth is hard to come by. Minor Issues: - You have set B=200 for all the experiments for MDC-MI. Do the results change when B is increased or decreased ?? Does variance go down or the bias itself changes with B - This would be a useful insight to have. - There is typo in the ELBO lower bound equation in page 11. There are other minor typos. Please correct for it. - Since the paper is about estimate treatment effect from measurements of an unobserved confounder - it is important to cite - https://ftp.cs.ucla.edu/pub/stat_ser/r366-reprint.pdf from the causal DAG literature. - The covariance of X given Z for the DLVM model is not clear - It seems to say exp ( or some matrix vector products) * Identity. What does this mean ? - The feature dimensions seems to be set at 10 - so would we expect the same results in much higher dimensions - like say 100s for few tens of thousands of samples?? ********UPDATE after reading the rebutall,changes and the new experiment********** I appreciate the authors actually accepting that identifiability issues cannot be easily resolved even if one knows P*(Z|X). I recommend the authors to elaborate on this point in the camera ready version. However, showing that the proposed methods work on a real benchmark (semi-synthetic one used in Shalit et. al 2017) is commendable. However, I find that the MF method is competitive (almost all the time) with their method when Doubly robust estimators are used. But having matched an existing baseline (the MF method) that deals with confounders on real data and showing superior synthetic results and authors clarifying and toning down their theoretical claims, I am inclined to increase the score to Weak accept.