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.
41,654 characters · 5 sections · 29 citation commands
Optimal selection of the number of control units in kNN algorithm to estimate average treatment effects
\thispagestyle{empty}
\setcounter{page}{1} \pagenumbering{arabic}
Matching estimators are commonly used to estimate average treatment effects when treatment assignment is independent of outcome conditional on confounding variables Cameron05. However, it seems that optimal selection of the number of matching units ($k$) is an open question as most of the applied literature fixes $k$ small without any optimal criterion, often $k=1$ Abadie2016. Choice of $k$ is a trade-off between bias and variance. Small $k$ implies small bias, but high variance, and vice versa. So, the aim of this paper is to propose a simple approach to optimally select $k$. In particular, we focus on the $k$ nearest-neighbors (kNN) algorithm when used to estimate treatment effects: average treatment effect (ATE), and average treatment effect on treated (ATET).
Optimal selection of $k$ in kNN algorithm is done by cross validation when prediction or classification are the objectives Hastie2009. This is easily performed as the statistical object of interest is observable in-sample. On the other hand, average treatment effects are not observable anywhere. So, we follow similar arguments to Athey2015 proposing an unbiased function of observable variables to optimally choose $k$ through cross validation when the main inferential concern is treatment effects.
We apply our proposal to analyze ATET due to enrollment in 401(k) on amount and probability of positive accumulated net financial assets. We found significant ATET whose point estimates are approximately \$15 K and 19% given an optimal $k^*=19$ in both cases.
After this brief introduction, Section (ref) shows our proposal to optimally select $k$ in the kNN algorithm. Section (ref) displays results from some simulations exercises, and Section (ref) shows the results of our application regarding 401(k) enrollment. Section (ref) ends with some concluding remarks.
We will focus on the average treatment effect on treated as we will use this in our application. However, results associated with the average treatment effects are in the Appendix (see subsection (ref)). We build the matching units using $k$ nearest neighbors (kNN) with replacement, where $k$ is selected minimizing the mean squared error between a conditional unbiased estimator of the average treatment effect on the treated (ATET) and a matching estimator of the ATET.
We consider a framework where there is a random binary treatment ($D_i=\left\{0,1\right\}$) that is independent of the outcome variable ($Y_i$) conditional on independent observable variables ($\bm{X}_i\subset \mathbb{R}^P$), $i=1,2,\dots,N$. The average treatment effect on treated,
where $\bm{X}_i$ are control variables, $\left\{Y_i(0),Y_i(1)\right\}$ are the potential outcomes under different states, untreated and treated,
However, we do not simultaneously observe the same individual under both states. So, we propose to use the matching approach to build the synthetic controls, the general formula is Cameron05,
where $A_k^i(\bm{x})=\left\{j:||\bm{X}_i-\bm{X}_j||<||\bm{X}_i-\bm{X}_l||, j\in\left\{D_j=0\right\}, l\in\left\{D_l=0\right\}\right\}$ is the set of the $k$ closest untreated units to treated unit $i$ in terms of the Euclidean norm in covariates, $j,l=1,2,\dots,N_0$, $||\cdot ||$ is the Euclidean norm, $N_0$ is the sample size of untreated individuals ($D_i=0$), $N_1$ is the sample size of treated individuals ($D_i=1$), and $k$ is the cardinality of $A_k^i(\bm{x})$, that is, the number of controls for treated unit,
Therefore, the synthetic control for treated individual $i$ is built using an average of untreated nearest individuals (neighbors).
The choice of $k$ is a trade-off between bias and variance. Selecting just one neighbor minimizes bias, this is equal to zero if there is exact match ($k=1, \bm{x_i}=\bm{x_j}$), but implies high variability. On the other hand, a large amount of neighbors decreases variance but increases bias. In general, $\bm{x_i}=\bm{x_j}$ is an event of probability 0 for continuous covariates, then Abadie2006 show that matching estimators have an asymptotic bias. This asymptotic bias can be ignored if $N_1^{P_c/2}=O(N_0)$ where $P_c$ is the number of continuous covariates ($P_c/2>1$). The bias of the ATET is
where $\mu_0(\bm{x}_i)=\mathbb{E}\left[Y|\bm{X}=\bm{x},D=0\right]$. This can be estimated non-parametrically using a series expansion estimator to obtain $\hat{B}^{TM}$ Abadie2011.
Therefore by equations (ref) and (ref), the bias-corrected matching estimator is
Abadie2006 show under suitable assumptions (see Assumptions 1, 2', 3' and 4 in their paper) that
where $V^{E,T}=\frac{1}{N_1}\sum_{i=1}^N \left(D_i+(1-D_i)\frac{J_k^i}{k}\right)^2 \sigma^2(\bm{X}_i,D_i)$, $J_k^i=\sum_{j=1}^N\mathbbm{1}\left\{i\in A_k^j\right\}$ is the number of times unit $i$ is used as a match given $k$ matches per unit, $\sigma^2(\bm{X}_i,D_i)=\mathbb{V}(Y|\bm{X}_i=\bm{x}_i,D_i=d_i)$, which can be consistently estimated by $\hat{\sigma}^2(\bm{X}_i,D_i)=\frac{k}{k+1}\left(Y_i-\frac{1}{k}\sum_{l\in L_k^i}\right)$, $L_k^i(\bm{x})$ is the set of the $k$ closest units to unit $i$ in terms of the Euclidean norm in covariates such that have the same value for the treatment variable, and $V^{\Delta(\bm{X}),T}=\mathbb{E}\left[(\left[\mu_1(\bm{X}_i)-\mu_0(\bm{X}_i)\right]-{\Delta}^{T})^2|\bm{D}=1\right]$ is the variance of the conditional average treatment effect on treated.
In an ideal situation, we select $k$ in the regression nearest neighbor framework (kNN) minimizing the mean squared error of the average treatment effect on treated (equation (ref)), that is,
However, the average treatment effect on treated ($\Delta^{T}$) is not observed. Therefore, we follow similar arguments to Athey2015 proposing an unbiased function of potentially observable variables to define the ATET in program (ref).
We have the following assumptions:
We assume conditional independent and identically distributed random sample drawn from a population. This implicitly implies Assumption (ref) (SUTVA). Then, we define the individual treatment effect on treated
where $P\left[D=1\right]$ is the marginal probability of treatment, which is the same for every individual given random sampling, and $P(\bm{X}_i)=P\left[D_i=1|\bm{X}_i=\bm{x}\right]$ is the conditional probability of treatment given regressors $\bm{X}_i$, that is, the propensity score Rosenbaum1983.
We require estimators for $P\left[D=1\right]$ and $P(\bm{X}_i)$ to estimate $\Delta_i^{T*}$. We use $\frac{1}{N}\sum_{i=1}^N D_i$, $N=N_0+N_1$, as estimator for the marginal probability of treatment, and a logit model to estimate the conditional probability of treatment. Obviously, these estimators can be changed. Notice that this introduces an extra source of variability as we do not observe $P\left[D=1\right]$ and $P(\bm{X}_i)$ in observational data sets. This kind of issue is discussed by Athey2015 who point out that transformations like equation (ref) does not optimally use all available information, this may mean extra variability. Additionally, Abadie2016 show that ATE using matching based on the estimated propensity score is more efficient than matching based on the true propensity score. On the other hand, ATET using the estimated propensity score can be more o less efficient compared to the true propensity score. We performed some simulation exercises that show better inferential performance based on the estimated propensity score (available upon authors request). In addition, as our ATET (ATE) is asymptotically biased corrected, the mean squared error decomposition shows that variance explains the most part (75% on average), as expected. However, it remains a proportion to finite sample bias effect (25% on average).
Give that ${\Delta_i^{T*}}$ in equation (ref) is a conditional unbiased estimator for $\Delta^T$ (Proposition (ref)), we select $k$ in the kNN regression framework such that
where
and
where $G$ is the number of groups, $N_g$ is the sample size of each group, and $N_{1,g}$ is the sample size of the treated group in group $g$, $g=1,2,\dots,G$. Superscripts $Train$ and $Test$ refer to train and test data sets (see below).
Notice that $\widehat{\Delta_{i,g}^{T*}}$ is the sample version of the conditional unbiased estimator of the ATET $\Delta_i^{T*}$, so its average is a consistent estimator for ATET. On the other hand, $\hat{\Delta}_g^{BC}$ is bias-corrected version of the matching estimator Abadie2006.
kNN regression is a standard approach in the machine and statistical learning communities Hastie2009. So, we follow $k$-fold cross-validation to select the optimal $k$ Hastie2009. In particular, we randomly split the data set in $G$ roughly equal-sized groups such that we keep the proportion between treated and untreated individuals in each group, then we calculate $\frac{1}{N_g}\sum_{i=1}^{N_g}\widehat{\Delta_{i,g}^{T*}}$, $g=1,2,\dots,G$. We fit the model using $G-1$ groups, that is, we fix $k$, and identify the $k$ nearest neighbors in the control group at the training data set in terms of the Euclidean norm to individual $i$ in the treated group at the $G$-th left group (test data set), that is,
then we estimate $\frac{1}{k}\sum_{j\in A_k^i(\bm{x}_g^{Test})}Y_{j,g}^{Train}(0)$, and obtain $\hat{\Delta}_{i,g}^{TM}=Y_{i,g}^{Test}(1)-\frac{1}{k}\sum_{j\in A_k^i(\bm{x}_g^{Test})}Y_{j,g}^{Train}(0)$ and $\hat{\Delta}_{i,g}^{BC}=\hat{\mu}_0(\bm{x}_i)_{g}^{Test}-\frac{1}{k}\sum_{j\in A_k^i(\bm{x}_g^{Test})}\hat{\mu}_0(\bm{x}_j)_{g}^{Test}$ for $i\in\left\{D_{i,g}^{Test}=1\right\}$. So, we can calculate the mean squared error for the $G$-th left group. This is done for $g=1,2,\dots,G$ and different values of $k$ such that we average over the $G$ prediction errors given $k$, and select the $k$ which has the minimum average error (see program (ref)).
Given an optimal $k=k^*$, we use equation (ref) to estimate the average treatment effect on treated using the whole sample size. As we said, this estimator converges in distribution to a normal distribution, so we use the asymptotic result to build confidence intervals.
Algorithm (ref) summarizes our methodological proposal:
We set exactly the same simulation setting as Otsu2017 who proposed a weighted bootstrap to perform inference of matching estimators for ATE and ATET as it is well known that naive bootstrap produces invalid inference due to failing to reproduce the distribution of the number of times each unit is used as a match Abadie2008.
We show the simulation setting for exposition purposes:
where $\epsilon_i$, $v_i$, $\psi_i$ and $\zeta_i$ are mutually independent. The ATET as well as ATE is $\tau = 0.5$, $\gamma_1=0.15$, $\gamma=0.7$, $P=5$ and there are six different curves for $m_j(||\bm X_i||)$ (see Otsu2017 for details):
In addition, we perform other simulation exercises where the outcome variable is binary. In particular,
where $\epsilon_i\sim LG(0,1)$, $LG$ denotes a logistic distribution, $\beta=0.5$, $\tau= 0.5$, and $P(D_i=1)=P(\gamma_1+\gamma_2 ||\bm X_i||>v_i)=F_{LG}(\gamma_1+\gamma_2 ||\bm X_i||)$, $\gamma_1=0.15$, $\gamma_2=0.4$, $v_i\sim LG(0, 1)$ and $F_{LG}$ is the distribution function. Other components of the specification are as in the previous setting.
Observe that in this setting $ATE_i=\mathbbm{E}(Y_i(1)-Y_i(0)|\bm X_i=\bm x_i, Z_i = z_i)=P(Y_i(1)=1|\bm X_i=\bm x_i, Z_i = z_i)-P(Y_i(0)=0|\bm X_i=\bm x_i, Z_i = z_i)=F_{LG}(\tau + \beta m_j(||\bm X_i||))-F_{LG}(\beta m_j(||\bm X_i||))$, which has the same analytical expression for the ATET, but conditional on the relevant sample ($D_i=1$).
Our results are based on a sample size equal to 100 Otsu2017.\footnote{We perform other simulation exercises using larger sample sizes. Outcomes show similar results. Available upon authors request.} We have 1,000 replications for each simulation exercise.
First, we present the results for the average treatment effects on treated for the continuous outcome variable.\footnote{Results for the average treatment effects are in the Appendix, subsection (ref). Outcomes are similar to ATET.} Table (ref) displays the mean relative squared errors (MRSE),
where $\hat{\Delta}^{BC}_s(k)$ and $\hat{\Delta}^{BC}_s(k^*)$ are the ATET estimates using any $k$ and the optimal $k=k^*$, and $S_k$ is the number of times that $k$ was not optimal.
We see in Table (ref) that there are relative large advantages using our optimal $k^*$ as all figures are larger than 1.
Figure (ref) shows the average over 1,000 replications of the ratios between the length of the 95% confidence intervals associated with any number of neighbors ($k$, $x$ axis) and the optimal selection using our proposal. This is done for the six data generating settings ($m_j(||\bm X_i||), j=1,2,\dots 6$). We see in this figure that a small $k$ implies relative wide intervals, whereas large $k$ implies narrow intervals. This is expected as small $k$ implies high variability compared with large $k$. Observe that intervals lengths decreases very fast, and stabilize around $k=10$. On the other hand, we see in Figure (ref) that the type I error due to rejecting the null hypothesis $ATET=\tau$, given a significance level equal to 5%, is the lowest using $k=1$. Using $k$ minimizing the mean squared error gives a relatively small type I error compared with other $k$ selection.
Figure (ref) displays the histograms of the optimal $k$ in our six simulation settings. We see that optimal selection is dominated by small $k$. This could be evidence to support $k=1$, which is widely used in many applications. Although, it seems that $k=\left\{2,3\right\}$ are the most relevant, and large values of $k$ cannot be discarded, especially when using large sample sizes (we performed other simulations exercises that show this).
Results of the ATET for our simulation setting using a binary outcome are qualitatively similar to the continuous outcome. Table (ref) show the mean relative squared errors. Again, all figures are larger than 1, which means that using a non-optimal $k$ implies larger estimation errors.
Figure (ref) shows the relative 95% confidence intervals length ratios associated with $k$ and $k^*$. All six simulation settings look very similar with wider confidence intervals when $k$ is small, and narrower intervals when $k$ is large. There is stabilization around $k=10$. Figure (ref) shows type I errors using 5% as significance level. Larger $k$ implies larger errors, and using $k^*$ implies small type I errors in general, except when $k$ is very small. Small type I error can be explained by means of Figure (ref) where we see that optimal selection of $k$ is also based on small $k$. Here $k^*=\left\{2,3,4,5\right\}$ are particularly relevant.
In general, we can observe that our proposal implies relatively smaller estimation errors achieving a compromise between confidence interval length and type I error having as a result of small optimal $k$ most of the times.
To illustrate our proposal we estimate the average treatment effect on treated of 401(k) participation on accumulated net financial assets. 401(k) is a retirement plan where savings contributions are provided by employers deduced from employees' payment before taxes such that taxes on capital gains are avoided. Employees should achieve some eligibility criteria such that the plan does not have universal applicability in principle. The impact evaluation of 401(k) on net financial assets has been done previously by Benjamin2003,chernozhukov2004effects, Conley2012,chernozhukov2017double among others. We overcome the lack of random assignment following chernozhukov2017double arguments. Enrollment in a 401(k) plan is exogenous when controlling for observable variables which drive employment decisions, such as income, when the plan initially became available. Therefore, our dataset is from the 1991 Survey of Income and Program Participation composed by 9,915 observations.
We control for the same set of variables as chernozhukov2004effects. In particular, age (linear and squared), education (years), household size, income (seven levels), and binary variables indicating: defined benefit pension, home ownership, marital status (married), participation in IRA (individual retirement account) and two-earners status (both household heads contribute to household income). Details can be found in Benjamin2003.
Table (ref) shows unconditional mean difference tests of our outcome variables: accumulated net financial assets, and a binary variable indicating positive accumulated net financial assets. Employees enrolled in 401(k) plans have \$27.4 more accumulated net financial assets on average than employees who are not enrolled. The unconditional probability of positive net financial assets is 31% on average higher for the former. We can see that there are unconditional statistically significant differences.
Figure (ref) shows the ATET in the accumulated net financial assets using our optimal selection, $k^*=19$, and $k=1$, as a reference. We see similar ATET point estimates, which are approximately \$ 15K (similar outcomes are found by Conley2012 under the exogeneity assumption using a Bayesian approach for a local average treatment effects, LATE). However, the 95% confidence interval using $k^*$ is (\$10.9K,\$18.9K), whereas $k=1$ interval is (\$9.1K,\$19.9K), that is 35% larger. Results regarding the probability of positive accumulated net financial assets can be seen in Figure (ref). Both, $k^*=19$ and $k=1$ indicate statistically significant ATET whose point estimate is approximately 19.9% and 18.7% respectively. The 95% confidence intervals are (16.3%, 21.3%) and (17.7%, 21.3%) for $k=1$ and $k^*=19$. This means that using $k=1$ implies a 38.9% wider interval than using $k^*=19$. Observe that both exercises indicate the same optimal $k$. This means that diagnostic balance tests are based on the same control groups. Figure (ref) shows the love plot indicating that control groups based on $k^*=19$ satisfy mean balancing conditions.
We present a simple approach to obtain an optimal selection of the number of control units ($k$) to build the “synthetic” or “counterfactual” unit for each individual in the k nearest neighbors algorithm. Our approach is based on a simple unbiased reference estimator (individual treatment effect on treated for ATET and individual treatment effect for ATE) such that we select $k$ minimizing the mean squared error with respect to the reference estimator. Therefore, the optimal $k=k^*$ achieves a balance between confidence interval length and type I error.
Our application suggests that $k^*$ can be relatively large compared to the common practice of using $k=1$. This implies significant reductions regarding confidence interval lengths. Although, both choices show 5% statistically significant effects of 401(k) enrollment in the probability of positive accumulated net financial assets, and its amount.
A final remark is that using as reference estimator for treatments effects the one obtained using $k=1$ is an interesting alternative. We obtained same qualitative results, that is, a balance between interval length and type I error, and a considerable reduction in MRSE. Simulations exercises confirm this.