Estimation of Treatment Effects based on Kernel Matching

Matching is a crucial method for estimating causal effects, widely applied in various research fields, including policy evaluation (Sekhon,, 2004; Herron and Wand,, 2007), healthcare research (Rubin,, 1997; Christakis and Iwashyna,, 2003), and economics (Dehejia and Wahba,, 1999; Galiani et al.,, 2005; Abadie and Imbens,, 2006). The process of estimating causal effects fundamentally depends on imputing a potential outcome for an individual, calculated as a weighted average of the observed outcomes of their nearest neighbors from the relevant counterfactual set (Rosenbaum and Rubin,, 1983). Thus, causal effects are typically expressed as the difference between the imputed outcome and the actual observed outcome, similar to what would be expected from a randomized experiment. Various methods have been developed to assign matched subjects based on covariates, such as nearest available matching using propensity scores, Mahalanobis metric matching (Rosenbaum and Rubin,, 1985), full matching (Rosenbaum,, 1991; Hansen,, 2004), genetic matching (Diamond and Sekhon,, 2013), and matching using sufficient dimension reduction (Luo and Zhu,, 2020).

Among various matching methods, covariate matching (Abadie and Imbens,, 2006) can be challenging in high-dimensional settings, as it requires exact or near-exact matches across all covariates. This strict requirement often reduces the effective sample size, limiting the method’s reliability and suitability for high-dimensional data or complex distributions. Alternatively, propensity score matching (Abadie and Imbens,, 2016) depends on the correct specification of the propensity score model; misspecification can introduce bias and result in extreme weights. Beyond these matching methods, several alternative approaches have been proposed for causal effect estimation, including inverse probability weighting (IPW) (Hirano et al.,, 2003), subclassification (Rosenbaum and Rubin,, 1983), and doubly robust (DR) estimation (Bang and Robins,, 2005). The accuracy of the IPW estimator accuracy is highly sensitive to the correct specification of the propensity score model, with even minor misspecifications potentially leading to substantial errors due to its sensitivity to high variance from extreme weights (Imai and Ratkovic,, 2014). Doubly robust (DR) estimators, which combine propensity score and outcome regression models, can yield unbiased estimates if at least one model is correctly specified. However, DR estimators are sensitive to the accuracy of both models and may be complex to implement effectively.

In this paper, we address the above challenges by introducing a novel kernel-matching estimator and providing a detailed methodological framework, including discussions on the consistency and asymptotic normality of the estimator. Kernel matching, as outlined by Heckman et al., (1997) and Heckman et al., 1998b , is a non-parametric approach that applies kernel functions to weight control cases based on their similarity to treatment cases in the covariate space, facilitating the estimation of treatment effect (Handouyahia et al.,, 2013). The kernel matching method allows for increased flexibility and stability through adjustable bandwidth parameters. By balancing bias and variance through local smoothing, kernel matching produces more robust estimates than traditional one-to-one propensity score matching, especially in situations with limited control cases available for creating a matching cohort. It mitigates extreme bias or variance and reduces dependence on stringent model assumptions (Berg,, 2011).

Although more complex than basic matching techniques, kernel matching extends interval and nearest-neighbor matching by assigning weights to all control cases associated with each treatment case, giving greater emphasis to those with closer covariate values. Smith and Todd, (2005) offered an intuitive explanation of kernel matching and its generalizations, such as local linear (Heckman et al.,, 1997; Heckman et al., 1998a, ; Heckman et al., 1998b, ) and local quadratic matching (Ham et al.,, 2011). By smoothing over the covariate space, kernel matching reduces the risks associated with model misspecification, resulting in more reliable and unbiased treatment effect estimates. This approach is particularly advantageous in settings with high overlap between treatment and control cases, ensuring that comparisons are made between genuinely comparable groups.

