Causal assessment often relies on the framework of counterfactual inference. In this framework, each unit, i, has a potential outcome given that they received a treatment and a potential outcome given that they received a control – and , respectively. This framework seeks to contrast the outcome, Y for an individual under these two hypothetical states as shown in Eq. 1 (Rubin 1974).

The effect of the treatment on the outcome can then summarized by calculating population-level effect estimates, such as the average treatment effect (ATE), which is defined as the expected difference in outcomes (Eq. 2).

Estimating this requires access to the outcome for the state in which units were not assigned (i.e., and ). In practice, however, these true counterfactuals are never observed as a single population (or individual) cannot simultaneously be both treated and untreated. This is known as the ’fundamental problem of causal inference.’ Therefore, approximations that employ more than one population are used as a proxy for these unobserved states (Holland 1986). These approximations seek to construct populations such that the observed ATE, , equals the true ATE that would arise from a counterfactual population. In other words, we seek an

A decomposition of the ATE, demonstrates that a sufficient condition for unbiased estimation is that and (Kempthorne 1955). Within the counterfactual framework, this equality is central to the assumption of strong ignorability (Eq. 4) (Rosenbaum and Rubin 1983).

This assumption states a unit’s assignment to a treatment is independent of that unit’s potential outcomes, , and that treatment assignment is, therefore, ignorable. Causal claims borne from data that satisfy this requirement are regarded as unbiased as all confounding factors that could induce a dependence between are equally represented in the treatment and comparator arms (Rubin 1974). Consequently, this means that the distribution of features is the same in both arms and features are said to be balanced. Other assumptions, such as positivity and the Stable Unit Treatment Value Assumption (SUTVA), are also necessary and assumed to be true (Rubin 1980).

Matching and weighting are popular pre-analysis manipulations to approximate the unconditional form of strong ignorability in observational populations. These methods create pseudo-populations in which the assumption is met without need for further manipulation (Rubin 1973). This is opposed to methods of statistical adjustment, which occur peri-analysis, and approximate the conditional form of strong ignorability (Leger 1994). Arguably, the most common strategy for weighting is the inverse probability of treatment weighting (IPW) (Thoemmes and Ong 2016), though other methods include the direct minimization of imbalance (Gret- ton et al. 2009; Kallus 2016; Kallus 2017) or weighting by the odds of treatment, kernel weighting, and overlap weighting (Rosenblatt 1956; Hellerstein and Imbens 1999; Hazlett 2016; Li, Morgan, and Zaslavsky 2018; Kallus 2018b).

A commonality among these methods is that they implicitly or explicitly all specify a distribution function, P, that the expectation in Eq. 2 is taken with respect to. This distribution is often the distribution associated with the treated (), the controls (), or a combination thereof (e.g.

). This choice of distribution can lead to high variance effect estimates in circumstances where there are regions of poor overlapping support between the treated and untreated populations. An effect of this is often observed in the context of IPW analyses with instability due to propensity scores near zero or one. (Kang and Schafer 2007).

In this work, we instead construct an implicit distribution, P, that focuses on the regions of the sample space with significant overlap between the treated and untreated populations. Such a construction involves an inherent trade-off between coverage and variance. For example, mixture distributions that will be valid for a larger region of the sample space will also produce high variance estimates in the context of a fixed sample budget. In the context of infinite sample sizes and positivity, one could specify any distribution P without concern for effect estimate variance. The mixture distribution of the treated and untreated populations would be a reasonable choice given a goal of maximizing coverage. However, in real-world settings with limited data, positivity may not be present and ATE estimates over such a distribution may be high variance in practice and theoretically invalid. In such a setting, valid estimates can only be made for subpopulations with significant distributional overlap. We formulate an approach that constructs a distribution P for estimating Eq. 2 that both maximizes coverage and minimizes variance. Informally, P can be considered the distribution of a natural experiment where the choice of treatment, T, is independent of potential confounders, X.

We propose the Counterfactual -GAN (cGAN) that uses an adversarial approach to learn a distribution that trades off coverage and effect estimate variance for two or more observational study arms. This approach learns stable, feature balancing weights without reliance on the propensity score. The target distribution, P, is identified by minimizing the Pearson -divergence between P and the sampling distributions for each study arm. To our knowledge, this is the first such application of the Pearson -divergence. Because P is being compared to all study arms, this encourages coverage, while, as we will show, the -divergence inherently minimizes the variance of importance sampling estimates of the ATE.

