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.
46,042 characters · 11 sections · 40 citation commands
Heterogeneous Treatment Effects in Panel Data
Suppose we observe an outcome of interest across $n$ distinct units over $T$ time periods; such data is commonly referred to as panel data. Each unit may have been subject to an intervention during certain time periods that influences the outcome. As a concrete example, a unit could be a geographic region affected by a new economic policy, or an individual consumer or store influenced by a marketing promotion. Our goal is to understand the impact of this intervention on the outcome. Since the effect of the intervention might vary across individual units and time periods as a function of the unit-level and time-varying covariates, we aim to estimate the heterogeneous treatment effects. This is a key problem in econometrics and causal inference, enabling policymakers or business owners to make more informed decisions about which units to target for future interventions.
We develop a new method, the Panel Clustering Estimator (PaCE), for estimating heterogeneous treatment effects in panel data under general treatment patterns. The causal effects are modeled as a non-parametric function of the covariates of the units, which may vary over time. To estimate heterogeneous effects, PaCE splits the observations into disjoint clusters using a regression tree and estimates the average treatment effect of each cluster. In addition to the development of this novel method, we make the following theoretical and empirical contributions.
\paragraph{Regression tree with bounded bias ((ref))} We show that, subject to mild assumptions on the regression tree and a density assumption on the covariates, the bias introduced by approximating the treatment effect function with a piece-wise constant function, produced by the regression tree, decreases polynomially as the number of leaves in the tree increases.
\paragraph{Optimal rate for average treatment effect estimation ((ref))} We obtain guarantees for the recovery of average treatment effects for each cluster of observations that matches the optimal rate, up to logarithmic factors. This extends the convergence guarantee of farias2021learning to a setting with multiple treatments. Similar to farias2021learning, we allow for general intervention patterns.
\paragraph{Enhanced empirical accuracy with simple estimator} Using semi-synthetic data, we demonstrate that PaCE achieves empirical performance competitive with, and often superior to, alternative methods for heterogeneous treatment effect estimation, including double machine learning, causal forest, and matrix completion-based methods. Furthermore, PaCE offers a much simpler and more interpretable estimate, as we achieve this enhanced performance using a regression tree with at most 40 leaves.
Causal inference using panel data is often approached through the synthetic control framework, where treated units are compared to a synthetic control group constructed from untreated units abadie2010synthetic, abadie2003economic. However, for these methods, the treatment is restricted to a block. Matrix completion approaches offer a more flexible alternative, allowing for general treatment patterns matrix_completion, bai2021matrix, xiong2023large. The treated entries are viewed as missing entries, and low-rank matrix completion is utilized to estimate counterfactual outcomes for the treated entries. In most of the matrix completion literature, the set of missing or treated entries are assumed to be generated randomly chen2019inference. However, this assumption is not realistic in our panel data setting. Other works explore matrix completion under deterministic sampling chatterjee2020deterministic, farias2021learning, foucart2020weighted, klopp2017robust, liu2017new, and our work extends the methodology from farias2021learning.
Several methods for estimating heterogeneous treatment effects use supervised learning techniques to fit models for the outcome or treatment assignment mechanisms chernozhukov2017double, chernozhukov2017orthogonal, DR, XLearner, nie2021quasi. Although these methods are intended for cases where observations are drawn i.i.d. from one or multiple predefined distributions, they could still be applied to panel data. In contrast, our approach is specifically designed for panel data with underlying structure. By avoiding complex machine learning methods, our approach not only provides accurate results for panel data but also offers interpretable estimates.
The first step of our method involves grouping data into clusters that vary in their treatment effects using a regression tree. As in breiman2001random, we greedily search for good splits, iteratively improving the tree's representation of the true heterogeneous treatment effects. However, unlike supervised learning where the splitting criterion is directly computable, treatment effects are not readily observed. Therefore, we estimate the splitting criterion, an approach used in athey2016recursive, athey2019generalized, su2009subgroup, wager2018estimation. While methods like causal forest only use information within a given leaf to estimate the treatment effect for observations in that leaf wager2018estimation, our method leverages global information from the panel structure to determine the average treatment effect of each cluster, enabling us to achieve high accuracy with just one tree.
The second step involves estimating the average treatment effect in each cluster. For this step, we utilize the de-biased convex estimator introduced by farias2021learning. This estimator is applicable to panel data with general treatment patterns but is specifically designed for estimating average, rather than heterogeneous, treatment effects. It uses a convex estimator, similar to those used in gobillon2016regional, xiong2023optimal, and introduces a novel de-biasing technique that improves upon the initial estimates. Additionally, farias2021learning builds on previous works chen2020noisy, chen2021bridging to provide theoretical guarantees for recovering the average treatment effect with provably minimal assumptions on the intervention patterns. The results of farias2021learning are specific to the case where there is only one treatment. We build on their proof techniques to show convergence when the number of treatments is $O(\log n)$, and also nest this technique into the first step of our algorithm to estimate heterogeneous effects.
For any matrix $A$, $\left\lVertA\right\rVert$ denotes the spectral norm, $\left\lVertA\right\rVert_{\mathrm{F}}$ denotes the Frobenius norm, $\left\lVertA\right\rVert_\star$ denotes the nuclear norm, and $\|A\|_{\text{sum}}$ denotes the sum of the absolute values of the entries in $A$. For any vector $a$, $\left\lVerta\right\rVert$ denotes its Euclidean norm, and $\left\lVerta\right\rVert_{\infty}$ denotes the maximum absolute value of its components. We use $\circ$ to denote element-wise matrix multiplication. For a matrix subspace $\mathbf T$, let $P_{{\mathbf T}^{\perp}}(\cdot)$ denote the projection operator onto ${\mathbf T}^{\perp}$, the orthogonal space of $\mathbf T$. We define the sum of two subspaces $\mathbf T_1$ and $\mathbf T_2$ as $\mathbf T_1+ \mathbf T_2 = \left\{ T_1+T_2\mid T_1\in\mathbf T_1,T_2\in\mathbf T_2 \right\}$. We write $a\lesssim b$ whenever $a \leq c\cdot b$, for some fixed constant $c$; $a\gtrsim b$ is defined similarly.
PaCE is designed for panel data, where an outcome of interest is observed across $n$ distinct units for $T$ time periods. Let $O\in \mathbb R^{n\times T}$ be the matrix of these observed outcomes. Our objective is to discern how these outcomes were influenced by various treatments under consideration.
We assume that, in a hypothetical scenario where the treatments were not applied, the expected outcomes would be represented by an unknown low-rank matrix $M^{\star}\in \mathbb R^{n\times T}$. To incorporate the influence of the external factors, we introduce the covariate tensor $\mathbf X:=[X_1,\dots,X_p] \in \mathbb R^{p\times n\times T}$, where an element $X_{izt}$ signifies the $i$-th covariate of unit $z$ at time $t$. For convenience, we will write $X^{zt} := [X_{1zt},\dots,X_{pzt}]$ as the vector of relevant covariates for unit $z$ at time $t$. We consider $q$ distinct treatments, each represented by a binary matrix $W_1,\dots,W_q$ that encodes which observations are subjected to that treatment. For each treatment $i\in[q]$ applied to unit $z$ at time period $t$, the treatment effect is modeled as a non-parametric Lipschitz function of the covariates $\mathcal T_i^{\star}(X^{zt})$. We consider the treatment effects of distinct treatments to be additive.\footnote{This is a flexible framework because, as we allow for up to $\log(n)$ treatments, the interactions of individual treatments can be modeled as additional treatments.} Combining all of these elements, the observed outcome matrix $O$ can be expressed as the sum of $M^{\star}$, the combined effects of the treatments, and a noise matrix $E$. Mathematically,
where $\mathcal T_i^{\star}(\mathbf X)$ denotes the matrix where the element in $z$-th row and $t$-th column is $\mathcal T_i^{\star}(X^{zt})$. We note that our model implicitly assumes unconfoundedness, which is common in the literature imbens2015causal.
Our goal is to estimate the treatment effect function $\mathcal T^{\star}_i$ for each $i\in[q]$, having observed $O$, $\mathbf X$, and $W_1,\ldots,W_q$. Our methodology PaCE consists of two main steps. First, we build a specialized tree for each treatment $i\in[q]$, designed to partition the observations into disjoint clusters that differ in their treatment effects. Next, we estimate the average treatment effect of each cluster, thus obtaining estimates for the treatment effect functions that are piece-wise constant functions of the covariates. The following subsections will delve into the specifics of these two steps.
Initially, the entire dataset with $n\times T$ entries is in one single cluster. This cluster serves as the starting point for the iterative partitioning process, where we split the data into $\ell$ more granular clusters for each treatment that better represent the true heterogeneous treatment effects. We present our method for clustering observations in (ref). Below, we discuss its two main steps.
\paragraph{Step 1} We solve (the convex) Problem (ref) to obtain a rough estimate of the counterfactual matrix, represented by $\hat M + \hat m\mathbf 1^\top$, for a given clustering of observations. The regularization term $\left\lVertM\right\rVert_\star$ penalizes the rank of $M$. The term $m\mathbf 1^\top$ is introduced to center each row of matrix $\hat M$ without penalization from the nuclear norm term. The inclusion of this term deviates from the convex problem in farias2021learning, and we find that its inclusion improves empirical performance. The choice of $\lambda$ is discussed in our theoretical results.
\paragraph{Step 2} In the second step, we greedily choose a valid split to achieve a better approximation of the true heterogeneous treatment effects. A valid split is one that complies with predefined constraints, which will be discussed in our theoretical guarantees. While the provided algorithm selects splits that minimize the mean squared error (MSE), alternative objectives such as mean absolute error (MAE) may be used instead.
After partitioning all observations into $\ell$ clusters for each of the $q$ treatments, the next step in our algorithm is to estimate the average treatment effect within each cluster. This is done using an extension of the de-biased convex estimator introduced in farias2021learning. We first obtain an estimate of the average treatment effects by solving a convex optimization problem. This initial estimate is then refined through a de-biasing step.
For simplicity, we define separate treatment matrices for each cluster. Letting for $k := q\times\ell$, $\Tilde Z_1,\dots, \Tilde Z_k$ are defined as follows: $\Tilde Z_{(i-1)\times \ell + j} = W_i\circ C_j^i$ for $i\in[q], j\in[\ell]$. These are binary matrices that identify observations that are in a given cluster and receive a given treatment. We denote the true average treatment effect of the treated in each cluster $i\in[k]$ as $\tilde\tau^{\star}_i$, and let $\delta_i$ be the associated residual matrices, which represent the extent to which individual treatment effects differ from the average treatment effect. More precisely, for $i\in[k]$,
Finally, we define $\delta$ as the combined residual error, where $\delta := \sum_{i=1}^k \delta_i$, and $\hat{E}:= E + \delta$ is the total error.
We normalize all of these treatment matrices to have Frobenius norm 1. This standardization allows the subsequent de-biasing step to be less sensitive to the relative magnitudes of the treatment matrices. Define $Z_{i} := \Tilde Z_i/\|\Tilde Z_i\|_{\mathrm{F}}$. By scaling down the treatment matrices in this way, our estimates are proportionally scaled up. Similarly, we define a rescaled version of $\tilde \tau^{\star}$, given by $\tau^{\star}_i:= \tilde\tau^{\star}_i \|\Tilde Z_i\|_{\mathrm{F}}.$
First, we solve the following convex optimization problem with the normalized treatment matrices to obtain an estimate of $M^{\star}$ and $\tau^{\star}$:
Because we introduced the $m\mathbf 1^\top$ term in the convex formulation, which centers the rows of $\hat M$ without penalization from the nuclear norm term, $\hat M$ serves as our estimate of $M^{\star}$, projected to have zero-mean rows, while $\hat m$ serves as our estimate of the row means of $M^{\star}$. For notational convenience, let $P_{\mathbf 1}$ denote the projection onto $\left\{ \alpha \mathbf 1^\top \mid \alpha\in\mathbb R^n \right\}$. That is, for any matrix $A\in \mathbb{R}^{n\times T}$,
Hence $P_{{\mathbf 1}^\perp}(M^{\star})$ is the projection of $M^{\star}$ with zero-mean rows. Let $m^{\star} := \frac{M^{\star}\mathbf 1}{T}$ be the vector representing the row means of $M^{\star}$. Let $r$ denote the rank of $P_{{\mathbf 1}^\perp}(M^{\star}) $,\footnote{Note that we can show $r \leq \text{rank}(M^*)+1$, so $P_{{\mathbf 1}^\perp}(M^{\star})$ is a low-rank matrix as well.} and let $P_{{\mathbf 1}^\perp}(M^{\star}) = U^{\star}\Sigma^{\star}V^{\star\top}$ be its SVD, with $U^{\star} \in \mathbb R^{n\times r}, \Sigma^{\star} \in \mathbb R^{r\times r},V^{\star} \in \mathbb R^{T\times r}$. We denote by $\mathbf T^{\star}$ the span of the space $\left\{ \alpha \mathbf 1^\top \mid \alpha \in \mathbb R^n \right\}$ and the tangent space of $P_{{\mathbf 1}^\perp}(M^{\star})$ in the manifold of matrices with rank $r$. That is,
The closed form expression of the projection $P_{{\mathbf T}^{\star\perp}}(\cdot)$ is given in (ref).
The following lemma gives a decomposition of $\hat \tau - \tau^{\star}$ that suggests a method for de-biasing of $\hat \tau$. The difference between our error decomposition and the decomposition presented in farias2021learning stems from the fact that we introduced the $m\mathbf 1^\top$ term in (ref).
As noted in farias2021learning, since $D^{-1}\Delta^1$ depends solely on observed quantities, it can be subtracted from our estimate. Thus, we define our de-biased estimator $\tau^d$ as $\tau^d_i= \left( \hat{\tau}- D^{-1}\Delta^1 \right)_i/\|\tilde Z_i\|_{\mathrm{F}}$ for each $i\in[k]$. The rescaling by $\|\tilde Z_i\|_{\mathrm{F}}$ adjusts our estimate to approximate $\tilde\tau^{\star}$, rather than $\tau^{\star}.$
The de-biased $\tau^d$ serves as our final estimate. Specifically, if an observation's covariates $X^{zt}$ belong to the $j$-th cluster, then for each treatment $i\in[q]$, our estimate of $\mathcal T^{\star}_i(X^{zt})$ is given by $\hat{\mathcal T}_i(X^{zt})=\tau^d_{(i-1)\times \ell + j}$, which represents the average treatment effect of treatment $i$ associated with the cluster that $X^{zt}$ belongs to.
Our analysis determines the conditions for the convergence of our algorithm's estimates to the true treatment effect functions $\mathcal T_i^{\star}(\mathbf X)$ for each $i \in [q]$. Our proof incorporates ideas from farias2021learning, wager2018estimation.
The final treatment effect estimate given by PaCE has two sources of error. First, there is the approximation error, which stems from the approximation of the non-parametric functions $\mathcal T^{\star}_i$ by piece-wise constant functions. The second source of error, the estimation error, arises from the estimation of the average treatment effect of each cluster. Note that the deviation between our estimate and the true treatment effect is, at most, the sum of these two errors. We will analyze and bound these two different sources of error separately to demonstrate the convergence of our method.
We start by bounding the bias introduced by approximating $\mathcal T^{\star}_i$ by a piece-wise constant function. To guarantee consistency, we establish constraints on valid splits. We adopt the $\alpha$-regularity condition from Definition 4b of wager2018estimation. This condition requires that each split retains at least a fraction $\alpha \in (0,\frac{1}{2})$ of the available training examples and at least one treated observation on each side of the split. We further require that the depth of each leaf be on the order $\log \ell$ and that the trees are fair-split trees. That is, during the tree construction procedure in (ref), for any given node, if a covariate $j$ has not been used in the splits of its last $\pi p$ parent nodes, for some $\pi >1$, the next split must utilize covariate $j$. If multiple covariates meet this criterion, the choice between them can be made arbitrarily, based on any criterion.
Additionally, we assume that the proportion of treated observations with covariates within a given hyper-rectangle should be approximately proportional to the volume of the hyper-rectangle. This is a `coverage condition' that allows us to accurately estimate the heterogeneous treatment effect on the whole domain using the available observations. For the purposes of this analysis, we assume that all covariates are bounded and normalized to $[0,1]^p$.
(ref) is a broader condition than the requirement for the covariates to be i.i.d. across observations, as assumed in wager2018estimation. Our assumption only requires that the covariates maintain an even density. Indeed, (ref) in (ref) illustrates that while i.i.d. covariates satisfy the assumption with high probability, the assumption additionally accommodates (with high probability) scenarios where covariates are either constant over time and only vary across different units, or constant across different units and only vary over time.
In the following result, we demonstrate that, as the number of observations grows, a regression tree, subject to some regularity constraints, contains leaves that become increasingly homogeneous in terms of treatment effect.
Using the above theorem, we can make the following observation. For any $i\in[q]$, denote by $\tilde\tau^{\star,i}$ the regression tree function that is piece-wise constant on each cluster $C_j^i$ for $j\in[\ell]$, with a value of $\tilde\tau^{\star}_{(i-1)\ell+j}$ assigned to each cluster $j$. Then, for all covariate vectors $X$, we have $\left\lvert\mathcal T_i^{\star}(X) - \tilde\tau^{\star,i}(X) \right\rvert \leq \frac{2L_i\sqrt{ p }}{\ell^s}$ because the value assigned to each cluster by $\tilde\tau^{\star,i}$ is an average of treatment effects within that cluster. This demonstrates that it is possible to approximate the true treatment effect functions with a regression tree function with arbitrary precision. As the number of leaves $\ell$ in the tree $\mathbb T$ increases, the approximation error of the tree decreases polynomially. The theorem requires an upper bound on $\ell$, which grows with $\min(n,T)$ as the dataset expands, due to the definition of the margin or error $M$. For PaCE, we constrain $\ell$ to be on the order of $\log n$ (we require $k = O\left( \log n \right)$ in (ref)), thereby complying with the upper bound requirement.
We now bound the estimation error, showing that it converges as $\max(n, T)$ increases. For simplicity, we will assume that $n\gtrsim T$ so that we can write all necessary assumptions and statements in terms of $n$ directly. We stress, however, that even without $n\gtrsim T$, the statement of (ref) also holds by replacing all occurrences of $n$ in (ref), (ref), and (ref) with $\max(n,T)$. Proving this only requires changing $n$ to $\max(n,T)$ in the proofs accordingly.
Our main technical contribution is extending the convergence guarantee from farias2021learning to a setting with multiple treatment matrices and a convex formulation that centers the rows of $\hat M$ without penalization. We require a set of conditions that extend the assumptions made in farias2021learning to multiple treatment matrices.
Throughout the paper, we define the subspace $\mathbf Z := \text{span}\{Z_i,i\in[k]\} + \text{span}\{\alpha \mathbf 1^\top \mid \alpha \in \mathbb{R}^n\}$. We next define $D^{\star}\in \mathbb {R}^{k \times k}$ and $\Delta^{\star1} \in \mathbb {R}^{k}$ similarly to $D$ and $\Delta^1$ in (ref), except we replace all quantities associated with $\hat M$ with the corresponding quantity for $M^{\star}$. That is, $D^{\star}_{ij} = \left\langle P_{{\mathbf T}^{\star\perp}}(Z_i),P_{{\mathbf T}^{\star\perp}}(Z_j) \right\rangle$, and $\Delta^{\star1}_i = \left\langle Z_i, U^{\star} V^{\star\top} \right\rangle$ for $i,j\in[k]$. Now, we are ready to introduce the assumptions that we need to bound the estimation error.
(ref) and (ref) are the direct generalization of Assumptions 3(a) and 3(b) from farias2021learning to multiple treatment matrices. Proposition 2 of farias2021learning proves that their Assumptions 3(a) and 3(b) are nearly necessary for identifying $\tau^{\star}$. This is because it is impossible to recover $\tau^{\star}$ if $M^{\star} + \gamma Z$ were also a matrix of rank $r$ for some $\gamma \neq 0$. More generally, if $Z$ lies close to the tangent space of $M^{\star}$ in the space of rank-$r$ matrices, recovering $M^{\star}$ is significantly harder (the signal becomes second order and may be lost in the noise). In our setting with multiple treatment matrices, we further require that there is no collinearity in the treatment matrices; otherwise, we would not be able to distinguish the treatment effect associated with each treatment matrix. Mathematically, this is captured by (ref).
The assumptions involving $E$ are mild: they would be satisfied with probability at least $1-e^{-n}$ if $E$ has independent, zero-mean, sub-Gaussian entries with sub-Gaussian norm bounded by $\sigma$ (see (ref)), an assumption that is canonical in matrix completion literature farias2021learning. Both (ref) and (ref) were also used in farias2021learning. Because $\delta$ is a zero-mean matrix that is zero outside the support of $\sum_{i\in[k]}Z_i$, the assumption $\left\lVert\delta\right\rVert = O(\sigma\sqrt{n})$ is mild. In fact, it is satisfied with high probability if the entries of $\delta$ are sub-Gaussian with certain correlation patterns moon2016linear. Furthermore, we note that if there is only one treatment matrix $(k=1)$ or if all treatment matrices are disjoint, then (ref) is guaranteed to be satisfied. This is because, for each $i\in[k]$, $\delta_i$ is a zero-mean matrix that is zero outside of the support of $Z_i$. Consequently, in these scenarios, we simply have $\left\lvert\left\langle Z_i, \delta \right\rangle\right\rvert =\left\lvert\left\langle Z_i, \delta_i \right\rangle\right\rvert = 0$.
We denote by $\sigma_{\max}$ and $\sigma_{\min}$ the largest and smallest singular values of $P_{{\mathbf 1}^\perp}(M^{\star})$ respectively, and we let $\kappa := \sigma_{\max}/\sigma_{\min}$ be the condition number of $P_{{\mathbf 1}^\perp}(M^{\star})$.
Note that the bound becomes weaker as $\|\tilde Z_j\|_{\mathrm{F}}$ shrinks. In other words, more treated entries within a cluster results in a more accurate estimate for that cluster, as expected. Our rate in (ref) matches, up to $\log n$ and $r$ factors, the optimal error rate obtained in Theorem 1 of farias2021learning, which bounds the error in the setting with one treatment matrix. Our bound accumulated extra factors of $\log n$ because we allowed $k = O(\log n)$ treatment matrices, and we have an extra factor of $r$ because we assumed $\left\lVert\delta\right\rVert\lesssim \sigma\sqrt n$, instead of $\left\lVert\delta\right\rVert\lesssim \sigma\sqrt{n/r}.$ Proposition 1 from farias2021learning proves that their convergence rate is optimal (up to $\log n$ factors) in the “typical scenario" where $\sigma, \kappa, r = O(1)$ and $\sigma_{\min} = \Theta(n)$. We can conclude the same for our rate in (ref).
The convergence rate of PaCE is primarily determined by the slower asymptotic convergence rate of the approximation error. However, our result in (ref) is of independent interest, as it extends previous results on low-rank matrix recovery with a deterministic pattern of treated entries to allow for multiple treatment matrices.
To assess the accuracy of PaCE, we demonstrate its performance on semi-synthetic data. Using publicly available panel data as the baseline $M^{\star} + E$, we introduced an artificial treatment and added a synthetic treatment effect to treated entries to generate the outcome matrix $O$. Since the true heterogeneous treatment effects $\mathcal T^{\star}(\mathbf X)$ are known, we can verify the accuracy of various methods. Our results show that PaCE often surpasses alternative methods in accuracy, while using a tree with no more than 40 leaves to cluster observations. Thus, PaCE not only offers superior accuracy, but also provides a simple, interpretable solution.
\paragraph{Data} To demonstrate the effectiveness of our methodology, we use two publicly available U.S. economic datasets. The first source comprises monthly Supplemental Nutrition Assistance Program (SNAP) user counts by zip code in Massachusetts, spanning January 2017 to March 2023, obtained from the Department of Transitional Assistance (DTA) public records mass_dta_reports. SNAP is a federal aid program that provides food-purchasing assistance to low-income families usda_snap. The second data source comprises nine annual demographic and economic data fields for Massachusetts zip codes, U.S. counties, and U.S. states from 2005 to 2022, provided by the U.S. Census Bureau us_census_bureau_acs1_api.
We use this data to conduct experiments on panels of varying sizes. In the first set of experiments, we use the SNAP user count as the baseline, $M^{\star} + E$, with the zip code-level census data for Massachusetts as covariates. The panel size is 517 zip codes over 70 months. The second set of experiments uses the population below the poverty line (from the U.S. census data) in each state and U.S. territory as the baseline and the remaining census data fields as covariates, with a panel size of 52 states and territories over 18 years. The third set of experiments is similar to the second, but the unit is changed to counties, resulting in a panel size of 770 counties over 18 years. The time period is provided as an additional covariate for the methods that we benchmark against.
\paragraph{Treatment generation} To generate the treatment pattern, we varied two parameters: 1) the proportion $\alpha$ of units receiving the treatment $\alpha\in\left\{ 0.05, 0.25, 0.5, 0.75, 1.0 \right\}$ and 2) the functional form of the treatment---either adaptive or non-adaptive. In the non-adaptive case, a random $\alpha$ proportion of units receive the treatment for a random consecutive period of time. In the adaptive case, for each time period, we apply the treatment to the $\alpha/2$ proportion of units that show the largest absolute percentage change in outcome in the previous time period. The purpose of the adaptive treatment is to mimic public policies that target either low-performing or high-performing units.
The treatment effect is generated by randomly sampling two covariates and either adding or multiplying them, then normalizing the magnitude of the effects so that the average treatment effect is 20% of the average outcome. For the SNAP experiments, the treatment effect is added to the outcome to simulate an intervention increasing SNAP usage. For the state and county experiments, the treatment effect is subtracted from the outcome to represent an intervention that reduces poverty.
Each set of parameters is tested on 200 instances, with the treatment pattern and treatment effect being randomly regenerated for each instance.
\paragraph{Benchmark methods} We benchmark PaCE against the following methods: double machine learning (DML) chernozhukov2017double, chernozhukov2017orthogonal, doubly robust (DR) learner DR, XLearner XLearner, and multi-arm causal forest athey2019generalized.\footnote{Additional methods are shown in (ref). In the main text, we include the five best methods.} We implemented PaCE to split the observations into at most 40 clusters. The hyper-parameter $\lambda$ used in the convex programs (ref) and (ref) was tuned so that $\hat M$ had a rank of $r=6$, which was chosen arbitrarily and fixed across all experiments. During the splitting procedure in (ref), we utilize binary search on each covariate to find the best split. See (ref) for implementation details.
\paragraph{Evaluation} We represent the estimated treated effects by $\hat{\mathcal T}(\mathbf X)$, a matrix where the element in $z$-th row and $t$-th column is the estimate $\hat{\mathcal T}(X^{zt})$. Accuracy of the estimates is measured using the normalized mean absolute error (nMAE), calculated relative to the true treatment effect $\mathcal T^{\star}(\mathbf X)$:
In (ref), we present the proportion of instances where each method results in the lowest nMAE. We observe that PaCE performs best in the SNAP and state experiments, achieving the highest proportion of instances with the lowest nMAE. However, it often performs worse than DML and Causal Forest in the county experiments. We hypothesize that this is because the county census data, having many units relative to the number of time periods, `most resembles i.i.d. data' and has less underlying panel structure. However, the average nMAE achieved by PaCE on this data set is competitive with that of DML and Causal Forest, as shown in Table (ref). Given that Causal Forest and DML use hundreds of trees, PaCE, which uses only one regression tree, is still a great option if interpretability is desired.
In (ref), we show the average nMAE as we vary the proportion of units treated. We only show the results for the SNAP experiments with a non-adaptive treatment pattern, and we refer to (ref) and (ref) in (ref) for the results of additional experiments and for the standard deviations. We observe that while the performance of all methods improves as the proportion of treated entries increases, PaCE consistently achieves the best accuracy on this set of experiments.
This work introduces PaCE, a method designed for estimating heterogeneous treatment effects within panel data. Our technical contributions extend previous results on low-rank matrix recovery with a deterministic pattern of treated entries to allow for multiple treatments. We validate the ability of PaCE to obtain strong accuracy with a simple estimator on semi-synthetic datasets.
One interesting direction for future research is to investigate whether the estimates produced by PaCE are asymptotically normal around the true estimates, which would allow for the derivation of confidence intervals. Further computational testing could also be beneficial. Specifically, assessing the estimator's performance across a broader array of treatment patterns and datasets would provide deeper insights into its adaptability and limitations. Expanding these tests to include multiple treatments will help in understanding the estimator's effectiveness in more complex scenarios.