When choosing a suitable matching method, selecting the appropriate function of the covariates is crucial as the function may significantly affect the asymptotic behaviors of the resulting estimators (Heckman et al.,, 1997; Abadie and Imbens,, 2006). Heuristically, Rosenbaum and Rubin, (1983) suggested that using the original covariates is preferred for its simplicity and for ensuring unbiased matching at the population level, provided that the ignorability assumption (i.e., the potential outcomes are independent of treatment after controlling for the specific set of covariates) and the overlap assumption (i.e., the range of data is similar across treatment groups), are satisfied. However, Abadie and Imbens, (2006) pointed out that the matching approach can be vulnerable to the “curse of dimensionality”, rendering matching ineffective when the dimension of covariates is high. To address this issue, a natural choice is using the propensity score, which is the probability of a subject being a treatment case given its covariates, proposed by Rosenbaum and Rubin, (1985). This reduction in dimensionality makes it easier to find similar treatment and control cases. Similar to the original covariates, the propensity score also ensures the unbiasedness of matching under the ignorability assumption and the overlap assumption. Mathematically, suppose $X$ is the vector of covariates, the propensity score $p(X)$ can be defined as

$$p(X)=F(X^{\top}\beta),$$

where $\beta$ is the parameter vector, $F$ is a link function that maps $X^{\top}\beta$ into a probability in $(0,1)$ (e.g., a logit or probit function). Nevertheless, the estimation of the propensity score commonly relies on parametric models, including the link function $F$, which are susceptible to model misspecification. In high-dimensional settings, kernel matching enhances propensity score matching by smoothing across similar observations. The kernel function averages treatment effects among units with similar propensity scores, stabilizing estimates and reducing variance. High-dimensional data often suffer from sparsity, making exact matches between treatment and control cases difficult to find. Kernel matching addresses this issue by weighting nearby units based on their propensity scores, allowing for effective comparisons even when exact matches are rare. The kernel function introduces a smoothing effect, which acts as a form of regularization. This regularization helps control overfitting to noise, a common issue in high-dimensional data, making the treatment effect estimation more robust. In high-dimensional spaces, irrelevant covariates can introduce noise into the matching process. By focusing on propensity scores, which encapsulate the relevant information from all covariates, kernel matching reduces the impact of irrelevant covariates, making the method less sensitive to noise. Therefore, kernel matching based on propensity scores is a powerful tool for addressing the challenges of high-dimensional settings, as it effectively manages the complexities associated with having many covariates.

In this paper, we propose estimating treatment effects using kernel matching based on $F(X^{\top}\hat{\beta})$ with observational data. Through a Monte Carlo simulation study with different settings, we demonstrate the practical application of kernel matching and highlight its advantages over traditional estimation methods. This study underscores the importance of adopting robust matching techniques in causal inference and showcases the potential of kernel matching to enhance the validity and reliability of treatment effect estimates. Compared to traditional nearest-neighbor matching, kernel matching does not require pre-specifying the number of matched individuals (e.g., one-to-one matching). Instead, kernel matching adaptively selects suitable matches through the kernel function. By using more control cases through weighting, kernel matching can reduce the variance of the estimate. We show that the kernel matching method is a powerful tool in the arsenal of methods for causal inference in observational studies.

The rest of the paper is organized as follows. In Section 2, after introducing the mathematical notation and the concepts of treatment effect estimation, the proposed kernel matching estimator is developed. In Section 3, the asymptotic properties of the kernel-matching estimator are established. A Monte Carlo simulation study is used to evaluate the performance of the proposed method in Section 4. Section 5 provides a practical data analysis to illustrate the proposed methodology. Finally, some concluding remarks are provided in Section 6.

Following the framework of Rubin, (1974), the treatment effect on the potential outcomes can be defined as follows. For unit $i$ ($i=1,\ldots,N$), let the variable $W_{i}$ indicates whether an active treatment was received $(W_{i}=1)$ or not received $(W_{i}=0)$, for $W_{i}\in\{0,1\}$. The observed outcome is $Y_{i}=W_{i}Y_{i}(1)+(1-W_{i})Y_{i}(0)$, where $Y_{i}(1)$ and $Y_{i}(0)$ indicate the potential outcomes that the individual would have received active treatment or not, respectively. In other words, for unit $i$, the observed outcome $Y_{i}$ for treatment $W_{i}$ can be expressed as