This paper proceeds as follows: Section 2 defines the model and learning procedures, Section 3 presents an evaluation of this model through a simulation and an application to real-world clinical data, and finally, Section 4 discusses open issues, limitations of the model, and future work.

We introduce the Counterfactual (cGAN), an adversarial approach to feature balance in causal inference that is based on importance sampling theory. Using an adversarial approach based on variational minimization based on the f-GAN, we minimize the sum of the Pearson -divergences between a deep generative model and the sampling distributions from each arm of a study. We show that minimizing the -divergence is equivalent, up to a constant factor, to minimizing the variance of importance sampling estimates to be made in approximating quantities such as ATEs. Similar to other weighting approaches, this approach assumes SUTVA, positivity, and no unmeasured confounders. In the following, P is the constructed target distribution and the sampling distribution for each study arm.

Figure 1: Architecture of Counterfactual

Importance Sampling and the -divergence Importance sampling is a strategy for estimating expectations under an unknown target distribution given a known proposal distribution (Muller 1966). Though the importance sampling has broader usage than our application, we focused on the use of importance sampling for estimation of the average treatment effect (ATE) because of its close relationship with the divergence. The importance sampling weight is de-fined as a likelihood ratio: the likelihood of an observation under the target distribution, p(x) divided by the likelihood under the proposal distribution, q(x). Weighted expectations based on the proposal distribution approximate unweighted expectations from the target distribution at shown in Eq. 5.

Consider the units in an arm of an observational study as being samples from such a proposal distribution. One strategy for obtaining unbiased expectations of treatment effects is to identify importance sampling weights for each arm that approximate expectations from a shared target distribution. However, this problem is underspecified given that we could choose any target distribution with the correct support. In this work, we choose the target distribution that yields importance sampling approximations with smallest variance. Eq. 6 shows the form for the variance of importance sampling estimates where is the constant function. This choice is to make the formulation of the cGAN as outcome agnostic as possible. This form highlights its connection with the -divergence, which has a function form as shown in Eq. 7. This connection was previously noted in (Dieng et al. 2017). Therefore, the solution which minimizes the -divergence would also minimize the variance expectations for unknown outcomes. Of note, importance sampling is known to be a method that can produce high variance estimates, but since we will be minimizing the variance directly, this is less of a concern here.

Likelihood Ratio ATE Estimation Typically, the expectation in the ATE is taken with respect to the original feature distribution, q(x). Under cGAN-weighted data, expectations are taken with respect to the target distribution p(x). As such, calculations of the ATE from the cGAN are not equivalent to what many would classically consider the ATE, but rather, is an ATE with respect to the new, learned feature distribution. We call this new estimate the . This inequality is demonstrated in Equation 8. This set of equations shows that the typical ATE, , is not equivalent to the expectation that we estimate, the

Consider two distributions and that represent two arms of a study. It is possible to make unbiased estimates based on a single distribution, P, leveraging likelihood ratios/importance sampling weights as shown in Eq. 9.

(9) We will leverage an approach based on adversarial learning to simultaneously maximizes coverage, minimizes the variance defined in Eq. 6, and directly estimates likelihood ratios,

f-GAN The f-GAN framework provides a strategy for estimation and minimization of arbitrary f-divergences based on a variational divergence minimization approach (Nowozin, Cseke, and Tomioka 2016).

where T is a class of function such that is the function that characterizes the -divergence, f(u) = is the Fenchel conjugate of and P and Q are probability distributions with continuous densities, p(x) and q(x). T is typically a multi-layer neural network. This formulation lower bounds the -divergence based on functions T, P, and Q in such a way that unbiased noisy gradients of the lower bound can be easily obtained based on samples from P and Q. In addition, the variational function, T, has a tight bound for which is

equivalent to in the case of the -divergence. To respect the bounds of T that result in valid likelihood ratios, we represent T as a nonlinear transformation of an unbounded function . The likelihood ratio, , is easily derived from here and provides the importance sampling weights necessary for approximating expectations under p(x) as shown in Eq. 5.

