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.
49,546 characters · 12 sections · 94 citation commands
Synthetic Principal Component Design: Fast Covariate Balancing with Synthetic Controls
Estimating the average effects of a binary treatment is one of the main goals of empirical economic and political studies. Randomization in controlled trials is one of the golden rules for estimating average treatment effects (ATE). Suppose the treatment assignment procedure guarantees that the potential outcomes are independent of the treatment status. In that case, a simple difference-in-mean (i.e., average outcomes of the treated and control units) estimator becomes an unbiased estimator of ATE. Nevertheless, a fully randomized experiment may be affected by a significantly high variance in the final estimation. Such variance can be reduced via exploiting the feature information in the observed data. Naturally, we focus on the following question: can the observed covariates improve the statistical properties of the ATE estimators via experimental designing? rubin2008objective,kasy2016experimenters. This problem is referred to as covariate balancing, which restricts the randomization to achieve covariate balance between treatment groups efron1971forcing,morgan2012rerandomization.
Covariate balancing has been substantially explored in the literature. One approach is applying covariate balancing via an adequately designed propensity score imai2014covariate,zhao2019covariate, which requires that we have access to a reasonably large sample of experimental units. However, large samples of experimental units are not always available in practice. Therefore, if only a moderate number of units are accessible to the treatment, bansal2018gram solve an NP-hard combinatorial optimization problem to balance the empirical covariates. To address this issue, we aim to design a computationally feasible optimal design of experimental studies.
We specifically consider the optimal design of experiments for the synthetic control estimator as in abadie2003economic,abadie2010synthetic,abadie2015comparative, which becomes an attractive estimation procedure when only a small number of units can be exposed to the experiment. That is, the experiment designer can observe the pre-treatment panel outcome data for a number of units in a number of time periods. Synthetic control compares treated units with a weighted average of untreated units. The weights are determined via empirical fit on the observed pre-treatment outcome. As an example, consider the study of Proposition 99's effect presented in abadie2010synthetic. Proposition 99 is a large-scale anti-smoking legislation program that California implemented in 1988. The policy maker wants to estimate the effect of this piece of legislation. Synthetic control suggests the policy maker to estimate the counterfactual outcomes after 1988 by using the observed outcomes of the states without legislative restrictions. abadie2010synthetic produces the synthetic control by a combination of Colorado, Connecticut, Montana, Nevada, and Utah and find out that annual per-capita cigarette sales in California were about 26 packs lower than what they would have been by the year 2000. Beyond this application, the method of synthetic controls has been used in many other empirical policy evaluation problems, including legalized prostitution cunningham2018decriminalizing, corporate political connections acemoglu2016value, taxation kleven2013taxation, to just name a few. To realize the benefits offered by synthetic control, doudchenko2021synthetic,abadie2021synthetic considers how an optimally designed experiment can help the experiment designer to further reduce the variance. doudchenko2021synthetic,abadie2021synthetic propose an optimization approach to select the control group based on the observed pre-treatment outcome. Namely, the choice of the treated units aims to balance the weighted average of treated and untreated covariates. As such, the designer can choose the best non-negative weights. doudchenko2021synthetic proves that the underlying optimization problem is still NP-hard and abadie2021synthetic relaxes the optimization problem into a canonical Quadratic Constraint Quadratic Program (QCQP). Nevertheless, the resulting QCQP is rather computationally demanding and applicable algorithms (e.g., SDR(cf. Semidefinite Programming) luo2010semidefinite) are not guaranteed to reach a global optimum.
In this paper, we aim to design the first globally convergent optimization algorithm for the weighted covariate balancing formulation doudchenko2021synthetic. Moreover, the algorithm is practically efficient. We achieve this by first removing the so-called units cardinality constraint being treated in doudchenko2021synthetic. Surprisingly, we find out that this relaxation can be shown to be equivalent to a phase synchronization problem singer2011angular. Although phase synchronization is still an NP-hard problem zhang2006complex, many practical algorithms have been recently developed for this nonconvex problem, bandeira2017tightness,boumal2016nonconvex,liu2017estimation,zhong2018near. Moreover, as we will argue, this nonconvex problem is polynomial-time solvable under a suitable data generating process (i.e. average case complexity). Motivated by this line of research, we propose Synthetic principal component Design (SPCD), which optimizes the treatment decision via (a normalized variate of) the generalized power method with spectral initialization chen2021spectral. If we assume that the pre-trement data is sampled from the linear fixed-effect model studied in abadie2010synthetic,xu2017generalized,athey2021matrix,ferman2021properties and invoke the realizable assumptions considered in abadie2021synthetic, we can further establish its global optimization guarantee and statistical estimation guarantee for the proposed synthetic control procedure. To the best of our knowledge, this is the first computational paradigm of combinatorial optimization-based experiment design which enjoys a global optimization guarantee.
\paragraph{Paper organization} We organize our paper as follows. In Section (ref), we introduce the setup in doudchenko2021synthetic,abadie2021synthetic where the authors consider covariate balancing under the synthetic control setting and reformulate it as a phase synchronization problem. In Section (ref), we introduce the generalized power method and provide global a optimization guarantee under the linear factor model abadie2003economic,xu2017generalized,athey2021matrix,ferman2021properties with a realizable assumption abadie2021synthetic. In Section (ref), we apply our method to both simulated and real-world data sets. We end with some closing remarks in Section (ref).
\paragraph{Synthetic Control} Synthetic control abadie2003economic,abadie2010synthetic is one of the leading methods to estimate the causal effect of a binary treatment in the panel data setting. It constructs a weighted combination of untreated groups used as controls, to which the outcome of the treatment group is compared. This construction of the control group significantly improves the performance when there are limited units that can be exposed to the experiment designer. Statistical consistency and inference properties have been established for the synthetic control under the linear factor models xu2017generalized,athey2021matrix,ferman2021properties, which is also the starting point of our theory.
\paragraph{Covariate Balancing} To make the difference-in-mean estimator more precise, experimenters sometimes restrict the treatment assignment to achieve covariate balance between treatment groups, i.e., \[ \min_{\{D_i\}_{i=1}^n} \left |\left|\frac{1}{\sum \limits_{i}D_i}\sum_{i:D_i = 1}X_{i} -\frac{1}{\sum \limits_{i}(1-D_i)}\sum_{i:D_i = 0}X_{i}\right |\right|^2 \] where $\{X_i\in\mathbb{R}^n\}_{i=1}^n$ are the observed units and $\{D_i\in\{0,1\}\}_{i=1}^n$ denotes the treatment experiments the designer aims to optimize efron1971forcing,morgan2012rerandomization,imai2014covariate,harshaw2019balancing. Indeed, this is a 0-1 NP-hard combinatorial optimization problem. Various different methods to approximately solve the problem have been proposed, including rerandomization morgan2012rerandomization,kallus2018optimal,kallus2021optimality,li2018asymptotic, design particular propensity score imai2014covariate,zhao2019covariate and the recently proposed Gram-Schmidt walk bansal2019algorithm,harshaw2019balancing. In this paper, we follow the setting in doudchenko2021synthetic,abadie2021synthetic which extend the covariate balancing problem to the synthetic control setting. doudchenko2021synthetic has proved that this problem is still NP-hard and abadie2021synthetic further relaxes the optimization problem into a canonical form of QCQP.
\paragraph{Phase Synchronization} Phase synchronization aims to estimate $n$ unknown angles $\theta_i\in [0,2\pi],i\in[N]$ through noisy measurements of their offset $\theta_i-\theta_j$ mod $2\pi$ singer2011angular, which has received intense interests in areas such as time synchronization between distributed networks, giridhar2006distributed; ranking cucuringu2016sync; computer vision wang2013exact,martinec2007robust; and optics and inverse problems rubinstein2001reconstruction,alexeev2014phase,singer2011three. To globally address the non-convex optimization problem arising from the inverse problem, bandeira2017tightness provides the first global results for the SDR under certain data generating procedures. In this paper, we follow a line of research which takes advantage of generalized power methods boumal2016nonconvex,liu2017estimation,zhong2018near to solve the problem in a way that shares a similar global guarantee to that of the SDP relaxation, but is computationally much more amenable. It is worth noting that, although the optimization problem considered in our paper is the same as phase synchronization, the data generating process is different from the considered earlier and this leads to our problem-specific generalized power method.
In this section, we introduce the mathematical formulation of synthetic control (SC) and the optimal experiment design problem studied in doudchenko2021synthetic.
\paragraph{Problem setup} We aim to estimate the effect of a binary treatment under the panel data setting. Researchers have access to the outcome metric of interest $Y \in\mathbb{R}^{N\times T}$ for $N$ units during $T$ time periods. At time $T$, researchers are required to execute an experiment by assigning a binary treatment described by $D_i\in\{-1,1\}, i\in [N]$ based on the observed pre-treatment data. If $D_i=1$, then a treatment needs to be applied to unit $i$. After the treatment experiment, furthermore outcomes are observed for additional $S$ time periods $t= T+1,\cdots, T+S$. During this period, every unit $i\in[N]$ in each time period $t$ is associated with the following two random outcomes: $Y_{it}(-1)=\mu_{it}+e_{it}, \text{ and } Y_{it}(1)=Y_{it}(-1)+\tau_i,$ where $\mu_{it}$ is the base outcome, ${\tau}$ is the treatment effect aiming to estimate and $e_{it}$ is the zero mean i.i.d idiosyncratic noise with variance $\text{Var}(\epsilon_{it})=\sigma$. Once treatment $D_i$ is applied, the experimenter is able to realize $Y_{it}=\frac{(D_i+1)}{2}Y_{it}(1)+\frac{(1-D_i)}{2}Y_{it}(-1)$.
Estimating the treatment effect $\tau$ is quite challenging because once we implement a treatment on unit $j$ (i.e., $D_j=1$) and observe the outcome $Y_{j,T+1}(1)$, then counterfactual outcome $Y_{j,T+1}(-1)$ is not observable. With the pre-treatment observation $Y_{iT}$, synthetic control literature abadie2010synthetic,xu2017generalized constructs the counterfactual estimate for a treated unit $j$ from a weighted average of other units' outcomes: $\hat Y_{j,T+1}(-1)=\sum_{i:D_i=-1} w_iY_{i,T+1}$. The weights $w_i$ are learned from the pre-treatment observed data via minimizing $\sum_{t=1}^T(Y_{jt}-\sum_{i:D_i=0}w_iY_{it})^2$. Then the treatment effect of unit $j$ we estimate can be written as $\tau_j=Y_{j,T+1}-\hat Y_{j,T+1}(-1)$.
In this subsection, we consider the synthetic design objective function proposed for two-way fixed effect in doudchenko2021synthetic,abadie2021synthetic and reveal its hidden connection with the phase synchronization problem. doudchenko2021synthetic,abadie2021synthetic aim to design treatment assignments $\{D_i=\pm 1\}_{i=1}^N$ and weights $\{w_i\ge0\}_{i=1}^N$ for outcome experiments at time $T+1$. For we aims to estimate a two-way fixed effect where the treatment effects are homogeneous, we can consider the weighted average treatment effect on the treated (wATET) $\tau=\sum_{i=1}^N \frac{D_i+1}{2}w_i\tau_i$ instead bottmer2021design. wATET can be estimated as a difference in weighted means estimator $\hat\tau=\sum_{i:D_i=1}w_i Y_{i,T+1}-\sum_{i:D_i=-1}w_i Y_{i,T+1}$ with $\sum_{i:D_i=1}w_i=\sum_{i:D_i=-1}w_i=1$. Following doudchenko2021synthetic, upon the outcome model, the mean squared error of the difference-in-weighted-means estimator admits the decomposition $$ \mathbb{E}\left[(\hat \tau-\tau)^2|\{D_i,w_i\}_{i=1}^N\right]=\underbrace{\left(\sum_{i:D_i=1} w_i\mu_{i,T+1}-\sum_{i:D_i=-1}w_i\mu_{i,T+1}\right)^2}_{\text{weighted covariate balancing}}+\sigma\sum_{i=1}^N w_i^2. $$ The designer aims to design the experiment with a lowest expected mean square error. Thus, doudchenko2021synthetic proposed the following mixed-integer programming for experimenting design with Synthetic Control:
By making a further simplification of the problem (ref), we introduce a change of variable $W_i=w_iD_i$. For $w_i\ge 0$, then $D_i=\textnormal{sgn}(W_i)$ and $w_i=|W_i|$. At the same time, the constraint $\sum_{i:D_i=1} w_i=\sum_{i:D_i=-1} w_i=1$ is equivalent to $\mathbbm{1}^\top W=0$ and the objective function \[ \frac{1}{T}\sum_{t=1}^T\left(\sum_{i:D_i=1} w_iY_{it}-\sum_{i:D_i=-1}w_iY_{it}\right)^2+\sigma \sum_{i=1}^Nw_i^2 \] can be reformulated as $W^\top (YY^\top+\lambda I) W$, where $W=[w_1,\cdots,w_N]^\top$ and $\mathbbm{1}\in\mathbb{R}^N$ is the all one vector. Thus, (ref) can be recast into
Although the reformulation (ref) translates the problem into a compact matrix form, it is still a nonconvex problem due to the constraint $||W||_1=1$. To deal with the constraint $\mathbbm{1}^\top W=0$, we add an extra term $\lambda(\mathbbm{1}^\top W)^2$ to the objective function, where $\lambda$ is a pre-defined hyper-parameter. Although this penalty method cannot produce the exact global solution, we can still recover the sign of the global solution (see Theorem (ref)). Once the sign of the global solution is identified, the remaining effort of computing the magnitude reduces to solving a convex problem ((ref)).
The following theorem states that the problem is equivalent to another well-known NP-hard non-convex problem --- Phase Synchronization singer2011angular.
The proof of Theorem (ref) and Theorem (ref) is omitted in the main text due to page limit and is shown in Appendix (ref).
In this section, we propose the generalized power method with spectral initialization to solve our problem. The method is inspired by the lines of its early efforts for phase synchronization problems journee2010generalized,boumal2016nonconvex,liu2017estimation,zhong2018near.
Spectral relaxation singer2011angular is the first simple and efficient approach to solve the phase synchronization problem. singer2011angular relaxed the $N$ constraints $|x_i|=1,i\in[N]$ to $||x||_2^2=n$. Then the solution becomes the leading eigenvector. liu2017estimation,zhong2018near showed that the eigenvector estimator is almost close to the global optima under certain data generating process. Following these works, we take our initial guess of the optimal experiment to be $\textnormal{sgn}(v)$, where $v$ is the smallest eigenvector of the matrix $(Y Y^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)$ with $\alpha,\lambda>0$ as two pre-defined hyper-parameters.
To further improve the experiment assignment, we utilize the generalized power method journee2010generalized,luss2013conditional, which considers the linearization of the objective function at the current point and moves towards a minimizer of this linear function over the non-convex set $\mathcal{C}$. It is also referred as to Frank-Wolfe algorithm for non-convex problems in the literature. Thus, the generalized power method enjoys monotonic improvement (boumal2016nonconvex). When $g(x)$ is a quadratic function taking the form as $x^\top A x$, the method will become similar to the power methods for the eigenvector problem. The only difference is the normalization step. The power method normalizes the whole vector but the generalized power method normalizes the individual entries. The generalized power method can be also understood as projected gradient descent smith1994optimization. Indeed, the update \[ y^{t+1}=\textnormal{sgn}[((\frac{1}{\beta} YY^\top+\frac{\sigma}{\beta} I+\frac{\alpha}{\beta } \mathbbm{1}\mathbbm{1}^\top)^{-1}+I) y^{t}] \] can be understood as a projection step ($\textnormal{sgn}$) after a gradient descent update with step size $\frac{1}{\beta}$: \[ y^{t+1}=\textnormal{sgn}(y^{t}+\frac{1}{\beta}((YY^\top+\sigma I+\alpha \mathbbm{1}\mathbbm{1}^\top)^{-1})y^{t}). \] Thus the algorithm shares a sufficient ascent condition for each iteration (shown in Appendix (ref)). Our algorithm is called Synthetic principal component Design (SPCD) and it is summarized in Algorithm (ref).
\paragraph{Normalized Variant} In boumal2016nonconvex,zhong2018near,liu2017estimation, the global optimality result is highly dependent on the assumption that the top eigenvector of the iteration matrix lies in $\{-1,1\}^N$. However, in our setup, the top eigenvector of the iteration matrix is the smallest eigenvector of the covariance matrix which may not be a $\{-1,1\}^N$ vector. This is also the case appearing in the phase retrieval candes2015phase,chen2015solving and the degree corrected stochastic block model zhao2012consistency,jin2015fast. Inspired by the SCORE jin2015fast,jin2022improvements method for degree corrected stochastic block models, we further introduce a normalization step to the generalization power method and call the new algorithm Normalized SPCD (cf. NormSPCD, see (ref) for details). We use the diagonal component of the inverse covariance as an estimate of the true normalization component. In Appendix (ref), we show that NormSPCD can be interpreted as a Riemannian gradient descent with a specific metric. Empirical results show that it is better than the original GPW in Figure (ref). This normalization technique may be of independent interest in other applications.
In this subsection, we provide the global optimization guarantee for the (normalized) generalized power method. bandeira2017tightness,boumal2016nonconvex,liu2017estimation,zhong2018near have shown that phase retrieval is globally solvable under certain generative models. We will show that generalized power method can globally converge under certain data generating processes, which are quite different from the ones assumed in the previous works. Following abadie2021synthetic, we consider a realizable linear factor model (also referred to as "interactive fixed-effects model") abadie2010synthetic,xu2017generalized,athey2021matrix,ferman2021properties, which has already been commonly employed in the literature as a benchmark model to analyze the properties of synthetic control estimators amjad2018robust,li2020statistical. Recently, shi2021assumptions justify the linear assumption from an independent causal mechanism viewpoint. The linear latent factor model is stated in the following assumption. .
To obtain the global optimality result, we further make the following realizable assumption that there is only one realizable experiment (zero error experiment) in population.
Under Assumptions (ref) and (ref), we can show the following global optimality result and the proof is shown in the Appendix (ref).
This section report the numerical tests of our algorithm. In subsection (ref), we validate our algorithm on the latent factor model. In subsection (ref), we demonstrate our algorithm on two real world datasets. Both experiments have shown the effectiveness of our proposed experiment design algorithm in terms of the the root-mean-square error (RMSE), where the squared differences between the true values of the treatment effects and the respective estimates are computed for each treatment period and averaged.
We first introduce a simplified implementation of (Norm)SPCD, which although not guaranteed optimum but efficient, simple and effective in practice. In the simplified implementation, we don't solve the convex program ((ref)) exactly, but using $w=\frac{2(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast}{||(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast||_1}$ to approximate instead. From (ref), we know that once the optimal design profile $y^\ast$ is obtained, then $w$ is the optimal design weight. Notice that we don't exactly globally solve the problem ((ref)) in the simplified implementation, although we obtained the right experiment profile $y^\ast$ (Theorem (ref)). The weight we obtained here is the solve of the penalized approximation, but empirically it works good. The whole process is described in Algorithm (ref). In all the experiment in this paper, we use this simplified implementation.
We generate data from the linear factor model (also referred to as an "interactive fixed-effects model", xu2017generalized,athey2021matrix). The outcome $Y$ comes from $$Y_{it}= v_t^\top\gamma_i+\tau \frac{D_{it}+1}{2}+e_{it},\,\forall i\in[N],\forall t\in[T+S].$$ where $\gamma_i$ is a vector of latent unit factor of dimension $L$ generated as a standard Gaussian. We simulate the treatment effect $\tau$ with ground turth $1$. The idiosyncratic noise is sampled from Normal $(0,\sigma^2)$ with $\sigma=1$. Fixing the test time period $S=10$ and the number of units $N=10$, we simulate different pairs of $(L,T)$ selection. We simulated three different choices of the time factor vector $v_t$
Synthetic principal component design (SPCD) selects the treated units whose features are representative of the whole aggregate market of interest abadie2021synthetic. At the same time, SPCD surpasses SC in every setting in Figure (ref), (ref) and (ref).
To exam our algorithm on real data, we follow arkhangelsky2019synthetic,doudchenko2021synthetic and utilize the US Bureau of Labor Statistics and the Anti-Smoking Legislation data to examine the validity of our algorithm. Besides synthetic control (SC), which randomly selects one unit to implement the treatment, we also implement an additional random baseline, which randomly select units as control/group with probability $1/2$. Through estimating the average treatment effect on the treated, we compare SC and the random assignment baseline with SPCD in terms of the RMSE. The final result is shown in Table (ref). SPCD surpasses SC a large margin on both of the datasets.
The Abadie–Diamond–Hainmueller Smoking Data. abadie2010synthetic uses SC to study the effects of Proposition 99, a large-scale anti-smoking legislation program that California implemented in 1988. \footnote{The Abadie–Diamond–Hainmueller Smoking data is first organized by abadie2010synthetic. In this paper, we use the organized version by arkhangelsky2019synthetic at \url{https://github.com/synth-inference/synthdid/blob/master/data/california_prop99.csv} which drop the data of minimum wage laws, gun laws to abortion laws in the original data and only considers the smoking outcome data.} To simulate the bias of SC and SPCD on this application, following athey2021matrix, we consider observations for 38 states (excluding California due to Proposition 99) from 1970 through 2000. We regard the first $T$ year as pre-treatment periods to produce the design and use the last $31-T$ years as post-treatment periods to test the performance of the treatment assignment. The final result is shown in Table (ref) and Figure (ref). Our design surpasses the random design by a large margin on most of the selection of time $T$. We also compare our method with the rerandomization design morgan2012rerandomization,kallus2018optimal,li2018asymptotic in Figure (ref), which shows that our algorithm is still better than 500000 times of rerandomization.
One typical design produced by our algorithm is shown in Figure (ref). The experiment design for different pre-treatment length $T$ is shown in Figure (ref). The plots show that our selection of the control group is robust to different pre-treatment time period and has the ability to represent all different geographic, demographic, racial, and social structure of states in the United State.
US Bureau of Labor Statistics. We also apply our algorithm on the unemployment rate of 50 states in 40 months from the US Bureau of Labor Statistics (BLS). \footnote{{The BLS Statistics data is available from the BLS website. In this paper, we use the organized version by arkhangelsky2019synthetic at \url{https://github.com/synth-inference/synthdid/blob/master/experiments/bdm/data/urate_cps.csv}. We thank arkhangelsky2019synthetic's authors carefully organize the data and open source it on github.} } We run 50 simulations such that each simulation utilizes a 20-by-$T+S$ matrix sampled from the original 50-by-40 dataset. More specifically, we randomly select 20 units and use the first $T$ time period to select the synthetic design and synthetic weight. The remaining $S$ time periods are the consecutive months that follow. In our experiment, we fix $S=5$ and run both experiment for $T=5,10$. The final result in terms of the RMSE is shown in Table (ref).
In this paper, we consider the optimal experiment design problem for synthetic control, where an NP-hard weighted covariate balancing problem is needed to be solve. Surprisingly, we reformulate the problem as a phase synchronization problem and propose a fast spectral initialized (normalized) generalized power method to address the resulting optimization problem efficiently. In face of a realizable linear factor model, we provide the first global optimization results for experiment design. Empirically, our method surpasses the original Synthetic Control and random design a large margin in terms of RMSE on various datasets. Additionally, the newly proposed normalization incorporated in GPM may have separate applications in degree correlated stochastic block models.
\printbibliography