$$Y_{i}=\begin{cases}Y_{i}(0)&\text{ if }W_{i}=0,\\ Y_{i}(1)&\text{ if }W_{i}=1.\end{cases}$$

In this paper, we are interested in the average treatment effect for the whole population (ATE)

$$\tau=E[Y(1)-Y(0)],$$

and the average treatment effect for the treated (ATT)

$$\tau_{t}=E[Y(1)-Y(0)\mid W=1].$$

The assumptions listed below are essential for establishing limit properties and will remain part of the underlying assumptions throughout this article. These assumptions are commonly found in the existing literature for causal inference (see, for example, Rosenbaum and Rubin,, 1983; Heckman et al., 1998b, ; Abadie and Imbens,, 2016) and are considered standard assumptions.

Note that Assumption 1 is also known as “strong ignorability” (Rosenbaum and Rubin,, 1983), which implies that $X$ is sufficient for eliminating the impact of confounding factors (i.e., “unconfoundedness”). Hahn, (1998) established bounds for the asymptotic variance and explored asymptotically efficient estimation methods under Assumption 1(i).

Here, we define the proposed kernel-matching estimators for the ATE and ATT. Following the notation defined above, for $N\in\mathbb{N}$, suppose $\{Y_{i},W_{i},X_{i}\}_{i=1}^{N}$ are independent draws from the population with the distribution of $\{Y,W,X\}$. Let $N_{1}$ and $N_{0}$ be the number of treated and control units, respectively. The individual treatment effect is $\tau_{i}=Y_{i}(1)-Y_{i}(0)$. For every unit in the study, we can only get one of the potential outcomes, $Y_{i}(0)$ and $Y_{i}(1)$, and the other is unobservable or missing.

Let $K(\cdot)$ be a kernel function, and $h$ be the bandwidth parameter that converges to zero as $n$ goes to infinity. The kernel-matching estimators based on the propensity score $p(X)$ for the missing potential outcomes can be expressed as

$$\widehat{Y}_{i}(0)=\begin{cases}Y_{i}&\text{ if }W_{i}=0,\\ \frac{\sum_{j:W_{j}=0}^{N}Y_{j}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}{\sum_{j:W_{j}=0}^{N}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}&\text{ if }W_{i}=1,\end{cases}$$

and