The Counterfactual The cGAN builds on importance sampling theory and extends the f-GAN framework to learn feature balancing weights through an adversarial training process. Previously, (Tao et al. 2018) have explored importance weights from critics of divergence-based GAN models. However, unlike this method and other f-GANs where there is a generator, G and a single variational function, the cGAN employs dual training from at least two variational functions (Figure 1).

Consider a set of A treatments, each associated with one of A populations, or arms of a study. Each population contains units and are drawn from an unknown and population-specific distribution . Based on the connection between the -divergence and the variance of importance sampling estimates outlined above, our objective is to identify a target distribution that minimizes the -divergence to all populations being compared: . This is the sum of the divergences between the generator and the unweighted treatment arms. It is minimized when p(x) equals and is directly proportional to the sum of the variances of importance sampling estimates under the target distribution, P, with proposals, . Because of the constant in Eq. 6, minimizing the -divergence is equivalent to minimizing a normalized variance which weighs each population equally regardless of the number of units and the magnitude of the treatment effect,

Algorithm 1: Minibatch stochastic gradient cGAN optimization

As a byproduct of minimizing this divergence, we will also identify a set of importance weights, , for each unit in each population that allows estimation of expectations from the same target distribution, P, thus satisfying the unconditional form of strong ignorability. Using these importance weights, expectations can be approximated as where , where

is an normalizing constant, p is the den- sity of the shared target distribution, is the density of the proposal distribution, and . Note that our strategy eliminates the need to explicitly evaluate and as the likelihood ratio is estimated directly by the f-GAN. If desired, expectations can also be approximated using the sample-importance-resampling (SIR) algorithm where samples approximately distributed according to p can be simulated by drawing samples from the weighted empirical distribution Freitas, and Gordon 2001).

The objective function for the cGAN is shown in Eq. 11 and is closely related to the objective defined in (Nowozin, Cseke, and Tomioka 2016). parameterizes the generative model and parameterizes the variational model for each treatment arm, a. In our experiments, are neural networks that mirror discriminators in the traditional GAN framework and is a neural networks that mirrors the generator. Note that the generator in the original f-GAN framework is usually . In our case, to achieve the desired directionality of the -divergence, the empirical distribution

must be and the generator must be P.

Importance weights can be computed based on the fact that the bound in Eq. 10 is tight for

where . We can therefore, approximate the desired importance weights as described in Eq. 5 as for all and . Ultimately, the ATE can be estimated between any two treatment arms according to Eq. 9. For example, the ATE between arms 1 and 2 could be estimated as

Practical Considerations In the original GAN and f-GAN formulations the gradients for the generator is replaced with a related gradient that significantly speeds convergence of the model. Because our objective is minimization of the true -divergence rather than perfect distributional matching, we do not employ this loss function trick but instead apply the gradient as derived from the loss function in Eq. 11.

Although it is the case that the domain of the Fenchel conjugate for the -divergence is R, we constrained it to which produces valid likelihood ratios.

Gradient descent-based optimization of GANs is a notedly difficult task (Mescheder, Geiger, and Nowozin 2018; Arjovsky and Bottou 2017; Gulrajani et al. 2017). Though many methods are proposed to stabilize training, we have found it sufficient to employ a set of algorithmic heuristics: (i) standardization of our data by the joint mean and variance over all A populations prior to training; (ii) periodically recentering the distribution of each discriminator to a noisy estimate of the mean of the generator distribution. This recentering is accomplished by setting the value of a vector that is added to the input of the discriminators.

The approach for minibatch stochastic gradient descent for the cGAN is shown in Algorithm 1. The objective function F (Eq. 11) is optimized by minimizing with respect to the parameters of the generator and maximizing with respect to the parameters of the discriminators.

Related Work Causal inference with observational data has a rich literature that cuts across many disciplines (Thrus- field 2017; Rubin 1973; Rubin 1974; Pearl 2000) including machine learning (Johansson, Shalit, and Sontag 2016; Kallus 2018a; Shalit, Johansson, and Sontag 2017; Ratkovic 2014; Schwab, Linhardt, and Karlen 2018). More specifi-cally there have been several approaches to applying adversarial networks for counterfactual inference (Kallus 2018a; Yoon, Jordon, and van der Schaar 2018). However, most existing methods for counterfactual inference are not directly comparable to the cGAN, as we aim to identify the most appropriate counterfactual distribution given the available data and maximize feature balance whereas most methods evaluate ATE estimation or ITE estimation directly.

