Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
85,651 characters · 17 sections · 37 citation commands
A Convexified Matching Approach to Imputation and Individualized Inference
Keywords--- Missing value imputation, optimal coupling, matching, propensity score, individualization, entropic regularization.
One central topic in applied econometric research is to assess the effects of policy interventions reliably. On the one hand, to analyze nonexperimental data, advanced econometric estimates are devised to provide granular counterfactual answers by leveraging specific structural models. In an influential paper in 1986 lalonde1986evaluating, LaLonde questioned whether such sophisticated econometric estimates are credible. He compared them with the experimental benchmarks on some coarse summary---for example, the average treatment effect (ATE)---and concluded unfavorably. On the other hand, when experimental data, as in randomized controlled trials, are available, elementary statistical estimates can determine the ATE and eliminate confounding explanations. The ATE estimate, however, does not answer whether the treatment works for an individual. In modern applications such as personalized medicine and online marketing, the treatment effects vary across individuals; the treatment might be beneficial for some individuals but ineffective for others liu2016there. One potentially costly approach is to conduct individualized experiments in different time windows, known as N-of-1 trials hill1961principles,liang2023randomization. Therefore, it is desirable to develop nonexperimental methods that conform to the coarse estimates---for example, those for ATE or the average treatment effect on the treated (ATT)---at the aggregate level, while delivering more granular, individualized inference.
Individualized inference is about imputing missing counterfactual outcomes---the outcomes that would have been observed had the subjects received an opposite treatment---and their uncertainty quantifications. Three popular approaches are available in the literature: matching, regression imputation, and synthetic control. This paper proposes a new convexified matching method by integrating these three approaches. The central quantity we play is a coupling matrix between the treated and control data clouds, which resembles matching and is used to synthesize the counterfactual outcomes as in classic nonparametric regression.
Matching rubin2006matched,rosenbaum2020modern is a widely adopted way to impute counterfactual outcomes. For each subject, the counterfactual outcome is estimated by identifying units in the opposite treatment group similar to the subject in covariates and then averaging their outcomes. Due to its simplicity, nearest neighbor matching, which finds subjects closest under a suitable distance between the covariates, is favored in practice. Exact or approximate matching with high-dimensional covariates can be difficult, and thus propensity score matching was proposed as a remedy rosenbaum1983central. There are more sophisticated matching methods to improve balance or forbid certain matches, often referred to as the optimal matching rosenbaum2020modern,zubizarreta2023handbook; they require solving combinatorial optimization problems by network optimization techniques or mixed integer programming, which can be computationally expensive. However, even analyzing the aggregate ATE estimates resulting from combinatorial optimization problems is highly nontrivial, and these estimates may differ from simple ATE estimates based on propensity score reweighting rosenbaum1987model,hirano2001estimation,hirano2003efficient. As in matching, we solve for a coupling matrix---a convex relaxation to the combinatorial constraints---without access to outcomes, separating the design and analysis phases in an observational study. For a large-scale matching problem with $N$ units, we devise entropic-regularized algorithms from optimal transport villani_2003,peyre2019computational to solve convexified matching, calling a matrix scaling subroutine with a number of times nearly independent of $N$; see Section (ref). In contrast, in a conventional matching problem, solvability within polynomial time of $N$ may not be feasible.
Regression imputation is another approach to estimating counterfactual outcomes. The idea is simple: under the ignorability or unconfoundedness assumption rosenbaum1983central, the underlying functions that map the covariates to the responses under the treatment and the control are identifiable and estimable based on observational data. It translates a causal inference problem into a regression problem. Parametric or nonparametric regression techniques estimate the conditional response function under the treatment and the control, given the covariates, and in turn find the conditional average treatment effect (CATE), see hahn1998role,heckman1997matching,heckman2005structural,macurdy2011flexible. At the aggregate level, semiparametric efficient estimation of ATE leveraging regression imputation techniques is studied in hahn1998role. Admittedly, if the quantity of interest is at the aggregate-level---a one-dimensional parameter such as ATE or ATT---the specific nonparametric regression technique matters less, as long as it estimates the function reasonably well. There has been a fruitful line of literature on the use of flexible machine learning methods in estimating CATE belloni2014inference,athey2016recursive,chernozhukov2017double,farrell2015robust,farrell2021deep,wager2018estimation, in place of classic Nadaraya-Watson (NW) kernel nonparametric regression. Albeit naive, let us use the NW estimator as an example to illustrate a shared feature by regression imputation and matching: to compute a counterfactual outcome for a treated unit, NW finds local neighbors in the control group and uses convex weights based on the nonparametric kernel to synthesize the outcome. We will adopt this aggregation by convex weights feature into our convexified matching, whereas our coupling weights are optimized globally. Our coupling weights are determined to optimize imputation quality and individual uncertainty quantification; see Section (ref).
Hybrid methods were proposed to combine matching/propensity score weighting rosenbaum1987model,hirano2001estimation,hirano2003efficient,li2018balancing with regression imputation techniques to efficiently estimate ATE, notably the doubly robust estimators laan2003unified,cattaneo2010efficient,farrell2015robust,farrell2021deep,chernozhukov2022locally. The estimator compares favorably as it is valid when either the propensity score function or the regression function is consistent, robust to bias due to misspecification. Similar in spirit, bias correction using regression to improve the nearest neighbor matching was studied in abadie2011bias,abadie2006large.
From a different vein, synthetic control abadie2003economic,abadie2010synthetic,abadie2021using provides a fresh look at counterfactual imputation: a convex combination of control units may synthesize a subject under treatment better than any one of the control unit. Mathematically, a convex combination of data may simultaneously reduce the approximation error or bias, and, at the same time, reduce the variance driven by the idiosyncratic errors. This convex relaxation idea also eases the computation for the optimization problem in matching. We adopt this convex combination idea to approximate the overall covariate information under treatment by convex combinations of the covariate information under control, defined based on a data set-to-data set coupling; see Section (ref).
We develop a convexified matching method that solves an optimal coupling, which in turn defines convex combinations for missing value imputation. Unlike combinatorial optimization problems, where the optimization variable must be an extreme point of a certain polytope to represent an assignment, we optimize over the whole polytope as the variable represents the weights of the convex combination. Moreover, we add an entropic regularization to the optimization, where the strength of regularization controls the bias and variance tradeoff. The resulting formulation is a smooth convex optimization problem confined to the polytope that can be solved efficiently. Lastly, by properly specifying the constraints, we can pair with propensity score weighting rosenbaum1987model,hirano2001estimation,hirano2003efficient estimators. We can guarantee that the aggregate summary of individual treatment effects coincides with the desired ATE or ATT estimates.
We introduce an inference procedure to construct individual confidence intervals for the estimated counterfactual outcomes based on the convexified matching method. Providing a credible confidence interval around the individual treatment effect can help decide whether to adopt the treatment for each individual. The width of the confidence interval is determined by the approximation error in covariate balancing and the entropy of the convex weights, which comprise the objective function of the optimization. Therefore, how we formulate the convexified matching is directly targeted at individualized inference. Entropic regularization plays a crucial role in inference, namely, controlling the width of the individual confidence intervals, and in computation, namely, designing fast iterative algorithms to solve the optimization.
In summary, we integrate important features from matching, regression imputation, and synthetic control. We impute counterfactual outcomes by convex combinations defined based on an optimal coupling. The coupling is a convex relaxation to optimal matching and can be solved efficiently using iterative matrix scaling subroutines called the Sinkhorn algorithm sinkhorn1967diagonal,cuturi2013sinkhorn. We estimate granular individual treatment effects while maintaining a desirable aggregate-level summary by properly constraining the coupling. We construct transparent, individual confidence intervals for the estimated counterfactual outcomes, where the optimization objective controls the width of the confidence intervals.
\paragraph{Notation} We denote by $\|\cdot\|_2$ and $\langle \cdot, \cdot \rangle$ the Euclidean norm and inner product, respectively. For a square matrix $A$, let $\mathrm{tr}(A)$ denote its trace. For matrices $A, B$ of the same dimension, $\langle A, B \rangle:= \mathrm{tr}(A^\top B)$ is the Frobenius inner product. For any integer $n \in \mathbb{N}$, let $1_n = (1, \ldots, 1) \in \mathbb{R}^n$ denote the vector whose entries are all $1$. Let $\Delta_n := \{a \in \mathbb{R}_+^n : \sum_{i = 1}^{n} a_i = 1\}$ denote the simplex of probability vectors and let $\Delta_n^+ := \{a \in \Delta_n : a_i > 0 ~~ \forall i\}$ denote the interior of $\Delta_n$.
We first introduce Kernel Synthetic Coupling (KSC), a convexified matching method for missing value imputation. Later in Section (ref), we build individual confidence intervals around the imputed values. To streamline the exposition, we cast the KSC method following the notations in the potential outcome framework neyman1923,rubin1974estimating for a binary treatment---a leading example for missing value imputation.
We consider $N$ units indexed by $1, \ldots, N$, where each unit is associated with two potential outcome variables $Y_i(1), Y_i(0) \in \mathbb{R}$ under treatment and control, respectively, and a covariate vector $x_i \in \mathbb{R}^d$. We can observe only one of $Y_i(1), Y_i(0)$ depending on the treatment assignment $Z_i = 1$ or $Z_i = 0$, denoting that unit $i$ is treated or not. Accordingly, the observed outcome $Y_i$ of unit $i$ is given as $Y_i = Z_i Y_i(1) + (1 - Z_i) Y_i(0)$. We let $\mathcal{T} = \{i \in [N] : Z_i= 1\}$ and $\mathcal{C} = \{i \in [N] : Z_i= 0\}$ denote the sets of indices of the treated units and control units, respectively, so that $\mathcal{T} \cup \mathcal{C} = \{1, \ldots, N\}$. Also, let $N_t:= |\mathcal{T}|$ and $N_c:= |\mathcal{C}|$ denote the numbers of the treated and control units. The individual treatment effect of unit $i$ is defined as $\tau_i = Y_i(1) - Y_i(0)$, which is unobservable as only one of $Y_i(1), Y_i(0)$ can be observed.
Without loss of generality, we focus on imputing the missing potential outcomes of the treated units; counterfactual outcomes of the control units can be similarly obtained. We estimate $Y_j(0)$ by some $\widehat{Y}_j(0)$ for each treated unit $j$. This, in turn, allows us to estimate the individual treatment effect $\tau_j$ by $Y_j - \widehat{Y}_j(0) = Y_j(1) - \widehat{Y}_j(0)$. The proposed convexified matching combines two ideas. First, we want to find a match between the data cloud in the treated group and that in the control using the covariate information. Second, as in the synthetic control method, the matching quality is determined by how well it approximates the covariate of each treated unit by a convex combination of the covariates of control units. As a result, this procedure involves finding a matrix, which we call a coupling, indexed by a pair $(i, j)$ for $i \in \mathcal{C}$ and $j \in \mathcal{T}$ denoting the weight assigned to the covariate of control unit $i$ to approximate the covariate of treated unit $j$. We solve the following optimization to obtain the optimal coupling:
The first term in the objective function (ref) is the average of the squared approximation errors for covariate balancing. In the first constraint (ref), for each treated unit $j$, the synthetic weights $(\tfrac{\pi_{i j}}{1 / N_t})_{i \in \mathcal{C}}$ sum to $1$ and thus the term $\sum_{i \in \mathcal{C}} \tfrac{\pi_{i j}}{1 / N_t} x_i$ in (ref) is a convex combination to approximate $x_j$.
The second term in (ref) is an entropic regularization with regularization parameter $\lambda > 0$. Increasing the parameter $\lambda$ encourages the solution $\pi$ to be uniform, namely, $\pi_{i j} = \tfrac{1}{N_c \cdot N_t}$ for all $(i, j)$'s, which essentially corresponds to increasing the number of neighbors in the nearest neighbor matching.
The second constraint (ref) enforces that all control units contribute equally and thus seeks a matching between treated and control data sets; this constraint can be generalized to incorporate inverse propensity scores, to be shown in the next section. If this constraint is dropped, the optimization can be decoupled into $N_t$ independent optimization of the same type, each with one treated unit $j \in \mathcal{T}$, also solvable by our optimization algorithm detailed in Section (ref). We shall show in a second that this constraint will be crucial. It enables the convex program to synthesize granular information while enforcing, at the coarse level, that the answers agree with typical estimators for the ATE or ATT.
Finally, the optimization variable $\pi \in \mathbb{R}_+^{N_c \times N_t}$ is a matrix whose entries are indexed by $(i, j)$ for $i \in \mathcal{C}$ and $j \in \mathcal{T}$, instead of the usual indexing by $i \in \{1, \ldots, N_c\}$ and $j \in \{1, \ldots, N_t\}$, to denote that $\pi_{i j}$ is the weight between control unit $i$ and treated unit $j$. We call $\pi \in \mathbb{R}_+^{N_c \times N_t}$ satisfying the constraints (ref) and (ref) a coupling, which generalizes doubly stochastic matrices to non-square matrices with prescribed row and column sums.
Upon obtaining the solution, denoted by $\widehat{\pi}$, we impute the counterfactual outcome $Y_j(0)$ of treated unit $j$ by the convex combination of the outcomes of control units with the weights $(N_t \widehat{\pi}_{i j})_{i \in \mathcal{C}}$:
which in turn leads to the individual treatment effect estimate $\widehat{\tau}_j := Y_j - \widehat{Y}_j(0)$. By the constraint (ref), the average of the estimated individual treatment effects coincides with the difference in the mean outcomes between the treatment and control groups denoted as $\widehat{\tau}^{\mathsf{DiM}}$:
Therefore, the proposed convexified matching allows for estimating individual treatment effects while maintaining the desired ATT estimate $\widehat{\tau}^{\mathsf{DiM}}$.
The resulting optimization is a smooth convex optimization problem. To see this, notice that the first term of the objective function (ref) is the following quadratic function of $\pi$:
where $K_{c c} \in \mathbb{R}^{N_c \times N_c}$, $K_{c t} \in \mathbb{R}^{N_c \times N_t}$, and $K_{t t} \in \mathbb{R}^{N_t \times N_t}$ are the Gram matrices whose entries are inner products of the covariates, namely, $K_{c c} = (\langle x_i, x_{i'} \rangle)_{i, i' \in \mathcal{C}}$, $K_{c t} = (\langle x_i, x_j \rangle)_{(i, j) \in \mathcal{C} \times \mathcal{T}}$, and $K_{c, t} = (\langle x_j, x_{j'} \rangle)_{j, j' \in \mathcal{T}}$. This quadratic function is convex as the Gram matrix $K_{c c}$ is positive semidefinite. The second term---the entropy term---is convex. The constraint set resulting from (ref) and (ref) is a convex polytope consisting of couplings. In Section (ref), we introduce and analyze efficient algorithms to solve this optimization problem based on a simple iterative matrix scaling procedure called the Sinkhorn algorithm sinkhorn1967diagonal.
Lastly, we discuss several existing methods and concepts related to the proposed method. abadie2021penalized considers multiple treated units ($N_t > 1$) and proposes running synthetic controls separately for treated units with a different regularization function. The proposed method differs from running synthetic controls separately because of the constraint (ref). If there is only one treated unit ($N_t = 1$) and no constraint (ref), the proposed method is equivalent to the standard synthetic control with entropic regularization as in hainmueller2012entropy which proposes entropic regularization for covariate balancing. Also, notice that the first term of (ref) is upper bounded under the constraint (ref) as follows:
where the inequality uses Jensen's inequality. This upper bound analytically shows that the average approximation error based on convex combinations is smaller than that based on pairwise distances, which aligns with the motivation mentioned in the previous section. Replacing the first term of the objective function (ref) by the above upper bound---while maintaining the constraints (ref), (ref)---leads to the optimal transport problem villani_2003 with entropic regularization cuturi2013sinkhorn.
We present two extensions of the simple formulation introduced in Section (ref). The first is extending the main insights to cover the nonlinear case when the potential outcomes can be nonlinear functions of covariates, using a kernel trick shawe2004kernel,steinwart_christmann_2008. The second is incorporating notions of propensity scores rosenbaum1983central into our synthetic coupling formulation.
\paragraph{Kernel Synthetic Coupling}
We can kernelize the presented simple formulation to extend to nonlinear cases. The kernel trick is widely used in machine learning to modify methods that originated in the linear setting to accommodate nonlinearity, such as ridge regression, support vector machines, and principal component analysis. Recall that our objective (ref) relies on the Gram matrices as shown in (ref). Therefore, to kernelize the simple formulation, we can replace the Gram matrices with the kernel matrices, namely, letting $K_{c c} = (k(x_i, x_{i'}))_{i, i' \in \mathcal{C}}$, $K_{c t} = (k(x_i, x_j))_{(i, j) \in \mathcal{C} \times \mathcal{T}}$, and $K_{t t} = (k(x_j, x_{j'}))_{j, j' \in \mathcal{T}}$ for some kernel function $k \colon \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$. The kernelized formulation is essentially the same convex optimization problem. It can be tackled by the algorithms we introduce in Section (ref).
Applying the above kernel trick is equivalent to replacing the first term of (ref) with the approximation error in the reproducing kernel Hilbert space (RKHS) $\mathcal{H}$ associated with the kernel $k$, namely, $\frac{1}{2 N_t} \sum_{j \in \mathcal{T}} \|\phi_{x_j} - \sum_{i \in \mathcal{C}} \tfrac{\pi_{i j}}{1 / N_t} \phi_{x_i}\|_\mathcal{H}^2$, where $\|\cdot\|_\mathcal{H}$ is the norm in the RKHS $\mathcal{H}$ and $\phi_{x} = k(x, \cdot)$ denotes the canonical feature of $x \in \mathbb{R}^d$. In other words, we generalize the matching of covariates in the Euclidean space to the matching of embedded canonical features in the RKHS $\mathcal{H}$, which leads to the following kernelized formulation:
where $v = (v_j)_{j \in \mathcal{T}} \in \Delta_{N_t}^+$ and $w = (w_i)_{i \in \mathcal{C}} \in \Delta_{N_c}^+$ are strictly positive probability vectors that generalize the constraints (ref) and (ref) to allow treated and control units to have different weights.\footnote{Whenever $v = \frac{1}{N_t} 1_{N_t}$, we will drop the dependency on $v$ and abbreviate as $\mathrm{KSC}^{w}_\lambda$.} With the optimal coupling $\widehat{\pi}$, we impute using convex combination same as in (ref), a classic idea in nonparametric regression.
\paragraph{Propensity Score} Using arbitrary weights through $v = (v_j)_{j \in \mathcal{T}} \in \Delta_{N_t}^+$ and $w = (w_i)_{i \in \mathcal{C}} \in \Delta_{N_c}^+$ allows for various average treatment effect estimators, particularly including those incorporating propensity scores. To estimate the ATT, we can set $v = \frac{1}{N_t} 1_{N_t}$ as before, and the average of the estimated individual treatment effects of the treated leads to
Now, let $\widehat{p} \colon \mathbb{R}^d \to (0, 1)$ be a suitable estimator---often estimated by logistic regression---of the propensity score, $x \mapsto \mathbb{P}(W = 1 \, | \, X = x)$. If the weights $w$ are chosen based on the propensity scores as
the aggregate effects match the normalized Inverse Probability Weighting (IPW) estimator hirano2001estimation, that is,
Namely, the proposed KSC method is guaranteed to yield individual treatment effect estimates that conform to (ref) at the aggregate level.
Similarly, solving $\mathrm{KSC}^{v,w}_\lambda$ in (ref) by setting
we obtain an coupling $\widehat{\pi}$ that matches the IPW estimator for ATE
Now, we demonstrate our imputation method, KSC, using a real-world data set. This section aims to concisely present the features of KSC and compare it with other imputation methods; further details and robust evaluations are provided in Section (ref) together with the inference procedure, to be presented in the next section.
We apply the proposed method to evaluate the National Supported Work (NSW) demonstration program, first analyzed by lalonde1986evaluating. This data set has been widely studied in the program evaluation literature. The NSW program was conducted from 1975 to 1978 in the United States, which aimed at providing a job training program to disadvantaged workers, where the treatment, namely, the training, was randomly assigned to some of them. The data we consider here is a specific subset of the experimental data used in dehejia1999causal, which consists of $N = 445$ subjects with $N_t = 185$ treated units and $N_c = 260$ control units. The outcome of interest $Y_i$ is the post-treatment earnings recorded in 1978. Each individual is associated with six variables: age, years of education, and four indicator variables denoting whether the individual is black, Hispanic, married, and a high school dropout. Following dehejia1999causal, the pretreatment outcomes, the earnings in 1974 and 1975, and two indicator variables denoting whether these pretreatment earnings are zero, are included as covariates. Accordingly, we have $x_i \in \mathbb{R}^d$ with $d = 10$.
Figure (ref)(a) and Figure (ref)(b) show the scatter plots of the imputed counterfactual outcomes of the treated units, where the $x$-coordinate is $\widehat{Y}_j(0)$ imputed by the proposed method, while the $y$-coordinate is the imputed counterfactual outcomes by the nearest neighbor matching with replacement or the regression imputation method. The proposed method relies on the linear kernel following the simple formulation in Section (ref) with $\lambda \in \{0.001, 0.01\}$, the nearest neighbor matching uses $5$ and $10$ nearest neighbors, and the regression method fits a linear regression model to the observed outcome using all the covariates and the treatment indicator. At first glance, the scatter plots visually suggest that all three imputation methods provide comparable imputed counterfactual outcomes. One crucial difference is that the average of the imputed counterfactual outcomes and, thus, the estimated individual treatment effects vary across the methods. The ATT estimate computed by the proposed method always matches---regardless of the choice of $\lambda$---the mean difference between the treatment and control groups, which is roughly 1794.3, as shown in (ref)---the unbiased estimator of the ATT in the experimental setup. In contrast, the $k$-nearest neighbor matching with $k = 5$, with $k = 10$, and the regression imputation method produce the ATT estimates 2030.5, 1776.6, and 1706.2, respectively. While the average of the imputed counterfactual outcomes does not depend on the choice of $\lambda$, we can see that it determines the level of individualization. The histograms along the $x$-axix show that the imputed counterfactual outcomes based on the proposed method with $\lambda = 0.001$ are more dispersed than those based on $\lambda = 0.01$; the latter is based on larger regularization leading to more uniform weights, which is similar to the fact that $k = 5$ is more dispersed---along the $y$-axis---than $k = 10$ in the nearest neighbor matching.
Since lalonde1986evaluating, there has been a debate on the credibility of the average treatment effect estimation using observational control groups. Since such control groups can be significantly different from the control group of the experimental data, the resulting ATE estimates can deviate much from the experimental benchmark, say, the mean difference of 1794.3 computed earlier, which has been the central argument of lalonde1986evaluating. Later, several methods dehejia1999causal,dehejia2002propensity have been proposed to recover the experimental benchmark estimate using observational control groups. Here, we focus on one such method based on the normalized inverse probability weighting (IPW) estimator (ref) and apply the KSC method with the weights based on the propensity scores. To this end, we consider a control group from the Panel Study of Income Dynamics (PSID) data, consisting of 2490 units, with which we estimate the propensity scores by logistic regression. The resulting ATT estimate by $\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$ is 2579.7, which is far from the benchmark of 1794.3, rooted in the stark difference between the experimental and observational control groups. To remedy this, we trim the PSID data by only taking the units whose propensity scores are in $[0.05, 0.95]$, yielding 214 units. The resulting ATT estimate by $\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$ is then 1748.0, which resembles the experimental benchmark. Taking this trimmed PSID data as a control group, we apply the proposed method with the weights $w_i$ based on the propensity scores as in (ref) and compare with the matching and regression methods. The results are shown in Figure (ref)(c) and Figure (ref)(d). Again, though visually similar in distributions, the ATT estimates differ across the methods, where the proposed method with the propensity score-based weights provides the ATT estimate 1748.0, regardless of $\lambda$. In contrast, the $k$-nearest neighbor matching with $k = 5$, with $k = 10$, and the regression method produce the ATT estimates 1567.3, 1003.0, and 1065.5, respectively. This example illustrates that the proposed imputation method allows for producing granular information, namely, the individual treatment effects, from a coarse level estimate by $\widehat{\tau}^{\mathsf{IPW}}_{\mathsf{ATT}}$. The coarse level estimate is closer to the experimental benchmark estimate than inherently imputation-based methods such as the nearest neighbor matching or regression.
Constructing confidence intervals for the estimated individual treatment effects is conducive to reliable decision-making. This section introduces and analyzes the method to construct individual confidence intervals based on the proposed convexified matching. To this end, we build a confidence interval for the missing counterfactual outcome $Y_j(0)$ for each treated unit $j \in \mathcal{T}$. We first introduce appropriate model assumptions. Then, we elucidate the construction of individual confidence intervals, followed by a numerical example, and discuss the coverage and efficiency of the proposed method.
\paragraph{The Model} Let $\mathcal{H}$ be the reproducing kernel Hilbert space (RKHS) associated with a kernel function $k \colon \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$, where $\langle \cdot, \cdot \rangle_\mathcal{H}$ denotes the RKHS inner product, $\|\cdot\|_\mathcal{H}$ denotes the RKHS norm, and $\phi_{x}:= k(x, \cdot) \in \mathcal{H}$ is the corresponding reproducing kernel function. We postulate the following model for the potential outcomes: for each $i = 1, \ldots, N$, letting $x_i \in \mathbb{R}^d$ be the covariate of unit $i$, we have
where $f_0, f_1 \in \mathcal{H}$ are some functions in the RKHS and $\sigma_0, \sigma_1 \ge 0$. Let $\mathbf{X}:= [x_1, \ldots, x_N]$ denote the fixed design matrix and let $\mathbf{Z} := [Z_1, \ldots, Z_N]$ denote the possibly random treatment assignment vector.
Under this setup, for each treated unit $j \in \mathcal{T}$, the goal is to build an individual confidence interval for its missing counterfactual outcome $Y_j(0)$. More specifically, we aim to construct a confidence interval $[\widehat{L}_j, \widehat{U}_j]$ that contains the conditional expectation $f_0(x_j) = \operatorname*{\mathbb{E}}[Y_j(0) \, | \, x_j] = \operatorname*{\mathbb{E}}[Y_j(0) \, | \, \mathbf{X}, \mathbf{Z}]$ with a desired level of confidence as follows:
namely, the probability that the interval $[\widehat{L}_j, \widehat{U}_j]$ contains the conditional expectation $f_0(x_j)$ is at least $1 - \alpha$, where the probability is conditional on the design matrix $\mathbf{X}$ and the treatment assignments $\mathbf{Z}$.
\paragraph{Construction} We consider the following Kernerlized Synthetic Coupling (KSC) discussed in Section (ref). To construct an individual interval for the counterfactual control outcome of each treated unit, we solve (ref) with uniform weights $v = \frac{1}{N_t}$ on the treated, which we rewrite in the matrix form as follows: let $(w_i)_{i \in \mathcal{C}}$ be positive numbers such that $\sum_{i \in \mathcal{C}} w_i = 1$,
Let $\widehat{\pi}$ be the solution of the above optimization problem and define a probability transition matrix $\widehat{P} = (\widehat{p}_{i j})_{(i, j) \in \mathcal{C} \times \mathcal{T}}$, where $\widehat{p}_{i j} = N_t \widehat{\pi}_{i j}$ so that $\sum_{i \in \mathcal{C}} \widehat{p}_{i j} = 1$ for any $j \in \mathcal{T}$. For each treated unit $j$, recall the imputed value $\widehat{Y}_j(0)$ defined in (ref), for which we construct a confidence interval $[\widehat{L}_j, \widehat{U}_j]$, where $\widehat{L}_j$ and $\widehat{U}_j$ are defined as follows:
Here, $z_{1 - \frac{\alpha}{2}}$ is the $(1 - \frac{\alpha}{2})$-quantile of the standard normal distribution and $\widehat{\theta}, \widehat{\sigma}_0$ are suitable estimates of the norm $\|f_0\|_\mathcal{H}$ and the standard deviation $\sigma_0$, respectively.
We apply kernel ridge regression to the control group data to estimate $\|f_0\|_\mathcal{H}$ and $\sigma_0$. We can estimate $f_0 \in \mathcal{H}$ by $\widehat{f}_0(x) = \sum_{i \in \mathcal{C}} \widehat{\beta}_i k(x_i, x)$, where $\widehat{\beta} = (K_{c c} + \rho I_{N_c})^{-1} Y_c \in \mathbb{R}^{N_c}$ with a kernel matrix $K_{c c}$ defined in Section (ref), a vector of control outcomes $Y_c = (Y_i)_{i \in \mathcal{C}} \in \mathbb{R}^{N_c}$, and a ridge regularization parameter $\rho > 0$. Then, we estimate $\|f_0\|_\mathcal{H}$ by
and the variance $\sigma_0^2$ by
Before diving into each component of the construction, we first look at a simple numerical example to illustrate the performance of the proposed method.
\paragraph{Numerical Example} We consider a set of fixed one-dimensional covariates $x_1, \ldots, x_N \in [0, 1]$ with $N = 500$ and we randomly choose $N_t = 200$ treatment units. Let $\mathcal{H}$ be the RKHS associated with the Gaussian kernel $k(x, x') = \exp(-\gamma \cdot |x - x'|^2)$, where $\gamma = 2.5$, and suppose $f_0(x) = k(0.5, x)$, which yields $\|f_0\|_\mathcal{H} = (k(0.5, 0.5))^{1 / 2} = 1$. Then, we generate the potential outcomes by $Y_i(0) = f_0(x_i) + \varepsilon_{0, i}$, where $\varepsilon_{0, 1}, \ldots, \varepsilon_{0, N}$ are independently drawn from $N(0, \sigma_0^2)$. As $x_1, \ldots, x_N$ and $Z_1, \ldots, Z_N$ are fixed, the randomness is solely from the residuals $\varepsilon_i$'s, which we sample 1000 times.
Figure (ref) shows the constructed confidence intervals $[\widehat{L}_j, \widehat{U}_j]$ for each treated unit $j$ for different values of $\sigma_0$. For each $\sigma_0 \in \{0.1, 1, 3\}$, we solve (ref) for $\lambda \in \{0.1, 0.01, 0.001\}$ with uniform weights, run the kernel ridge regression to obtain $\widehat{\theta}, \widehat{\sigma}_0$, and produce $[\widehat{L}_j, \widehat{U}_j]$ which is shown as blue shaded regions in the plots. We compare this interval with an “oracle” interval $[L_j^\ast, U_j^\ast]$, which---as we will explain in the next section--- is guaranteed to have an exact $1 - \alpha$ coverage, namely,
By construction, if the estimates $\widehat{\theta}, \widehat{\sigma}_0$ are accurate enough, the interval $[\widehat{L}_j, \widehat{U}_j]$ is supposed to contain the oracle interval $[L_j^\ast, U_j^\ast]$, which can be verified in Figure (ref). For each $\sigma_0$, decreasing $\lambda$ results in smaller approximation errors, namely, the bias correction reflected in $\left(K_{t t} + \widehat{P}^\top K_{c c} \widehat{P} - 2 K_{c t}^\top \widehat{P} \right)_{j j}^{1/2}$, and the estimated counterfactual outcomes (solid blue curves) that are more wiggly as the obtained weights have more variations. As we will see, $[\widehat{L}_j, \widehat{U}_j]$ is wider than $[L_j^\ast, U_j^\ast]$ by the twice of the bias correction, namely,
We can verify from Figure (ref) that this gap between the two intervals gets smaller as $\lambda$ decreases, consistent with the above observation that the bias decreases.
Figure (ref) plots the coverage of $[\widehat{L}_j, \widehat{U}_j]$ and $[L_j^\ast, U_j^\ast]$, namely, the probability that they contain $f_0(x_j)$ estimated using the 1,000 samples. When $\lambda = 0.1$, we can see that the coverage of $[\widehat{L}_j, \widehat{U}_j]$ is almost $1$, meaning that the interval is too conservative, which results from the fact that the bias term is not small enough compared to the variance of the residuals. As $\lambda$ decreases, we can see that the coverage of $[\widehat{L}_j, \widehat{U}_j]$ gets closer to $1 - \alpha = 0.95$, which is consistent with the fact that the bias decreases. Particularly, when $\sigma_0 = 3$ and $\lambda = 0.001$, the coverage of $[\widehat{L}_j, \widehat{U}_j]$ is closest to the coverage of the oracle interval $[L_j^\ast, U_j^\ast]$, suggesting that the bias is dominated by the variance. Meanwhile, for $(\sigma_0, \lambda) = (1, 0.01)$ and $(\sigma_0, \lambda) = (3, 0.1)$, the coverage of $[\widehat{L}_j, \widehat{U}_j]$ is larger than the coverage of the oracle interval $[L_j^\ast, U_j^\ast]$ by approximately a constant, which suggests that the bias term and the variance term are nearly balanced. In the next section, we provide theoretical guidance to choose the regularization parameter $\lambda$ to balance the bias and variance tradeoff.
Now, we dive into each component in constructing individual confidence intervals to dissect the coverage and efficiency. We discuss how the entropic regularization parameter $\lambda$ drives the bias and variance tradeoff, as seen in the numerical simulations. Furthermore, we will show that optimizing the entropic-regularized convexified matching program (ref) with the desired $\lambda$ chosen above is closely connected to optimizing the average of the squared lengths of the individual confidence intervals $\frac{1}{N_t} \sum_{j \in \mathcal{T}} \mathrm{Len}_j^2$.
\paragraph{Point Estimate} For each treated unit $j \in \mathcal{T}$, the center of the confidence interval $[\widehat{L}_j, \widehat{U}_j]$ is the imputed counterfactual outcome $\widehat{Y}_j(0)$:
This weighted average estimator leverages all the information in $\mathbf{X}$ to anchor an optimal coupling (or matching) between the treated and control units for synthesizing the controls. It uses convex combinations of the control units to synthesize the granular counterfactual outcomes for the treated units, while balancing the aggregate level statistics, the ATT. We remark that the optimal coupling $\widehat{\pi}$---and thus $\widehat{P} = N_t \widehat{\pi}$ ---is calculated based solely on $\mathbf{X}$ and $\mathbf{Z}$, not the outcomes $Y_i$'s. Recall the following bias-variance decomposition of the point estimate $\widehat{Y}_j(0)$:
\paragraph{Bias-Awareness} The point estimate could be biased as
where the last equality uses the reproducing property of the RKHS. Since $f_0$ is unknown and estimating such a function uniformly on the support of the covariate vector can be difficult, we rely on a bias-awareness correction by invoking the Cauchy-Schwarz inequality:
This explains subtracting and adding the term $\widehat{\theta} \cdot \left(K_{t t} + \widehat{P}^\top K_{c c} \widehat{P} - 2 K_{c t}^\top \widehat{P} \right)_{j j}^{1/2}$ in (ref) and (ref), respectively, where $\widehat{\theta}$ estimates the unknown quantity $\|f_0\|_\mathcal{H}$.
\paragraph{Variance} Under the model,
If we knew the bias $\mathcal{B}_j$, we could construct the following oracle confidence interval $[L_j^\ast, U_j^\ast]$:
Due to (ref), the oracle interval $[L_j^\ast, U_j^\ast]$ is guaranteed to have the exact $1 - \alpha$ coverage, namely, (ref) holds. However, since $\mathcal{B}_j$ is unknown, we rely on the aforementioned bias-awareness correction to construct the confidence interval $[\widehat{L}_j, \widehat{U}_j]$, obtained by inflating the oracle interval using (ref).
\paragraph{Efficiency: Interval Length} We now explain how the width of the confidence interval is related to our convexified matching objective. As a result, the explanation will shed light on why the entropic regularization serves as a tuning parameter to trade off bias and variance for overall interval length efficiency. First, we need the following simple facts on the Kullback-Leibler divergence, Hellinger distance, and $\chi^2$ distance.
The above fact will help us relate the entropy term and the standard Euclidean norm of the weights. For a fixed $j \in \mathcal{T}$, let us consider that $(\widehat{p}_{i j})_{i \in \mathcal{C}} \in \Delta_{N_c}$, the $j$-th column of $\widehat{P}$, is delocalized in the sense $\max_{i \in \mathcal{C}} \widehat{p}_{i j} \le (M/2-1)/ N_c$ with some constant $M > 4$. The interval length is as follows:
From $\widehat{\pi}_{i j} = \widehat{p}_{i j} / N_t$ and Proposition (ref), we deduce the following bounds:
Averaging over $j \in \mathcal{T}$, and barring the constant $M$, the upper and lower bounds on the right-hand sides remind us of the entropic regularized convexified matching objective in (ref), with the regularization parameter
The above tradeoff confirms the empirical findings: for a high signal-to-noise problem, to balance bias and variance, we need a smaller $\lambda$ for better approximation error; for a low signal-to-noise problem, a larger $\lambda$ is preferred to reduce variance. Such a phenomenon confirms the observations in Figure (ref). Searching for a coupling minimizing the objective of convexified matching is closely related to minimizing the interval length, provided the regularization parameter $\lambda$ is appropriately chosen. Therefore, optimizing the entropic regularized convexified matching program (ref) with the desired $\lambda$ chosen above is closely connected to optimizing the average of the squared lengths of the individual confidence intervals $\frac{1}{N_t} \sum_{j \in \mathcal{T}} \mathrm{Len}_j^2$.
This section introduces and analyzes optimization algorithms to solve the convexified matching problem, KSC defined in (ref) with uniform weights $v = \frac{1}{N_t}$ on the treated, as shown in (ref). Notice that it is an instance of the constrained nonlinear optimization problem: for $n, m \in \mathbb{N}$, $a \in \Delta_n^+$, $b \in \Delta_m^+$, $\lambda > 0$, and $g \colon \mathbb{R}^{n \times m} \to \mathbb{R}$,
where $\Pi_{a, b} = \{\pi \in \mathbb{R}_+^{n \times m} : \pi 1_m = a ~ \text{and} ~ \pi^\top 1_n = b\}$ and $h \colon \mathbb{R}_+^{n \times m} \to \mathbb{R}$ is the entropy function defined by $h(\pi) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \pi_{i j} \log \frac{\pi_{i j}}{e}$. The convexified matching amounts to the case where $g$ is a convex quadratic function, in which (ref) is a convex program taking the form (ref). It can be solved using the interior point method in $\mathrm{poly}(nm)$ time, possibly slow in practice for coupling matrices with $nm$ large.
In this section, we devise a faster algorithm that requires an almost dimension-free number of oracle calls to a subroutine---a simple iterative matrix scaling procedure called the Sinkhorn algorithm sinkhorn1967diagonal. When $\lambda$ is sufficiently large, the algorithm compares favorably to the interior point method, requiring only a $\log(1/\varepsilon)$ number of Sinkhorn subroutines. Later, we modify and extend the algorithm to cover all $\lambda \in \mathbb{R}_+$ by connecting to the steepest descent method under the Kullback-Leibler divergence. The convergence is admittedly slower for small $\lambda$: $\log(n m)/\varepsilon$ number of Sinkhorn subroutines is required, but it still compares favorably to the interior point method.
\paragraph{Additional notation} For a matrix $A \in \mathbb{R}^{n \times m}$, let $\|A\|_{\mathrm{F}}$ denote its Frobenius norm, $\|A\|_1 := \sum_{i = 1}^{n} \sum_{j = 1}^{m} |A_{i j}|$ denote its $L^1$ norm, and $\|A\|_\infty := \max_{1 \le i \le n} \max_{1 \le j \le m} |A_{i j}|$ denote its $L^\infty$ norm. Let $h$ keep its definition as the entropy function, and let $\nabla h(A)$ denote the gradient of $h$, which is $\log(A)$, the entrywise logarithm of $A$, provided all entries of $A$ are positive. For any $a \in \Delta_n$ and $b \in \Delta_m$, let $\Pi_{a, b}^{+} := \{A \in \Pi_{a, b} : A_{i j} > 0 ~~ \forall i, j\}$. Let $\Delta_{n, m} = \{A \in \mathbb{R}_+^{n \times m} : \sum_{i = 1}^{n} \sum_{j = 1}^{m} A_{i j} = 1\}$. For vectors $u, v \in \mathbb{R}^n$, let $\frac{u}{v}$ denote the entrywise division provided all the entries of $v$ are nonzero.
When $g$ is convex and $\lambda > 0$, the strong convexity of $h$ on $\Pi_{a, b}$ implies that (ref) admits a unique minimizer. Moreover, as $h$ prevents the minimizer from having a zero entry, (ref) admits a unique minimizer on $\Pi_{a, b}^{+}$, which is the interior of $\Pi_{a, b}$. It turns out that this unique minimizer of (ref) is the fixed point of an operator related to the entropic regularized optimal transport problem cuturi2013sinkhorn, as we shall show in Proposition (ref). Later, in Theorem (ref), we derive convergence to the minimizer via the iterative matrix scaling subroutine. This subroutine is referred to as the Sinkhorn algorithm sinkhorn1967diagonal, which we introduce below.
For simplicity, Algorithm (ref) states the Sinkhorn algorithm such that the column sum of $P^{(k)}$ matches with $b$, namely, $(P^{(k)})^\top 1_n = b$ for any $k \ge 0$. Then, to assess the accuracy of the scaled matrix $P^{(k)}$, we may only need to check a suitable discrepancy between $a$ and the row sum of $P^{(k)}$, that is, $P^{(k)} 1_m$.
The following proposition shows that when $g$ is convex, the unique minimizer of (ref) is the fixed point of an operator given by the composition of $\Phi_\lambda$ and $\nabla g$.
Based on the established connection between the optimality condition of (ref) and the Sinkhorn algorithm, we propose two algorithms to solve (ref) in the following sections.
By Proposition (ref), finding the fixed point of $T_\lambda$ is equivalent to solving (ref). Therefore, we propose solving (ref) by the fixed-point iterations based on the operator $T_\lambda \colon \mathbb{R}^{n \times m} \to \Pi_{a, b}^{+}$, which is simply choosing the initial point $\pi^{(0)} \in \Pi_{a, b}$ properly and iterating $\pi^{(k + 1)} \leftarrow T_\lambda(\pi^{(k)})$ for $k \ge 0$, where $T_\lambda(\pi^{(k)}) = \Phi_\lambda(\nabla g(\pi^{(k)}))$ is approximated by $\mathrm{Sinkhorn}(e^{-\nabla g(\pi^{(k)}) / \lambda}, a, b)$ as explained in Definition (ref). Algorithm (ref) summarizes this procedure.
We show that $T_\lambda$ is a contraction if the gradient of $g$ is Lipschitz and $\lambda$ is larger than the Lipschitz constant. In such a case, Algorithm (ref) converges to the fixed point quickly. When $g$ is a quadratic function, we can write the Lipschitz constant in terms of the $L^\infty$ norm of the Hessian matrix, which is independent of the input dimension $n m$. We do not require $g$ to be convex for this result; $T_\lambda$ admits a unique fixed point due to the contraction argument, which is independent of Proposition (ref). The role of Proposition (ref) is to translate this convergence result into the convergence to the minimizer of (ref) for convex $g$ by equating the fixed point to the minimizer.
The fixed-point algorithm may not converge for sufficiently small $\lambda$. To fill in the gap when $\lambda \in [0, \| H\|_{\infty})$ left unanswered by Theorem (ref), we employ the steepest descent method under Bregman divergence to solve (ref) as follows: let $f = g + \lambda h$,
where $D_h(\pi, \pi^{(k)}) = h(\pi) - h(\pi^{(k)}) - \langle \nabla h(\pi^{(k)}), \pi - \pi^{(k)} \rangle$ is the Bregman divergence, which equals the Kullback-Leibler divergence on $\Pi_{a, b}$, and $\tau_k > 0$ is a suitable step size. With the entropic regularization $h$, we implement the steepest descent over the polytope $\Pi_{a, b}$ under the Kullback-Leibler divergence. We emphasize that this can be solved using the Sinkhorn algorithm as well since (ref) is equivalent to
The above can be approximated by applying the Sinkhorn algorithm to $\nabla g(\pi^{(k)}) + (\lambda - \tau_k^{-1}) \nabla h(\pi^{(k)})$. Algorithm (ref) summarizes this procedure. Note that letting $\tau_k = \lambda^{-1}$ for all $k \ge 0$ in Algorithm (ref) recovers Algorithm (ref).
We show that Algorithm (ref) outputs a sequence that converges to the minimum of (ref) provided the step size is below a certain threshold.
In this section, we numerically investigate our imputation and individualized inference method using the NSW data and compare it with other methods, continued from Section (ref).
We conduct a robustness check of the imputation results in Section (ref) on the NSW data by specifying different nonlinear kernels. We consider the RBF kernel and the polynomial kernel of degree two. The results are shown in Figure (ref) and Figure (ref). Though the overall patterns are similar to those using the linear kernel, we can observe some differences. First, for the experimental data, the imputed counterfactual outcomes via the RBF kernel with $\lambda = 0.001$ show the largest variability. In Figure (ref)(a), there are two treated units whose imputed counterfactual outcomes are above 15000, which do not occur in the results via the linear or polynomial kernel; as a result, the right tail of the histogram along the $x$-axis is longer than the other cases. For the PSID data, we can see that using the polynomial kernel leads to the most spread out imputed outcomes, visualized by the right tails of the histograms along the $x$-axis in (a) and (b) of Figure (ref) which are longer than the other cases. We point out that unlike the NN-matching method and our KSC method, whose counterfactual distributions are supported on $\mathbb{R}_+$ even for the nonexperimental PSID data, the regression imputation method presents a wider support for imputed values, even with negative outcomes as seen across Figures (ref), (ref) and (ref).
We apply the method presented in Section (ref) to construct confidence intervals around the imputed counterfactual outcomes. Figure (ref) (a) shows the results for the experimental data using the linear kernel with $\lambda \in \{0.001, 1\}$, which plots the estimated individual treatment effects (ITEs), $\widehat{\tau}_j = Y_j - \widehat{Y}_j(0)$, against the logarithm of earnings in 1975, where the confidence intervals are shown as error bars around the estimated ITEs. Here, only 74 treated units (out of 185 treated units) with positive earnings in 1975 are shown. We pick this index merely for visualization of individualization. For comparison, we also show the blue points representing the ITEs based on imputing all the treated units with the mean of the control group outcomes, along with the confidence intervals based on the standard error of the mean; the case where there is no individualization in the imputed control outcomes, and thus no meaningful individual confidence intervals addressing the bias. Let us first focus on the point estimates of the ITEs. We can see that most of the estimated ITEs fall in the interval $(-10000, 10000)$, seemingly uncorrelated with the earnings in 1975, while there are seven individuals whose ITEs are above 10000. These seven units with the largest ITEs are likely to be the ones who benefit most from the treatment, standing out from the rest of the treatment group. The confidence intervals help us decide how sure we are about these seemingly significant ITEs. We can confidently say that the top two ITEs are extremely large compared to the rest of the ITEs since the lower ends of their individual confidence intervals are larger than the upper ends of the rest of the ITEs, regardless of the choice of $\lambda$. Now, looking at the confidence intervals, rather than only looking at the point estimates, the four units whose ITEs are between 10000 and 30000 are less distinguishable from the rest. The error bars do help to suggest which units benefit the more here. For the third largest ITE, around 20000, the conclusion may vary depending on $\lambda$. If we read the confidence intervals corresponding to $\lambda = 1$ (red, dashed), this unit seems to be more significant than the rest, while this conclusion is less assertive for $\lambda = 0.001$ as the confidence interval of this unit overlaps with the interval $(-15000, 20000)$ which contains the rest of the confidence intervals. In other words, using $\lambda = 0.001$ leads to the most conservative conclusion regarding who benefits the most from the treatment as the confidence intervals are wider. Unlike the confidence intervals obtained by the proposed method, the intervals based on the mean imputation (shown in blue) are extremely narrow, which does not take into account individual heterogeneity, overlooking a potentially large bias; the confidence interval here is useless for individual decision making and may result in overly optimistic conclusions, for instance, essentially the majority benefits significantly from the program. (b) and (c) of Figure (ref) repeat the same analysis using the RBF kernel and the polynomial kernel of degree two, respectively. The overall patterns are similar to those using the linear kernel, but the confidence intervals tend to be wider. Accordingly, even the second largest ITE above 30000 is not as distinguishable from the rest as in the linear kernel case.
The results based on the PSID data are shown in Figure (ref). One notable difference is that the intervals are generally much wider than those of the NSW experimental data. One possible explanation is that the PSID data, even after trimmed, differs significantly from the experimental data, leading to the matching with more conservative error estimates. Accordingly, the confidence intervals are wider, which makes it harder to make a decisive conclusion with confidence.
Lastly, we briefly comment on the choice of $\rho$, the regularization parameter of the kernel ridge regression mentioned in Section (ref). Increasing $\rho$ leads to a smaller estimate of $\|f_0\|_\mathcal{H}$, leading to a smaller bias correction. Increasing $\rho$ will typically result in a larger estimate of $\sigma_0$, implying a large variability. The correct choice of $\rho$ depends on the underlying signal-to-noise ratio for the data-generating process in the control outcomes, namely, the complexity of the function vs.\ the amount of noise. In all examples of this section, we choose $\rho$ using 5-fold cross-validation.
We proposed a convexified matching method for missing value imputation and individualized inference, integrating favorable features from optimal matching, regression imputation, and synthetic control. We impute counterfactual outcomes based on convex combinations of observed outcomes, defined by an optimal coupling between the treated and control data sets. Finding an optimal coupling is a convex relaxation to the combinatorial optimal matching problem, for which we propose efficient algorithms based on matrix scaling. Unlike existing imputation methods, we begin with a desirable aggregate-level summary and estimate granular-level individual treatment effects by properly constraining the coupling so that the estimated individual effects are consistent with the aggregate summary. We provided a method to construct individual confidence intervals for the estimated counterfactual outcomes, along with a simulation study to demonstrate the effectiveness of our method. The simulation confirms our theoretical insight on the level of entropic regularization needed. We establish that entropic regularization plays a crucial role in both efficiency in inference and computation: first, the convexified matching objective is gauged for minimizing the width of the individual confidence intervals, trading off bias and variance; second, the entropic regularization enables us to design fast algorithms. We demonstrated the empirical performance of our method on the NSW data, using both experimental and nonexperimental data sets, with a comparison to existing methods and robustness checks.
\printbibliography