$$\widehat{Y}_{i}(1)=\left\{\begin{array}[]{lll}\frac{\sum_{j:W_{j}=1}^{N}Y_{j}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}{\sum_{j:W_{j}=1}^{N}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}&\text{ if }&W_{i}=0,\\ Y_{i}&\text{ if }&W_{i}=1.\end{array}\right.$$

The proposed kernel-matching estimator for the ATE is defined as

$$\begin{split}\hat{\tau}_{N}=&\frac{1}{N}\sum_{i=1}^{N}\left(\hat{Y}_{i}(1)-\hat{Y}_{i}(0)\right)\\ =&\frac{1}{N}\sum_{i=1}^{N}\left(2W_{i}-1\right)\left(Y_{i}-\frac{\sum_{j:W_{j}=1-W_{i}}^{N}Y_{j}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}{\sum_{j:W_{j}=1-W_{i}}^{N}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}\right),\end{split}$$

and the proposed kernel-matching estimator for the ATT is defined as

$$\begin{split}\hat{\tau}_{t}=\frac{1}{N_{1}}\sum_{i=1}^{N}W_{i}\left(Y_{i}-\frac{\sum_{j=1}^{N}(1-W_{j})Y_{j}K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}{\sum_{j=1}^{N}(1-W_{j})K\left(\frac{p(X_{j})-p(X_{i})}{h}\right)}\right),\end{split}$$

where $N_{1}=\sum_{i=1}^{N}W_{i}$ is the number of treatment cases in the sample.

In observational studies, to estimate the propensity scores, following Rosenbaum and Rubin, (1983), a generalized linear model $p(X)=F(X^{\top}\beta)$ can be used. The logit and probit functions are commonly used parametric link functions $F(\cdot)$.

Let $\beta^{*}$ be the true value of the propensity score model parameter vector, such that $p(X)=F(X^{\top}\beta^{*})$. The estimator based on matching the true propensity score is given by $\widehat{\tau}_{N}(\beta^{*})$, and $\widehat{\tau}_{N}(\hat{\beta})$ is a kernel matching estimator based on the estimated propensity score $F(X^{\top}\hat{\beta})$, where $\hat{\beta}$ is the maximum likelihood estimator of $\beta$ based on the observed data $\{Y_{i},W_{i},X_{i}\}^{N}_{i=1}$, i.e.,

$$\hat{\beta}=\arg\max_{\beta}\ln L\left(\beta\mid W_{1},X_{1},\ldots,W_{N},X_{N}\right)$$

with the log-likelihood function

		$\displaystyle\ln L\left(\beta\mid W_{1},X_{1},\ldots,W_{N},X_{N}\right)$
		$\displaystyle\quad=\sum_{i=1}^{N}W_{i}\ln F\left(X_{i}^{\prime}\beta\right)+\left(1-W_{i}\right)\ln\left(1-F\left(X_{i}^{\prime}\beta\right)\right).$

Next, we define a type of matching estimator where the unknown parameter $\beta$ is estimated using the maximum likelihood method. The kernel-matching estimators based on the $F(X^{\top}\hat{\beta})$ for the ATE can be expressed as

$$\begin{split}\hat{\tau}_{N}(\hat{\beta})&=\frac{1}{N}\sum_{i=1}^{N}(\hat{Y}_{i}(1)-\hat{Y}_{i}(0))\\ &=\frac{1}{N}\sum_{i=1}^{N}\left(2W_{i}-1\right)\left(Y_{i}-\frac{\sum_{j:W_{j}=1-W_{i}}Y_{j}K\left(\frac{F(X_{j}^{\prime}\hat{\beta})-F(X_{i}^{\prime}\hat{\beta})}{h}\right)}{\sum_{j:W_{j}=1-W_{i}}K\left(\frac{F(X_{j}^{\prime}\hat{\beta})-F(X_{i}^{\prime}\hat{\beta})}{h}\right)}\right),\end{split}$$

and for ATT, the estimator can be expressed as

$$\begin{split}\hat{\tau}_{t,N}(\hat{\beta})=\frac{1}{N_{1}}\sum_{i=1}^{N}W_{i}\left(Y_{i}-\frac{\sum_{j=1}(1-W_{j})Y_{j}K\left(\frac{F(X_{j}^{\prime}\hat{\beta})-F(X_{i}^{\prime}\hat{\beta})}{h}\right)}{\sum_{j=1}(1-W_{j})K\left(\frac{F(X_{j}^{\prime}\hat{\beta})-F(X_{i}^{\prime}\hat{\beta})}{h}\right)}\right).\end{split}$$

In this section, we establish that the proposed kernel-matching estimators for the ATE and ATT are consistent and asymptotically normally distributed, with a convergence rate of $N^{1/2}$.

Let $g_{w}(p)=\mathbb{E}(Y(w)|p(X)=p)$, under Assumptions 1 and 2, we provide the following theorem to show that the estimators $\hat{\tau}$ and $\hat{\tau}_{t}$ are consistent estimators of $\tau$ and $\tau_{t}$, respectively, under mild conditions.

The proof of Theorem 1 is given in the Appendix.

Next, we elucidate the formal result for the asymptotic normality of the proposed estimators.

The proof of Theorem 2 is given in the Appendix.

In the following theorems, we will provide the large sample properties of $\hat{\tau}_{N}(\hat{\beta})$ and $\hat{\tau}_{t,N}(\hat{\beta})$, including the consistency in Theorem 3 and the asymptotic normality in Theorem 4.

The proofs of Theorem 3 and Theorem 4 are presented in the Appendix.

In this section, a Monte Carlo simulation study is used to evaluate the performance of the proposed estimators and compare the performance with those of existing estimators. The simulation results confirm the effectiveness of the kernel matching estimators based on the estimated propensity score in estimating ATE and ATT. For comparative purposes, the following existing ATE and ATT estimation methods are considered:

• •

  Adabie-Imbens matching estimators based on the original covariates (Abadie and Imbens,, 2006) (denoted as “Covariate”).

• •

  Adabie-Imbens matching estimators based on true propensity score (Abadie and Imbens,, 2006) (denoted as “True PS”).

• •

  Adabie-Imbens matching estimators based on estimated propensity score (Abadie and Imbens,, 2006, 2016) (denoted as “Estimated PS”).

• •

  Normalized inverse probability weighting (IPW) (Hirano et al.,, 2003; Imbens,, 2004) estimators for ATE and ATT given by

  $$\hat{\tau}_{\mathrm{IPW}}^{\mathrm{ATE}}=\sum_{i=1}^{N}\frac{W_{i}\cdot Y_{i}}{\hat{p}\left(X_{i}\right)}\left(\sum_{i=1}^{N}\frac{W_{i}}{\hat{p}\left(X_{i}\right)}\right)^{-1}-\sum_{i=1}^{N}\frac{\left(1-W_{i}\right)\cdot Y_{i}}{1-\hat{p}\left(X_{i}\right)}\left(\sum_{i=1}^{N}\frac{\left(1-W_{i}\right)}{1-\hat{p}\left(X_{i}\right)}\right)^{-1}$$

  and

  $$\hat{\tau}_{\mathrm{IPW}}^{\mathrm{ATT}}=\frac{1}{N_{1}}\sum_{i:W_{i}=1}Y_{i}-\sum_{i:W_{i}=0}Y_{i}\cdot\frac{\hat{p}(X_{i})}{1-\hat{p}(X_{i})}\left(\sum_{i:W_{i}=0}\frac{\hat{p}(X_{i})}{1-\hat{p}(X_{i})}\right)^{-1},$$

  respectively (denoted as “IPW”).

• •

  Doubly robust estimators (Robins et al.,, 1994; Bang and Robins,, 2005) for ATE and ATT given by

  	$\displaystyle\hat{\tau}_{\mathrm{DR}}^{\mathrm{ATE}}$	$\displaystyle=$	$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left\{\left[\frac{W_{i}Y_{i}}{\hat{p}\left(X_{i}\right)}-\frac{W_{i}-\hat{p}\left(X_{i}\right)}{\hat{p}\left(X_{i}\right)}\hat{\mu}\left(1,X_{i}\right)\right]\right.$
  			$\displaystyle\left.\qquad-\left[\frac{\left(1-W_{i}\right)Y_{i}}{1-\hat{p}\left(X_{i}\right)}-\frac{W_{i}-\hat{p}\left(X_{i}\right)}{1-\hat{p}\left(X_{i}\right)}\hat{\mu}\left(0,X_{i}\right)\right]\right\}$

  and

  $\displaystyle\hat{\tau}_{\mathrm{DR}}^{\mathrm{ATT}}=\frac{1}{N_{1}}\sum_{i=1}^{N}\left[W_{i}Y_{i}-\frac{\left(1-W_{i}\right)Y_{i}\hat{p}\left(X_{i}\right)+\hat{\mu}\left(0,X_{i}\right)\left(W_{i}-\hat{p}\left(X_{i}\right)\right)}{1-\hat{p}\left(X_{i}\right)}\right],$

  respectively (denoted as “DR”).

• •

  The kernel-matching estimators for ATE and ATT presented in Section 2 using the Gaussian kernel function with a bandwidth $h$ of 0.01 (denoted as “Proposed”). In preliminary simulations (results not shown), we also tested the Epanechnikov kernel function and various bandwidth values, finding that the Gaussian kernel with $h=0.01$ generally provided better performance for the kernel-matching estimator.

To simulate the control and treatment cases with covariates, the following five settings are considered:

• •

  Setting I:

  • –

    Two independent covariates, $X_{1}$ and $X_{2}$, follow the uniform distribution in the interval $[-1/2,1/2]$.

  • –

    The potential outcomes are generated based on the following models

    $$Y(0)=3X_{1}-3X_{2}+\epsilon_{0}$$

    and

    $$Y(1)=5+5X_{1}+X_{2}+\epsilon_{1},$$

    where $\epsilon_{0}$ and $\epsilon_{1}$ are random errors that follow the standard normal distribution independent of $(W,X_{1},X_{2})$.

  • –

    The treatment variable, $W$, is related to $(X_{1},X_{2})$ through the propensity score with a logistic function

    $$\Pr(W=1|X_{1}=x_{1},X_{2}=x_{2})=\frac{\exp(x_{1}+2x_{2})}{1+\exp(x_{1}+2x_{2})}.$$

  • –

    The ATE and ATT estimators use a single match ($M=1$).

• •

  Setting II:

  • –

    Two independent covariates, $X_{1}$ and $X_{2}$, follow the uniform distribution in the interval $[-1/2,1/2]$.

  • –

    The potential outcomes are generated based on the following models

    $$Y(0)=10X_{1}+\epsilon_{0}$$

    and

    $$Y(1)=5-10X_{1}+\epsilon_{1},$$

    where $\epsilon_{0}$ and $\epsilon_{1}$ are random errors that follow the standard normal distribution independent of $(W,X_{1},X_{2})$.

  • –

    The treatment variable, $W$, is related to $X_{2}$ through the propensity score with a logistic function

    $$\Pr(W=1|X_{1}=x_{1},X_{2}=x_{2})=\frac{\exp(2x_{2})}{1+\exp(2x_{2})}.$$

  • –

    The ATE and ATT estimators use a single match ($M=1$).

• •

  Setting III:

  • –

    Four covariates, the vector $X=(X_{1},X_{2},X_{3},X_{4})$ follows a four-dimensional normal distribution $\mathcal{N}\left(0,\mbox{\boldmath$\Sigma$}\right)$, where the variance-covariance $\Sigma$ is given by

    $\displaystyle\mbox{\boldmath$\Sigma$}=\frac{1}{3}\begin{pmatrix}1&-1&0&0\\ -1&2&0&0\\ 0&0&1&-1\\ 0&0&-1&2\\ \end{pmatrix},$

    i.e., $X_{1}$ and $X_{2}$ are correlated, $X_{3}$ and $X_{4}$ are correlated, and $(X_{1},X_{2})$ and $(X_{3},X_{4})$ are independent.

  • –

    The potential outcomes are generated based on the following models

    $$Y=0.4+0.25\sin(8(X^{\top}\mbox{\boldmath$\theta$})-5)+0.4\exp(-16(4(X^{\top}\mbox{\boldmath$\theta$})-2.5)^{2})+W+\epsilon,$$

    where $\mbox{\boldmath$\theta$}=(1,1,1,1)^{\top}$, $\epsilon$ is a random error that follows the standard normal distribution independent of $X$.

  • –

    The treatment variable, $W$, is related to $X=(X_{1},X_{2},X_{3},X_{4})$ through the propensity score with a logistic function

    	$\displaystyle\Pr(W=1$	$\displaystyle|$	$\displaystyle X_{1}=x_{1},X_{2}=x_{2},X_{3}=x_{3},X_{4}=x_{4})$
    		$\displaystyle=$	$\displaystyle\frac{\exp(x_{1}-3x_{2}+2x_{3}+x_{4})}{1+\exp(x_{1}-3x_{2}+2x_{3}+x_{4})}.$

  • –

    The ATE and ATT estimators use a single match ($M=1$).

• •

  Setting IV:

  • –

    Four covariates, the vector $X=(X_{1},X_{2},X_{3},X_{4})$ follows a four-dimensional normal distribution $\mathcal{N}\left(0,\mbox{\boldmath$\Sigma$}\right)$, where the variance-covariance $\Sigma$ is given by

    $\displaystyle\mbox{\boldmath$\Sigma$}=\frac{1}{3}\begin{pmatrix}1&-1&0&0\\ -1&2&0&0\\ 0&0&1&-1\\ 0&0&-1&2\\ \end{pmatrix}.$

  • –

    The potential outcomes are generated based on the following models

    $$Y=0.15+0.7\left[\frac{\exp(\sqrt{2}(X^{\top}\mbox{\boldmath$\theta$}))}{1+\exp(\sqrt{2}(X^{\top}\mbox{\boldmath$\theta$}))}\right]+W+\epsilon,$$

    where $\epsilon$ is a random error that follows the standard normal distributions independent of $X$.

  • –

    The treatment variable, $W$, is related to $X=(X_{1},X_{2},X_{3},X_{4})$ through the propensity score with a logistic function

    	$\displaystyle\Pr(W=1$	$\displaystyle|$	$\displaystyle X_{1}=x_{1},X_{2}=x_{2},X_{3}=x_{3},X_{4}=x_{4})$
    		$\displaystyle=$	$\displaystyle\frac{\exp(x_{1}+x_{2}+x_{3}+x_{4})}{1+\exp(x_{1}+x_{2}+x_{3}+x_{4})}.$

  • –

    The ATE and ATT estimators use a single match ($M=1$).

• •

  Setting V:

  • –

    Three independent covariates $X_{1}\sim\mathcal{N}(0,1)$, $X_{2}\sim\text{Bernoulli}(0.5)$, and $X_{3}\sim\mathcal{N}(0,1)$.

  • –

    The potential outcomes are generated based on the following models

    $$Y(0)=\sin(1+X_{1}+0.5X_{2})+\epsilon_{0}$$

    and

    $$Y(1)=\sin(2+X_{1}+0.5X_{2})+\epsilon_{1},$$

    where $\epsilon_{0}$ and $\epsilon_{1}$ are random errors that follow the standard normal distributions independent of $X$.

  • –

    The treatment variable, $W$, is related to $X=(X_{1},X_{2},X_{3})$ through the propensity score with a logistic function

    	$\displaystyle\Pr(W=1$	$\displaystyle|$	$\displaystyle X_{1}=x_{1},X_{2}=x_{2},X_{3}=x_{3})$
    		$\displaystyle=$	$\displaystyle\frac{\exp(2x_{1}+x_{2}-3x_{3}+3)}{1+\exp(2x_{1}+x_{2}-3x_{3}+3)}.$

  • –

    The ATE and ATT estimators use a single match ($M=1$).

We generated datasets with sample sizes of $N=200,500$, and $1000$, and conducted $1000$ Monte Carlo simulations for each scenario to estimate causal effects. The data generation processes for Settings I and II follow the approach of Abadie and Imbens, (2016), while those for Settings III and IV are based on the method of Frölich, (2004) with slight modifications. Similarly, the process for Setting V is adapted from the approach of Luo and Zhu, (2020), also with slight modifications.

We compare the point estimation of ATE and ATT of the six estimators (Covariate, True PS, Estimated PS, Proposed, IPW, and DR) in terms of the simulated biases (Bias), standard deviations (SD), and root mean square errors (RMSE). For comparative purposes, we also compute the 95% confidence intervals based on different estimators using bootstrap methods. Specifically, following Abadie and Imbens, (2006, 2016), we obtain the 95% confidence intervals based on the “Covariate”, “True PS” and “Estimated PS” using normal approximation with standard errors obtained by using 400 bootstrap samples. For the ” “Proposed”, “IPW” and “DR” estmators, the 95% confidence intervals are obtained using the percentile bootstrap method (Efron and Tibshirani,, 1993) with 400 bootstrap samples. We compare the interval estimation of ATE and ATT of the six estimators in terms of the simulated average widths (AW) and the coverage probabilities (CP). Simulation results are reported for the six estimators under Settings I–V in Tables 1–5 for ATE estimation and in Tables 6–10 for ATT estimation.

From Tables 1–10, we observe that the biases and RMSEs of the proposed estimator decrease as the sample size increases, which is consistent with the theoretical asymptotic results discussed in Section 3. For point estimation of both ATE and ATT, the proposed kernel matching estimator demonstrates comparable bias and RMSE performance across most settings, though it does not achieve the smallest bias and MSE in every scenario. From Tables 1–2, we observe that the proposed estimator for ATE outperforms the matching estimators based on covariates, true propensity score, and estimated propensity score in Settings I and II in terms of RMSE. From Tables 3–4, we observe that the proposed estimator for ATE shows a lower bias than the covariate-based matching estimator in Settings III and IV. This consistent performance highlights the effectiveness of the proposed estimator across varied settings, even when not uniformly outperforming all competitors.

In terms of interval estimation, the proposed estimator achieves coverage probabilities close to the nominal 95% level, suggesting strong reliability of the proposed method. While IPW estimator offers good point estimation, it struggles with interval estimation accuracy, failing to maintain 95% coverage in Settings IV and V. Similarly, the DR estimator performs well for point estimation in Settings I–IV but falters in Setting V. The proposed estimator’s interval estimation remains stable, often surpassing IPW and DR in coverage accuracy, particularly in settings with extreme true propensity scores. Based on these simulation results, DR, IPW, and the proposed estimators are generally recommended; however, the proposed estimator is especially favorable for both point and interval estimation in complex scenarios like Settings IV and V.

In this section, we apply the estimation procedures to analyze data from the National Supported Work (NSW) Demonstration, a program funded by both federal and private sources in the mid-1970s. This initiative provided 6–18 months of work experience for individuals facing economic and social challenges before enrollment. Our analysis demonstrates the impact of the NSW Demonstration labor training program on the earnings of these individuals. To evaluate the effectiveness of the proposed estimator, we estimate the average treatment effect on the treated (ATT) using three control groups. We use two datasets: one drawn from the experimental sample in LaLonde, (1986) and the other from Westat’s Matched Current Population Survey–Social Security Administration File (CPS-3) Dehejia and Wahba, (1999) ¹¹1The data can be obtained from the website at https://users.nber.org/$\sim$rdehejia/nswdata2.html..

In this dataset, the outcome variable (Re78) represents each individual’s real earnings in 1978, and the treatment variable (Treat) indicates whether the individual enrolled in the labor training program. The ten potential confounding variables include age (Age), years of schooling (Educ), race indicators for Black (Black) and Hispanic (Hispanic), marital status (Married), high school diploma status (Nodegr), real earnings in 1974 (Re74) and 1975 (Re75), and indicators for zero earnings in 1974 (U74) and 1975 (U75). This subset comprises 185 individuals in the treated group, 260 individuals in the control group reported by LaLonde, (1986), and 429 individuals from the CPS-3 control group. Table 11 presents the summary statistics for the NSW dataset and the $p$-values of the $t$-tests for comparing the means between the treated group and the two control groups from the two samples. Table 11 shows that the treated and control groups are well-balanced in terms of mean values for most variables, with the exception of the variable Nodegr. For the nine remaining covariates, we cannot reject the null hypothesis that the means are equal, indicating similar averages across treated and control groups. However, when comparing treated units from the experimental data with control units from the non-experimental CPS-3 data, significant mean differences appear in eight out of the ten covariates. The only covariates for which the null hypothesis of equal means cannot be rejected are Hispanic and education.