In contrast to representational learning approaches and some GAN approaches, our approach does not rely on a pre-defined outcome to identify matched cohorts. The approach outlined in (Kallus 2018a) is the most similar in spirit to our approach but differs in that our objective directly minimizes the variance of expectations that might be used in ATE estimation, whereas (Kallus 2018a) minimizes a bound on the variance of the average treatment effect on the treated. As a result, there is no need for a regularizer, to perform crossvalidation to select an appropriate level of regularization, or perform a constrained optimization over weights.

To evaluate the cGAN, including its utility in practice, we present results of a simulation and applications to real-world medical data.

Simulation

To evaluate the cGAN when the ground truth is known, we applied the model on simulated data of two populations/treatment arms, A = 2. Each population was comprised of two subpopulations. Each subpopulation contained 10 features, drawn from a randomly generated multivariate normal distribution with a normal-Wishart prior distribution. Population 1 was composed of an equal number of samples (N=1000) from subpopulation A and subpopulation B; and Population 2 was composed of an equal number of samples from subpopulation A and subpopulation C (N=2000). By construction, subpopulation A is a latent population associated with a natural experiment, since it is part of both Population 1 and 2.

Because our simulation deliberately constructs populations from a shared subpopulation distribution (A), we would expect points generated from this subpopulation to have higher weights. Intuitively, the variance of importance sampling estimates should be small for both treatment groups (a = 1 and a = 2) if the learned target distribution, is one that overlaps both populations maximally while excluding density unique to one group.

To better demonstrate how the cGAN supports counterfactual reasoning, we have additionally conducted an analysis of the average treatment effect (ATE) for our experiment with simulated data. We simulated a continuous outcome according to the subpopulation of origin – Pop 1A Gaussian (60, 1); Pop 1B Gaussian (40, 1); Pop 2A Gaussian (-10, 1); Pop 2C Gaussian (10, 1). Under this outcome function, the estimate of average treatment effect (ATE) under the mixture distribution (of Pop 1 and Pop 2) is 50. When estimating the ATE under the overlapping subpopulation distribution – those from Pop 1A and Pop 2A – the ATE is 70. We applied weights from the cGAN and comparators to the simulated outcomes to assess the ability of the weighting methods to estimate one of the two ATEs. In addition, we also calculated the effective sample size (ESS), , using the Kish Method (Kish 1965). The ESS may be used to determine the quality of a Monte Carlo approximations of importance sampling. The calculation of can be found in the equation below, wherein w are the weights.

ness of ATE, and (iii) the ESS, a variety of comparator methods were implemented in addition to the cGAN . They include binary regression propensity score; generalized boosted modeling of propensity scores (McCaffrey, Ridgeway, and Mor- ral 2004); covariate-balancing propensity scores (Imai and Ratkovic 2013); non-parametric covariate-balancing propensity scores (Fong, Hazlett, and Imai 2018); entropy balancing weights (Hainmueller 2011); empirical balancing calibration weights (Chan, Yam, and Zhang 2016); optimization-based weights (Keele and Zubizarreta 2014).

Results The results of our simulation is summarized in Figures 2. In the left hand-side of the Figure, the columns show the marginals of three pairs of continuous features. Row (i) shows the raw data, colored by which population units were drawn from. Row (ii) shows the same data as above, but coloring by subpopulation to highlight the overlapping distribution. Row (iii) shows a set of samples from the generator after training colored in blue. Row (iv) depicts the original data from Row (i) with the opacity of data points reflecting the importance weights. The right-hand side of the Figure shows the distribution of weights by subpopulation. Note that, in both Populations 1 and 2, the mean weights of units from subpopulation A have weights near , which is the uniform weight when 2000 units are in each population. Units from other subpopulations have near negligible weights, and would not meaningfully contribute to expectations in 9.

In the left-most figure, as you move down any column of feature pairs, it is apparent that points from the overlapping subpopulation A are both captured by the generator and assigned higher weights. This is confirmed by plotting the weights of data points by subpopulation (right-hand side of 2). Weights from subpopulations 1A and 2A are substantially higher than those from subpopulations 1B and 2C.

The results of this simulation further demonstrate that the ATE estimate from cGAN-weighted data is less biased than estimates from other weighting methods, given their respective targets. By construction, the causal effect of the comparable subpopulations is 70. cGAN-weighted data produced an ATE of 70.01. We see similarly good performance when inspecting the ESS. The cGAN has an ESS of 3870.

Figure 2: Simulation Results. Left: Select features (i) by population of origin; (ii) with subpopulation A highlighted; (iii) samples from the generator; (iv) opacity adjusted by weight. Right: Weights by subpopulation

Given that there are 4000 units that are comparable across the two arms (each subpopulation contains 2000 units), this is an appropriate estimate (Table 2).

unweighted 50.03 8000 IPW 92.00 6551 clipped IPW 87.24 6997 binary regression PS 92.00 6551 generalized boosted modeling PS 84.51 7207 covariate balancing PS 91.83 6686 non-parametric covariate balancing PS 37.65 11

optimization-based weights 52.07 114 cGAN 70.01 3870

Table 2: Results of Simulation. The average treatment effect and effective sample size (ESS) after application of weighting methods from the Counterfactual -GAN and comparators.

Clinical Data

We additionally applied the cGAN to an experiments using real-world clinical data from a large, academic medical center. For this experiment, we constructed the treatment and comparator cohorts according to the protocol and indication of a published randomized clinical trial. The experiment compares sitagliptin and glimepiride in elderly patients with Type II Diabetes Mellitus (N=144 per arm) (Hartley et al. 2015). We present the 37 most frequent clinical measurements from the electronic health record.

We evaluate the ability of the cGAN to improve feature balance by comparing the Absolute Standardized Difference of Means (ASDM) between the treatment and comparator cohorts under different weighting methods. the ASDM is a popular method of assessing cohort similarity, with a lower metric corresponding to improved feature balance. The ASDM is presented for the cGAN and the comparator weighting methods mentioned in the simulation. Under the clipped-IPW procedure, propensity scores greater than 90th percentile and less than 10th percentile are assigned to the values of the percentiles at 90th and 10th, respectively (Cole and Hern´an 2008).

Results The ASDM for the clinical cohorts is presented in Figure 3. These findings are summarized by the mean ASDM over all features, under the varying weighting methods in Table 3. cGAN improved mean ASDM from the unweighted cohort and improved feature balance the most among all evaluated methods. Note that this task is particularly challenging due to the high dimensionality of the data and small study size.

The results of this experiment can be found in Figure 3 and Table 3. They demonstrate that cGAN-weighting achieves better feature balance than comparator methods.

Figure 3: Absolute standardized difference of the means (ASDM) of real-world clinical features after application weighting methods from the Counterfactual -GAN and comparators.

In this paper, we introduce the Counterfactual -GAN. It is a deep generative model for feature balance that minimizes the variance of importance sampling estimates of treatment effects. We leverage the f-GAN framework for estimating the -divergence and likelihood ratios necessary for achieving this. The experiments presented here suggest that cGAN is an effective method of learning feature balancing weights to support counterfactual inference. If we assume that all potentially confounding variables are observed, the superiority of cGAN in learning balancing weights, suggests that ATE borne from cGAN-weighted cohorts would be less biased than those estimates generated from traditional weighting methods. The application of the model to real-world EHR data, demonstrates that this method could provide an alternative means to causal estimation from observational data when the assumptions of no unobserved confounding, positivity, and SUTVA are met. Our experiments suggest that the flexibility

of our framework produces improved feature balance relevant for valid causal estimates. This method does, however, come with limitations. Training of the model is completed via backpropogation. As such, this method is only suitable for fully differentiable functions. Therefore, matching based on a combination of discrete and continuous values poses a challenge. In addition, GANs are well known for their instability and lack of objective measures for convergence. This work shares those limitations. In future work, we will explore an extension of the cGAN which accommodates discrete data types and overcome the many current limitations of GANs.

Acknowledgments This research is supported by grants R01LM009886-10 and T15LM007079 from The National Library of Medicine.

[Arjovsky and Bottou 2017] Arjovsky, M., and Bottou, L. 2017. Towards principled methods for training generative adversarial networks.

unweighted 0.1103 IPW 0.0876 clipped IPW 0.0631 binary regression PS 0.0625 generalized boosted modeling PS 0.0749 covariate balancing PS 0.0681 non-parametric covariate balancing PS 0.0596 entropy balancing 0.0524 empirical balancing calibration weights 0.0524 optimization-based weights 0.0536 cGAN 0.0364

Table 3: Results of Application to Clinical Data. Absolute standardized difference of the means (ASDM) of real-world clinical features after application weighting methods from the Counterfactual -GAN and comparators.

[Chan, Yam, and Zhang 2016] Chan, K. C. G.; Yam, S. C. P.; and Zhang, Z. 2016. Globally Efficient Nonparametric Inference of Average Treatment Effects by Empirical Balancing Calibration Weighting. Technical report, University of Washington.

[Cole and Hern´an 2008] Cole, S. R., and Hern´an, M. A. 2008. Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 168(6):656– 64.

[Dieng et al. 2017] Dieng, A. B.; Tran, D.; Ranganath, R.; Paisley, J.; and Blei, D. M. 2017. Variational Inference via Upper Bound Minimization. In NIPS.

[Doucet, Freitas, and Gordon 2001] Doucet, A.; Freitas, N.; and Gordon, N. 2001. An Introduction to Sequential Monte Carlo Methods. In Sequential Monte Carlo Methods in Practice. New York, NY: Springer New York. 3–14.

[Fong, Hazlett, and Imai 2018] Fong, C.; Hazlett, C.; and Imai, K. 2018. Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements. The Annals of Applied Statistics 12(1):156– 177.

[Gretton et al. 2009] Gretton, A.; Smola, A.; Huang, J.; Schmittfull, M.; Borgwardt, K.; Sch¨olkopf, B.; Candela, J.; Sugiyama, M.; Schwaighofer, A.; and Lawrence, N. 2009. Covariate shift by kernel mean matching. Dataset Shift in Machine Learning, 131-160 (2009).

[Gulrajani et al. 2017] Gulrajani, I.; Ahmed, F.; Arjovsky, M.; Dumoulin, V.; and Courville, A. 2017. Improved Training of Wasserstein GANs. In NIPS.

[Hainmueller 2011] Hainmueller, J. 2011. Entropy Balancing for Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in Observational Studies. Political Analysis 16:25–46.

[Hartley et al. 2015] Hartley, P.; Shentu, Y.; Betz-Schiff, P.; Golm, G. T.; Sisk, C. M.; Engel, S. S.; and Shankar, R. R. 2015. Efficacy and Tolerability of Sitagliptin Compared with Glimepiride in Elderly Patients with Type 2 Diabetes Mellitus

and Inadequate Glycemic Control: A Randomized, DoubleBlind, Non-Inferiority Trial. Drugs & Aging 32(6):469–476.

[Hazlett 2016] Hazlett, C. 2016. Kernel Balancing: A Flexible Non-Parametric Weighting Procedure for Estimating Causal Effects. SSRN.

[Hellerstein and Imbens 1999] Hellerstein, J. K., and Imbens, G. W. 1999. Imposing moment restrictions from auxiliary data by weighting. Review of Economics and Statistics 81:1– 14.

[Holland 1986] Holland, P. W. 1986. Statistics and Causal Inference. JASA 81(396):945–960.

[Imai and Ratkovic 2013] Imai, K., and Ratkovic, M. 2013. Covariate balancing propensity score. Technical report, Harvard University.

[Johansson, Shalit, and Sontag 2016] Johansson, F. D.; Shalit, U.; and Sontag, D. 2016. Learning representations for counterfactual inference.

[Kallus 2016] Kallus, N. 2016. Causal inference by minimizing the dual norm of bias: Kernel matching & weighting estimators for causal effects. CEUR Workshop Proceedings 1792:18–28.

[Kallus 2017] Kallus, N. 2017. A Framework for Optimal Matching for Causal Inference. In AISTATS, 372–381.

[Kallus 2018a] Kallus, N. 2018a. Deepmatch: Balancing deep covariate representations for causal inference using adversarial training.

[Kallus 2018b] Kallus, N. 2018b. Optimal a priori balance in the design of controlled experiments. J. R. Statist. Soc. B 80(1):85–112.

[Kang and Schafer 2007] Kang, J. D. Y., and Schafer, J. L. 2007. Demystifying Double Robustness: A Comparison of Alternative Strategies for Estimating a Population Mean from Incomplete Data. Statistical Science 22(4):523–539.

[Keele and Zubizarreta 2014] Keele, L., and Zubizarreta, J. 2014. Optimal Multilevel Matching in Clustered Observational Studies : A Case Study of the School Voucher System in Chile. arXiv preprint arXiv:1409.8597 1–37.

[Kempthorne 1955] Kempthorne, O. 1955. The randomizaton theory of experimental inference. Journal of American Statistics 50:946–967.

[Kish 1965] Kish, L. 1965. Survey Sampling. In Survey Sampling. New York: Wiley. chapter Chapter 14.

[Leger 1994] Leger, A. S. S. 1994. Statistical Models in Epidemiology, volume 48. Oxford University Press.

[Li, Morgan, and Zaslavsky 2018] Li, F.; Morgan, K. L.; and Zaslavsky, A. M. 2018. Balancing Covariates via Propensity Score Weighting. Journal of the American Statistical Association 113(521):390–400.

[McCaffrey, Ridgeway, and Morral 2004] McCaffrey, D. F.; Ridgeway, G.; and Morral, A. R. 2004. Propensity Score Estimation With Boosted Regression for Evaluating Causal Effects in Observational Studies. Psychological Methods 9(4):403–425.

[Mescheder, Geiger, and Nowozin 2018] Mescheder, L.;

Geiger, A.; and Nowozin, S. 2018. Which training methods for gans do actually converge?

[Muller 1966] Muller, M. E. 1966. Review: J. M. Hammersley, D. C. Handscomb, Monte Carlo Methods ; Yu. A. Shreider, Methods of Statistical Testing/Monte Carlo Method, volume 37. Springer Netherlands.

[Nowozin, Cseke, and Tomioka 2016] Nowozin, S.; Cseke, B.; and Tomioka, R. 2016. f-GAN: Training Generative Neural Samplers using Variational Divergence Minimization. In NIPS, 271–279.

[Pearl 2000] Pearl, J. 2000. Causality. Cambridge, England: Cambridge University Press.

[Ratkovic 2014] Ratkovic, M. 2014. Balancing within the margin: Causal effect estimation with support vector machines. Department of Politics, Princeton University, Princeton, NJ.

[Rosenbaum and Rubin 1983] Rosenbaum, P. R., and Rubin, D. B. 1983. Assessing Sensitivity to an Unobserved Binary Covariate in an Observational Study with Binary Outcome.

[Rosenblatt 1956] Rosenblatt, M. 1956. Remarks on Some Nonparametric Estimates of a Density Function. The Annals of Mathematical Statistics 27(3):832–837.

[Rubin 1973] Rubin, D. B. 1973. The Use of Matched Sampling and Regression Adjustment to Remove Bias in Observational Studies. Biometrics 29(1):185.

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

[Rubin 1980] Rubin, D. 1980. Randomization analysis of experimental data: The Fisher randomization test comment. JASA 75(371):591–593.

[Schwab, Linhardt, and Karlen 2018] Schwab, P.; Linhardt, L.; and Karlen, W. 2018. Perfect match: A simple method for learning representations for counterfactual inference with neural networks.

[Shalit, Johansson, and Sontag 2017] Shalit, U.; Johansson, F. D.; and Sontag, D. 2017. Estimating individual treatment effect: generalization bounds and algorithms. In ICML.

[Tao et al. 2018] Tao, C.; Chen, L.; Henao, R.; Feng, J.; and Duke, L. C. 2018. Chi-square Generative Adversarial Network.

[Thoemmes and Ong 2016] Thoemmes, F., and Ong, A. D. 2016. A Primer on Inverse Probability of Treatment Weighting and Marginal Structural Models. Emerging Adulthood 4(1):40–59.

[Thrusfield 2017] Thrusfield, M. 2017. Observational studies. In Veterinary Epidemiology: Fourth Edition. Springer. 319– 338.

[Yoon, Jordon, and van der Schaar 2018] Yoon, J.; Jordon, J.; and van der Schaar, M. 2018. GANITE: estimation of individualized treatment effects using generative adversarial nets. In ICLR.