Matching Using Sufficient Dimension Reduction for Causal Inference

TL;DR

This paper introduces a novel matching method for causal inference using sufficient dimension reduction to create effective covariates, improving estimation accuracy under weaker assumptions.

Contribution

It proposes a nonparametric, sufficient dimension reduction-based matching approach that enhances causal effect estimation with weaker support conditions and flexible matching.

Findings

Method performs well in simulations

Real data analysis confirms effectiveness

Reduces bias compared to traditional methods

Abstract

To estimate casual treatment effects, we propose a new matching approach based on the reduced covariates obtained from sufficient dimension reduction. Compared to the original covariates and the propensity score, which are commonly used for matching in the literature, the reduced covariates are estimable nonparametrically under a mild assumption on the original covariates, and are sufficient and effective in imputing the missing potential outcomes. Under the ignorability assumption, the consistency of the proposed approach requires a weaker common support condition. In addition, the researchers are allowed to use different reduced covariates to find matched subjects for different treatment groups. We develop relative asymptotic results, and conduct simulation studies as well as real data analysis to illustrate the usefulness of the proposed approach.

Figures

Tables

Table 1. Table 1: Simulation Results for A ​ C ​ E 𝐴 𝐶 𝐸 ACE in Case 1

	I	I​I	I​I​I	I​V
	BIAS
Ambient	-0.1460	0.0929	0.4764	0.0571
Estimated PS	-0.0107	0.0317	0.0415	0.0154
True PS	0.0102	0.0073	0.0211	0.0013
Genetic Matching	-0.1572	0.0228	0.3116	0.0222
Active Set (Oracle)	-0.0138	0.0172	0.2005	0.0571
SDR (Oracle)	-0.0138	0.0006	0.0164	0.0014
Proposed	-0.0065	0.0029	0.0220	0.0016
	SD
Ambient	0.4258	0.0849	0.1388	0.0549
Estimated PS	0.6463	0.1051	0.1751	0.0567
True PS	1.1884	0.1856	0.2119	0.0537
Genetic Matching	0.3308	0.0637	0.1283	0.0554
Active Set (Oracle)	0.1583	0.0559	0.1054	0.0548
SDR (Oracle)	0.1583	0.0535	0.0917	0.0537
Proposed	0.1788	0.0551	0.1011	0.0525
	RMSE
Ambient	0.4499	0.1258	0.4962	0.0792
Estimated PS	0.6460	0.1097	0.1798	0.0588
True PS	1.1878	0.1856	0.2129	0.0537
Genetic Matching	0.3661	0.0677	0.3370	0.0597
Active Set (Oracle)	0.1588	0.0585	0.2265	0.0792
SDR (Oracle)	0.1588	0.0535	0.0931	0.0537
Proposed	0.1789	0.0552	0.1034	0.0525

Table 2. Table 2: Simulation Results for A ​ C ​ E 𝐴 𝐶 𝐸 ACE and A ​ C ​ E ​ T 𝐴 𝐶 𝐸 𝑇 ACET in Case 2

	$ACE$		$ACET$
	I${}^{\ast }$	I​I${}^{\ast }$	I${}^{\ast }$	I​I${}^{\ast }$
	BIAS
Ambient	-0.6805	0.0279	0.3132	0.0293
Estimated PS	-0.4273	-0.0609	-0.3044	0.0059
True PS	-0.0209	-0.0072	0.0166	0.0047
Genetic Matching	-0.6165	0.0125	0.0765	0.0148
Active Set (Oracle)	-0.0961	0.0036	0.0113	0.0056
SDR (Oracle)	-0.0961	-0.0003	0.0113	0.0002
Proposed	-0.0761	0.0021	0.0187	0.0032
	SD
Ambient	0.5151	0.0957	0.5796	0.1109
Estimated PS	0.8401	0.1312	0.7713	0.1419
True PS	2.1383	0.2997	1.1095	0.2066
Genetic Matching	0.4368	0.0763	0.4526	0.0852
Active Set (Oracle)	0.1822	0.0645	0.2712	0.0732
SDR (Oracle)	0.1822	0.0615	0.2712	0.0669
Proposed	0.2202	0.0653	0.2896	0.0729
	RMSE
Ambient	0.8533	0.0996	0.6586	0.1146
Estimated PS	0.9421	0.1446	0.8288	0.1419
True PS	2.1373	0.2996	1.1091	0.2066
Genetic Matching	0.6572	0.0623	0.3596	0.0681
Active Set (Oracle)	0.2059	0.0645	0.2713	0.0734
SDR (Oracle)	0.2059	0.0615	0.2713	0.0669
Proposed	0.2329	0.0653	0.2901	0.0729

Table 3. Table 3: Simulation Results for A ​ C ​ E 𝐴 𝐶 𝐸 ACE in Case 3

	A	B	C	D	E	F	G
	BIAS
Ambient	0.0551	0.0527	0.0508	0.0625	0.0586	0.0654	0.0570
Estimated PS	0.0075	0.0041	0.0021	0.0073	0.0124	0.0097	0.0149
True PS	0.0038	0.0020	0.0078	0.0046	0.0064	0.0041	0.0070
Genetic Matching	0.0343	0.0336	0.0246	0.0371	0.0353	0.0371	0.0321
Active Set (Oracle)	0.0247	0.0240	0.0208	0.0295	0.0287	0.0328	0.0247
SDR (Oracle)	0.0023	0.0019	0.0025	0.0027	0.0025	0.0025	0.0023
Proposed	0.0035	0.0037	0.0043	0.0050	0.0049	0.0055	0.0034
	SD
Ambient	0.0420	0.0431	0.0427	0.0445	0.0447	0.0440	0.0450
Estimated PS	0.0692	0.0721	0.0661	0.0757	0.0745	0.0738	0.0731
True PS	0.0943	0.0950	0.1093	0.0983	0.1015	0.1004	0.1177
Genetic Matching	0.0596	0.0590	0.0558	0.0600	0.0619	0.0574	0.0521
Active Set (Oracle)	0.0587	0.0593	0.0587	0.0601	0.0605	0.0625	0.0636
SDR (Oracle)	0.0327	0.0335	0.0329	0.0331	0.0337	0.0334	0.0339
Proposed	0.0365	0.0355	0.0356	0.0374	0.0372	0.0369	0.0357
	RMSE
Ambient	0.0693	0.0681	0.0663	0.0767	0.0737	0.0788	0.0726
Estimated PS	0.0695	0.0722	0.0661	0.0761	0.0755	0.0744	0.0745
True PS	0.0943	0.0950	0.1096	0.0984	0.1016	0.1004	0.1178
Genetic Matching	0.0687	0.0679	0.0609	0.0705	0.0712	0.0683	0.0612
Active Set (Oracle)	0.0637	0.0640	0.0623	0.0669	0.0670	0.0705	0.0682
SDR (Oracle)	0.0328	0.0336	0.0329	0.0331	0.0338	0.0334	0.0340
Proposed	0.0366	0.0356	0.0359	0.0377	0.0375	0.0373	0.0358

Equations

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

$$ACET=E[Y(1)-Y(0)|T=1].$$

$$Y(t)\perp \perp T|X,t=0,1$$

$$Y(t)\perp \perp T|R(X),$$

$${\hat{Y}}^{\text{mbox}(i)}(1-T^{\text{mbox}(i)})=\frac{1}{m}\sum _{j\in J_m(i)}{Y}^{\text{mbox}(j)}(T^{\text{mbox}(j)}),$$

$$W\perp \perp X|X^{\prime }\beta .$$

$$E[\{X-E(X)\}|X^{\prime }\beta ]={\Sigma }_X\beta ({\beta }^{\prime }{\Sigma }_X\beta {)}^{-1}{\beta }^{\prime }\{X-E(X)\},$$

$$Y(t)\perp \perp T|X^{\prime }{\beta }_t=a.$$

$$Y(t)\perp \perp X|X^{\prime }{\beta }_t,$$

$$f(y(t)|X^{\prime }{\beta }_t=a,T=t)$$

$$Y(0)\perp \perp X|(X^{\prime }{\beta }_0^{\text{mbox}D},T=0),Y(1)\perp \perp X|(X^{\prime }{\beta }_1^{\text{mbox}D},T=1),$$

$$Y(0)\perp \perp X|(X^{\prime }{\beta }_0,T=0).$$

$$f(\cdot |X=x,T=0)$$

$$Y(0)\perp \perp X|X^{\prime }{\beta }_0^{\text{mbox}D},$$

$$f(\cdot |X=x)=f(\cdot |X^{\prime }{\beta }_0=x^{\prime }{\beta }_0)=f(\cdot |X^{\prime }{\beta }_0=x^{\text{mbox}\ast \prime }{\beta }_0)=f(\cdot |X=x^{\text{mbox}\ast }).$$

$$E({\hat{\mu }}_X)-\mu =O(n^{\text{mbox}-2/p}),\text{mbox}and\mathrm{var}({\hat{\mu }}_X)=O_P(n^{\text{mbox}-\min \{1,6/p\}}),$$

$$E({\hat{\mu }}_{\pi })-\mu =O(n^{\text{mbox}-2}),\text{mbox}and\mathrm{var}({\hat{\mu }}_{\pi })=O_P(n^{\text{mbox}-1}).$$

$$E({\hat{\mu }}_r)-\mu =O(n^{\text{mbox}-2/\max \{r(0),r(1)\}}),\text{mbox}and\mathrm{var}({\hat{\mu }}_r)=O_P(n^{\text{mbox}-\min \{1,6/r(0),6/r(1)\}}).$$

$$n\left[\{1-(r+s)n^{\text{mbox}-1}{\}}^{\text{mbox}n}-(1-r{n}^{\text{mbox}-1}{)}^{\text{mbox}n}(1-s{n}^{\text{mbox}-1}{)}^{\text{mbox}n}\right]\to -rs\exp \{-(r+s)\},$$

$$n[\{1-(r+s)n^{\text{mbox}-1}{\}}^{\text{mbox}n}-\{1-(r+s)n^{\text{mbox}-1}+rs{n}^{\text{mbox}-2}{\}}^{\text{mbox}n}],$$

$$n\left[1-[1+rs{n}^{\text{mbox}-2}\{1-(r+s)n^{\text{mbox}-1}{\}}^{-1}{]}^{\text{mbox}n}\right]\to -rs.$$

$$\underset{n\to \infty }{\lim }(1-{ϵ}_n/n{)}^{\text{mbox}n}=e^{\text{mbox}rs},$$

$${\hat{\mu }}_X-\mu$$

$$E\{Y(1-T)|X^{\text{mbox}(i)}\}-E\{Y(1-T)|X^{\text{mbox}(J(i))}\}=(X^{\text{mbox}(i)}-X^{\text{mbox}(J(i))}{)}^{\prime }{g}_i^{\text{mbox}(i)}$$

$$+(X^{\text{mbox}(i)}-X^{\text{mbox}(J(i))}{)}^{\prime }{g}_2^{\text{mbox}(i)}(X^{\text{mbox}(i)}-X^{\text{mbox}(J(i))})+g_3^{\text{mbox}(i,J(i))}\Vert X^{\text{mbox}(i)}-X^{\text{mbox}(J(i))}{\Vert }^{\text{mbox}3},$$

$$X^{\text{mbox}(J(i))}-X^{\text{mbox}(i)}=O_P(n^{\text{mbox}-1/p}),E(X^{\text{mbox}(J(i))}-X^{\text{mbox}(i)})=O(n^{\text{mbox}-2/p}),$$

$$\mathrm{cov}\{(X^{\text{mbox}(J(i))}-X^{\text{mbox}(i)}{)}^{\text{mbox}\otimes k},(X^{\text{mbox}(J(j))}-X^{\text{mbox}(j)}{)}^{\text{mbox}\otimes l}\}=O(n^{\text{mbox}-1-2/p}).$$

$$f(u)$$

$$f(u,v)$$

$$P(\Vert X-a\Vert \le \Vert u\Vert n_{1-t}^{\text{mbox}-1/p})+P(\Vert X-b\Vert \le \Vert v\Vert n_{1-t}^{\text{mbox}-1/p})$$

$$-P(\Vert X-a\Vert \le \Vert u\Vert n_{1-t}^{\text{mbox}-1/p}\text{mbox}and\Vert X-b\Vert \le \Vert v\Vert n_{1-t}^{\text{mbox}-1/p})$$

$$\equiv I+II+III$$

$$f(u,v)=(1-n_{1-t}^{\text{mbox}-1})f(a+u{n}_{1-t}^{\text{mbox}-1/p})f(b+v{n}_{1-t}^{\text{mbox}-1/p})(1-r_n{n}_{1-t}^{-1}-s_n{n}_{1-t}^{-1}{)}^{\text{mbox}{n}_{1-t}-2}$$

Full text

Key words and phrases: Causal inference; Central subspace; Common support condition; Dimension reduction; Matching

1 Introduction

One of the commonly used approaches to the causal analysis of observational data is matching, which is a systematic way to find comparable treated and control subjects that have identical or approximate values on an appropriate function of the covariates (Rosenbaum and Rubin, 1983). The causal effects are then estimated from the matched dataset as if the treatment were randomly assigned. Given a function of the covariates, multiple methods have been developed in the literature to assign matched subjects, e.g., nearest available matching based on propensity scores, Mahalanobis metric matching (Rosenbaum and Rubin, 1985), full matching (Rosenbaum, 1991; Hansen, 2004), and genetic matching (Diamond and Sekhon, 2013). The latter two are data-adaptive in the sense that they both optimize a certain balance criterion based on the data.

Besides the choice of the matching methods, it is also crucial to select an appropriate function of the covariates that these methods are built on, as different choices can dramatically affect the asymptotic behaviors of the resulting estimators (Heckman et al., 1997; Abadie and Imbens, 2006). A natural choice of such functions is the original covariates, which has been commonly used in the literature (Rosenbaum and Rubin, 1983). In addition to its ease of use, it ensures the unbiasedness of matching in the population level, under the ignorability assumption and the strong common support condition reviewed in Section $2$ below. However, as pointed out in Abadie and Imbens (2006), this choice suffers from the “curse of dimensionality”, which makes matching ineffective when the dimension of the covariates is relatively large.

To address the above-mentioned issue, lower-dimensional functions of covariates have been used for matching in the literature, among which the most popular choice is the propensity score (e.g., Rosenbaum and Rubin, 1985), defined as the probability of a subject being treated given its covariates. Same as the original covariates, the propensity score also ensures the unbiasedness of matching under the ignorability assumption and the strong common support condition. However, the estimation of the propensity score commonly relies on parametric models, which is susceptible to model mis-specification. Recently, de Luna et al. (2011) considered the variable selection issue in causal inference, by assuming that a small subset of covariates, called the active set, is sufficient in predicting the outcomes. Using the active set for matching, the resulting estimator is asymptotically unbiased, and is effective in finite samples.

In all the existing matching procedures, the same function of covariates is used to find matches in different treatment groups. When the outcome regression pattern differs, which in de Luna et al. (2011), means varying active sets in different treatment groups, these sets will be merged into one set for subsequent matching. The merged active set will then always contain noisy information for both treatment groups. In this case, it is conceivable that the effectiveness of matching can be enhanced, if we allow for different functions of covariates to be used for different treatment groups.

In this article, we introduce a new set of functions of the covariates with the aid of sufficient dimension reduction (SDR, Li and Duan,1989; Li,1991). Unlike the existing approaches, the functions of covariates vary in different treatment groups to find matches. Under a mild assumption, the new functions enjoy the advantages of both the original covariates and the propensity scores. That is, they are easily estimable in a model-free manner, and are usually of low dimension in practice. In particular, their dimensions are always equal to or lower than the cardinality of the active set in de Luna et al. (2011). In addition, as seen later, the new functions help relax the strong common support condition, and thus enhance the applicability of the matching methods.

The structure of this article is as follows. In Section $2$, we will briefly review the potential outcomes framework and the literature of SDR, with an emphasis on sliced inverse regression (SIR, Li, 1991). We then propose a new matching estimator using SDR in Section $3$, and develop some relative asymptotic results in Section $4$. Section $5$ includes the details in implementation. Simulation studies and real data analysis are conducted in Section $6$ and $7$ to illustrate the proposed estimator. We will leave more discussion to Section $8$.

2 A review of literature

2.1 Potential Outcomes Framework

Matching can be explained under the commonly used potential outcomes framework (Splawa-Neyman et al., 1990; Rubin, 1974) that assumes the complete data to be $n$ copies of $(X,T,Y(0),Y(1))$, denoted as \{(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}},T^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}},Y^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(0), Y^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(1)),i=1,\ldots n\}. Here, X\in\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}} is the set of covariates with support $\Omega(X)$; $T$ is assumed to be the binary treatment assignment with $T=0$ being the control group and $T=1$ being the treatment group; and $Y(t)\in\mathbb{R}$ is the potential outcome under treatment $t$ and is observed only when $T=t$. The potential outcome framework assumes the existence of the hypothetical outcome under the treatment level that is not observed in the data. Different types of causal effects are defined depending on the research interest. For example, the average causal effect is defined as

[TABLE]

and the average causal effect among the treated is defined as

[TABLE]

To estimate these parameters from the observed data, researchers commonly regulate the missingness mechanism by the ignorability assumption (Rosenbaum and Rubin, 1983),

[TABLE]

in which $\perp\perp$ means independence. Under (1), Rosenbaum and Rubin (1983) introduced a family of measurable functions of $X$ called the balancing score, denoted by $R(X)$, which satisfies

[TABLE]

for $t=0,1$. Examples of $R(X)$ include the covariates $X$ and the aforementioned propensity score.

A natural interpretation of (2) is that if we find two subjects that share the same value of $R(X)$, then no matter whether they are from the same treatment group or not, their outcomes $Y(t)$ should have identical distributions. Consequently, if we manage to match each subject with those in the alternative treatment group according to their similarity in $R(X)$, and impute its missing potential outcome using the observed outcomes in the matched set, then we can treat the imputed dataset as observed completely, and estimate the aforementioned causal effects using the usual sample moments. This strategy is used in the Abadie-Imbens estimator (Abadie and Imbens, 2006), which imputes the missing potential outcome of subject $i$ by

[TABLE]

in which $J_{\scriptscriptstyle m}(i)$ is the matched set for subject $i$ with cardinality $m$. Clearly, $J_{\scriptscriptstyle m}(i)$ is determined by $R(X)$, and is crucial to the consistency of matching. For a fixed value of $m$, when a high-dimensional $R(X)$ is used, such as $X$, the subjects in $J_{\scriptscriptstyle m}(i)$ may differ dramatically in terms of $R(X)$. This may lead to a large bias in (3). When $R(X)$ fails to satisfy (2), which occurs if the propensity score is inconsistently estimated, then even if the subjects in $J_{\scriptscriptstyle m}(i)$ share the same value of $R(X)$ as subject $i$, their outcomes can be stochastically different and lead to a large bias again. These are the drawbacks of using the original covariates and the estimated propensity score for matching.

To ensure the consistency of matching using $R(X)$, the data must satisfy that $P(T=1|R(X))\in(0,1)$ almost surely, which, by Bayes’ Theorem, amounts to the coincidence of the support of $R(X)$ in the treatment and control groups. When $R(X)$ is either the covariates or the propensity score, the condition is further equivalent to the common support of the covariates across the treatment groups, which we refer to as the strong common support condition. The condition can be relaxed if other forms of $R(X)$ are used, for example, if we employ SDR to construct a low-dimensional balancing score $R(X)$, other than the propensity score, for matching. We review the basics of SDR in the following subsection.

2.2 Sufficient dimension reduction

SDR has been widely used to reduce the dimension of the covariates prior to data analysis. For the covariates $X$ and a response variable $W$, it assumes the existence of \beta\in\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p\times d}}}, where $d<p$, such that the reduced covariates $X^{\prime}\beta$ are sufficient for the subsequent modeling of $W$. That is,

[TABLE]

For identifiable parametrization, Cook (1998) introduced the central subspace, denoted by ${\mathcal{S}}_{\scriptscriptstyle W|X}$, as the column space spanned by $\beta$ in (4) with the smallest dimension $d$. Conditions required for existence of this space are fairly general, and are adopted throughout the article. If $\beta$ in (4) is further restricted so that each of its column must have exactly one nonzero component, then $X^{\prime}\beta$ becomes the active set in variable selection. Thus, by definition, the dimension of the central subspace must be equal to or less than the cardinality of the active set. In this sense, SDR can be more effective than variable selection in reducing the dimension of the covariates, subject to consistent estimation of the central subspace.

In the literature, multiple methods have been developed to estimate the central subspace, including sliced inverse regression (SIR, Li,1991), sliced average variance estimation (SAVE, Cook and Weisberg,1991) and the directional regression (Li and Wang, 2007), among which SIR is commonly used for its simplicity and finite sample effectiveness. For consistency, it requires the linearity condition on $X$:

[TABLE]

where $\Sigma_{\scriptscriptstyle X}$ is the covariance matrix of $X$. The condition includes the elliptical distribution of $X$ as a special case, and holds approximately for general high-dimensional covariates (Hall and Li, 1993). Thus, the linearity condition is not considered restrictive in applications. As no other model structure is assumed, SIR is termed model-free. More details can be found in Section $5$ and Li (1991).

A separate issue in SDR is to determine the dimension of the central subspace. In the literature, the sequential testing procedure (Li, 1991) has been shown effective for this purpose. For continuity of the context, we leave its detail to the Appendix A (Supplementary Material).

3 Matching based on SDR

From the discussions above, when the SDR assumption (4) holds in the data with the treated and control outcomes to be the response variable $W$, respectively, the reduced covariates $X^{\prime}\beta$ are low-dimensional, and are estimable nonparametrically had the outcomes been completely observed. Therefore, $X^{\prime}\beta$ will be desired functions of covariates to be matched on if they are balancing scores in (2). This motivation is justified in the following theorem.

Theorem 1**.**

For $t=0,1$, suppose that the central subspace ${\mathcal{S}}_{\scriptscriptstyle Y(t)|X}$ is $r(t)$-dimensional, where $r(t)<p$. Then under the ignorability assumption, its arbitrary basis matrix \beta_{\scriptscriptstyle t}\in\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p\times r(t)}}} generates a balancing score $X^{\prime}\beta_{\scriptscriptstyle t}$ that satisfies (2) for $t=0,1$.

Proof.

Without loss of generality, we denote $f(y(t)|X,T=t)$ as the density function of the conditional distribution $Y(t)|(X,T=t)$ with respect to a $\sigma$-finite measure $\nu$, for $t=0,1$. For any random element $R$, denote $\Omega(R|T=t)$ and $\Omega(R)$ as the support of $R|T=t$ and $R$, respectively. For any $a\in\Omega(X^{\prime}\beta_{\scriptscriptstyle t})$, if $a\not\in\Omega(X^{\prime}\beta_{\scriptscriptstyle t}|T=t)$, then $X^{\prime}\beta_{\scriptscriptstyle t}=a$ implies that $T=1-t$ almost surely, which automatically implies that

[TABLE]

Thus we only need to show (6) for any $a\in\Omega(X^{\prime}\beta_{\scriptscriptstyle t}|T=t)$. That is, $f(y(t)|X^{\prime}\beta_{\scriptscriptstyle t}=a,T=t)=f(y(t)|X^{\prime}\beta_{\scriptscriptstyle t}=a)$. By the definition of $\beta_{\scriptscriptstyle t}$, we have

[TABLE]

which means that $f(y(t)|X)=f(y(t)|X^{\prime}\beta_{\scriptscriptstyle t})$. For any $x\in\Omega(X|T=t)$ such that $X^{\prime}\beta_{\scriptscriptstyle t}=a$, $f(y(t)|X=x,T=t)=f(y(t)|X=x)=f(y(t)|X^{\prime}\beta_{\scriptscriptstyle t}=a)$, where the first equality is due to the ignorability assumption. Thus we have

[TABLE]

This completes the proof. ∎

Based on this theorem, if we can estimate the central subspaces ${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}$ and ${\mathcal{S}}_{\scriptscriptstyle Y(1)|X}$ consistently, say by ${\mathcal{S}}(\hat{\beta}_{\scriptscriptstyle 0})$ and ${\mathcal{S}}(\hat{\beta}_{\scriptscriptstyle 1})$, where ${\mathcal{S}}(\cdot)$ denotes the linear span of a matrix, then we can use the reduced covariates $X^{\prime}\hat{\beta}_{\scriptscriptstyle 0}$ in place of $X$ to find matches for each treated subject, and use $X^{\prime}\hat{\beta}_{\scriptscriptstyle 1}$ to find matches for each control subject.

Because we allow the two central subspaces to differ, we may end up using different criteria for different treatment groups in the Abadie-Imbens estimator. As mentioned before, this flexibility is meaningful, because it allows us to use the exactly useful information in the covariates without introducing noise.

Due to the missingness in the observed data, the central subspaces in Theorem 1 are not readily estimable by any aforementioned SDR method. The estimable central subspaces are those conditional on the treatment assignments, denoted by {\mathcal{S}}_{\scriptscriptstyle Y(0)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}} and {\mathcal{S}}_{\scriptscriptstyle Y(1)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}, where the superscript “$D$” refers to data. They are spanned by the basis matrices \beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}} and \beta_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}} that satisfy

[TABLE]

respectively. Fortunately, under the ignorability assumption and a mild additional condition, these estimable central subspaces are identical to their counterparts of interest, and thus can serve as the substitutes for the desired central subspaces. In the following, we denote $\Omega(R|T=t)$ and $\Omega(R)$ as the support of $R|T=t$ and the support of $R$, respectively, for any random element $R$ and $t=0,1$.

Theorem 2**.**

Under the ignorability assumption (1), for $t=0,1$, suppose that $\Omega(X^{\prime}\beta_{\scriptscriptstyle t}|T=t)=\Omega(X^{\prime}\beta_{\scriptscriptstyle t})$, then {\mathcal{S}}_{\scriptscriptstyle Y(t)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}={\mathcal{S}}_{\scriptscriptstyle Y(t)|X}.

Proof.

For simplicity, we only prove the theorem for $t=0$. The case for $t=1$ can be shown in the same manner. First, we show that {\mathcal{S}}_{\scriptscriptstyle Y(0)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}\subseteq{\mathcal{S}}_{\scriptscriptstyle Y(0)|X}, which is equivalent to:

[TABLE]

Following the notations in the proof of Theorem 1, we denote $f(\cdot|X)$ as the conditional density function of $Y(0)|X$ and $f(\cdot|X,T=0)$ as the conditional density function of $Y(0)|X$ when $T=0$. (7) is equivalent to that $f(\cdot|X=x,T=0)=f(\cdot|X^{\prime}\beta_{\scriptscriptstyle 0}=x^{\prime}\beta_{\scriptscriptstyle 0},T=0)$ for any $x\in\Omega(X|T=0)$. By the ignorability assumption, $f(\cdot|X=x,T=0)=f(\cdot|X=x)$, and by the definition of ${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}$, $f(\cdot|X=x)=f(\cdot|X^{\prime}\beta_{\scriptscriptstyle 0}=x^{\prime}\beta_{\scriptscriptstyle 0})$. Thus $f(\cdot|X=x,T=0)$ is measurable with respect to $x^{\prime}\beta_{\scriptscriptstyle 0}$, which means that

[TABLE]

Hence (7) holds. Conversely, to show that {\mathcal{S}}_{\scriptscriptstyle Y(0)|X}\subseteq{\mathcal{S}}_{\scriptscriptstyle Y(0)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}, note that it is equivalent to

[TABLE]

which holds if f(\cdot|X=x)=f(\cdot|X^{\prime}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}=x^{\prime}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}) for any $x\in\Omega(X)$. Because $\Omega(X^{\prime}\beta_{\scriptscriptstyle 0}|T=0)=\Omega(X^{\prime}\beta_{\scriptscriptstyle 0})$, there exists x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}\in\Omega(X|T=0) such that x^{\prime}\beta_{\scriptscriptstyle 0}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}. By the definition of ${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}$, we have

[TABLE]

Since x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}\in\Omega(X|T=0), by the ignorability assumption, f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}})=f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}},T=0), which, by the definition of {\mathcal{S}}_{\scriptscriptstyle Y(0)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}, further implies that f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}})=f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}},T=0). Thus f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}) is measurable with respect to x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}, which, similar to the above, implies that f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}})=f(\cdot|X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}). Because {\mathcal{S}}_{\scriptscriptstyle Y(0)|X}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}\subseteq{\mathcal{S}}_{\scriptscriptstyle Y(0)|X}, we have x^{\prime}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}=x^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*\prime}}}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}. Hence f(\cdot|X=x)=f(\cdot|X^{\prime}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}=x^{\prime}\beta_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle D}}}), which implies (8). This completes the proof. ∎

The condition in this theorem requires that no matter whether the subject is known to be treated, the support of the covariates is invariant in the directions that are informative to the treated outcome, and likewise for the control outcome. As the condition does not impose any restriction in the directions of covariates that are redundant to the outcomes, it relaxes the strong common support condition mentioned in Section $2$, and is referred to as the weak common support condition. A similar result for causal inference based on the outcome regressions can be found in Luo et al. (2016).

The theorem justifies that we can estimate the interested central subspaces using the observed data. Various SDR methods can be employed. We choose SIR as an example, and describe its implementation in detail in Section $5$. Given the resulting estimates, denoted by $\hat{\beta}_{\scriptscriptstyle 0}$ and $\hat{\beta}_{\scriptscriptstyle 1}$, respectively, we use $X^{\prime}\hat{\beta}_{\scriptscriptstyle 0}$ and $X^{\prime}\hat{\beta}_{\scriptscriptstyle 1}$ in the subsequent matching methods mentioned in the Introduction. As discussed above, the proposed matching approach outperforms the existing approaches based on either the original covariates or the propensity score in multiple ways. These advantages are further formulated in an asymptotic sense in Section $5$.

For the consistency of SIR, the covariates in each treatment group must satisfy the linearity condition (5), with $\Sigma_{\scriptscriptstyle X}$ being replaced with the conditional covariance matrix $\Sigma_{\scriptscriptstyle X|T}$. We allow $\Sigma_{\scriptscriptstyle X|T}$ to be stochastic, so that the proposed estimator is applicable, for example, when the covariates have different correlation structures in the treatment and control groups. When the reduced covariates in each group are elliptically distributed with a non-stochastic $\Sigma_{\scriptscriptstyle X|T}$ and equipped with an affinely invariant measure in the subsequent matching, it may be of interest to see whether the proposed approach is conditionally affinely invariant on the original covariates, which would then imply the conditionally equal percent bias reduction property (Rubin and Thomas, 1992, 1996). The answer depends on in which space such invariance is referred to: the approach is conditionally affinely invariant within the estimated central subspace but not within the true central subspace, as the estimated reduced covariates are invariant if the original covariates are rotated within the former, whereas the invariance does not hold if the rotations are in the latter.

4 Asymptotic study

Using the Abadie-Imbens estimator (3) on the reduced covariates, we now present some asymptotic results that illustrate the superiority of the proposed approach to those using the propensity score or the original covariates. For ease of presentation, we focus on the estimation of $ACE$. Similar results can be derived for $ACET$ with slight adjustments.

For simplicity, we fix the value of $m$ as the sample size grows for all the balancing scores. In the literature, its has been commonly realized that the asymptotic property of matching is difficult to tackle, if the estimated propensity score, rather than its true value, is used (Imbens and Rubin, 2015). Thus, in such studies, the true propensity score is assumed known a priori (Abadie and Imbens, 2006), unless a high-order kernel density function is used in kernel matching (Heckman et al., 1997). Following this convention, we assume both the propensity score and the central subspaces to be known a priori throughout the section. Although by doing this, we omit the effect of SDR estimation, the results in this section still help illustrate the advantage of using SDR for matching, and provide some guidelines when other SDR methods, such as the aforementioned SAVE and directional regression, are employed. In addition, as seen in the simulation studies, the effect of the SDR estimation is almost negligible in finite samples. For these reasons, we think that the results developed here are meaningful.

For regulation purpose, following Abadie and Imbens (2006), we adopt several assumptions as below:

Assumption 1**.**

The sample subjects are independent copies of $(X,T,Y(T))$. The distribution of $X$ is dominated by the Lebesgue measure on \mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}} with a compact support $\Omega(X)$. $E\{Y(t)|X\}$ is twice differentiable, and its Hessian matrix is Lipschitz continuous almost surely for $t=0,1$.

Assumption 2**.**

There exists $\tau>0$ such that the propensity score $\pi(X)$ satisfies that $\pi(X)\in(\tau,1-\tau)$ almost surely.

Assumption 1 is fairly general, and has been commonly adopted in the literature. Assumption 2 further regulates the strong common support condition in Section $2$. Similarly, we regulate the weak common support condition in the following, which relaxes Assumption 2.

Assumption 3**.**

Let $\pi(X^{\prime}\beta_{\scriptscriptstyle t})=P(T=1|X^{\prime}\beta_{\scriptscriptstyle t})$ for $t=0,1$. There exists $\tau_{\scriptscriptstyle r}>0$ such that $\pi(X^{\prime}\beta_{\scriptscriptstyle 0})<1-\tau_{\scriptscriptstyle r}$ and $\pi(X^{\prime}\beta_{\scriptscriptstyle 1})>\tau_{\scriptscriptstyle r}$ almost surely.

Under these assumptions, the following theorem shows the effect of the dimensionality of the balancing score on the convergence order of the resulting matching estimator. By conducting a more careful study on the variance of the estimator, it strengthens the result in Theorem 1 of Abadie and Imbens (2006) to allow a higher-dimensional balancing score for $\sqrt{n}$-consistent estimation.

Theorem 3**.**

Suppose that the ignorability assumption (1) and Assumption 1 hold. Let $\mu$ be $ACE$ and $\hat{\mu}_{\scriptscriptstyle X}$, $\hat{\mu}_{\scriptscriptstyle\pi}$ and $\hat{\mu}_{\scriptscriptstyle r}$ be the Abadie-Imbens estimator using the original covariates, the true propensity score, and the reduced covariates from the true central subspaces, respectively. Then,

(1) under Assumption 2,

[TABLE]

(2) under Assumption 3 and the linearity condition (5) for $X|T$,

[TABLE]

Proof.

We first show a mathematical property that will be used in the proof. Let $r,s>0$ be arbitrary positive real numbers. We show that, as $n\rightarrow\infty$,

[TABLE]

which is a complement of the well known result \lim_{\scriptscriptstyle n\rightarrow\infty}(1-rn^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1}}})^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n}}}=e^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-r}}}. To see why it holds, re-write the left-hand side as

[TABLE]

and divide it by \{1-(r+s)n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1}}}\}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n}}}, which converges to $\exp\{-(r+s)\}$. Then (9) is equivalent to

[TABLE]

Denote the left-hand side above by $\epsilon_{\scriptscriptstyle n}$. By simple algebra,

[TABLE]

which means that $\epsilon_{\scriptscriptstyle n}\rightarrow-(rs)$, and (9) holds. For simplicity, we prove the theorem for $\hat{\mu}_{\scriptscriptstyle X}$ and $m=1$, and denote $J_{\scriptscriptstyle 1}(i)$ by $J(i)$ for each subject $i$. The general case can be shown similarly. Denote $n_{\scriptscriptstyle 0}$ and $n_{\scriptscriptstyle 1}$ as the number of control and treated subjects, respectively. As we regard $T$ to be random, so are $n_{\scriptscriptstyle 0}$ and $n_{\scriptscriptstyle 1}$. By the law of large numbers, for $t=0,1$, $n_{\scriptscriptstyle t}n^{\scriptscriptstyle-1}\rightarrow(1-t)+(2t-1)P(T=1)$ in probability. Thus, for any $w\in\mathbb{R}$, O_{\scriptscriptstyle P}(n_{\scriptscriptstyle t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle w}}})=O_{\scriptscriptstyle P}(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle w}}}). For efficiency of presentation, without loss of generality, we treat each T^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}} as given, as well as $n_{\scriptscriptstyle 0}$ and $n_{\scriptscriptstyle 1}$, in the rest of the proof. Following equation (7) in Abadie and Imbens (2006), let $K(i)$ be the number of times subject $i$ is used to match the others. $\hat{\mu}_{\scriptscriptstyle X}$ can be decomposed as

[TABLE]

in which the first two terms on the right hand side are easily seen to be O_{\scriptscriptstyle P}(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/2}}}) with mean zero (Abadie and Imbens, 2006). The third term represents the bias in matching, which we denote by $B_{\scriptscriptstyle n}$. To show the result about $\hat{\mu}_{\scriptscriptstyle X}$, it suffices to show that E(B_{\scriptscriptstyle n})=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-2/p}}}) and \mathrm{var}(B_{\scriptscriptstyle n})=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-\min{1+2/p,6/p}}}}). By Assumption 1 and Taylor’s expansion,

[TABLE]

where g_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}} and g_{\scriptscriptstyle 2}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}} are the gradient and the Hessian matrix of $E\{Y(1-T)\mid X\}$ at X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}, and g_{\scriptscriptstyle 3}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i,J(i))}}} is a bounded random element. Thus, for the desired result about $E(B_{\scriptscriptstyle n})$ and $\mathrm{var}(B_{\scriptscriptstyle n})$, it suffices to show that, for an arbitrary pair of subjects $(i,j)$,

[TABLE]

in which $k,l\in\{1,2\}$, and v^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes 1}}}=v and vv^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes 2}}}=vv^{\mbox{\tiny{\sf T}}} for any real vector $v$. To show (4), we first suppose that both $i$ and $j$ belong to the same treatment group $T=t$. For any $a,b\in\Omega(X|T=t)$, let (u,v)=n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 1/p}}}(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(a))}}}-a,X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(b))}}}-b), in which $J(a)$ denotes the subject whose covariates are closest to a, and $J(b)$ likewise. From the proof of Theorem 1 in Abadie and Imbens (2006), let $f(\cdot)$ be the density function of random elements measurable with respect to \{X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}},i=1,\ldots,n\}, we have

[TABLE]

in which f^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}(a,u) is defined so as to make the equation holds. By Assumption $1$, f^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}(a,u) is bounded. As shown in Abadie and Imbens (2006), \{1-P(\|X-a\|\leq\|u\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}})\}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n_{\scriptscriptstyle 1-t}}}} converges to \exp[-2\pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p/2}}}\|u\|^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}f(a)/\{p\Gamma(p/2)\}]. Thus $u=O_{\scriptscriptstyle P}(1)$, which implies that X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(a))}}}-a=O_{\scriptscriptstyle P}(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}). Using the symmetry of \{1-P(\|X-a\|\leq\|u\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}})\}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n_{\scriptscriptstyle 1-t}}}} about the origin in \mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}, Abadie and Imbens (2006) further showed that E(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(a))}}}-a)=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-2/p}}}). By the compactness of $\Omega(X)$, such convergence is uniform for $a$, thus X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(i))}}}-X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}=O_{\scriptscriptstyle P}(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}) and E(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(i))}}}-X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}})=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-2/p}}}). Similar to (11), we further have

[TABLE]

in which P(\|X-a\|\leq\|u\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}\mbox{ or }\|X-b\|\leq\|v\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}) can be written as

[TABLE]

Obviously, $I\!I\!I=0$ for all large $n$. Let r_{\scriptscriptstyle n}=n_{\scriptscriptstyle 1-t}P(\|X-a\|\leq\|u\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}), and s_{\scriptscriptstyle n}=n_{\scriptscriptstyle 1-t}P(\|X-b\|\leq\|v\|n_{\scriptscriptstyle 1-t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/p}}}). Then

[TABLE]

for all large $n$, and $r_{\scriptscriptstyle n}\rightarrow r$ and $s_{\scriptscriptstyle n}\rightarrow s$, in which r=2\pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p/2}}}\|u\|^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}f(a)/\{p\Gamma(p/2)\} and s=2\pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p/2}}}\|v\|^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}f(b)/\{p\Gamma(p/2)\}. Since

[TABLE]

we have, for $k,l\in\{1,2\}$,

[TABLE]

in which \eta_{\scriptscriptstyle a,b}=\{1-r_{\scriptscriptstyle n}n_{\scriptscriptstyle 1-t}^{\scriptscriptstyle-1}-s_{\scriptscriptstyle n}n_{\scriptscriptstyle 1-t}^{\scriptscriptstyle-1}\}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n_{\scriptscriptstyle 1-t}}}}-(1-r_{\scriptscriptstyle n}n_{\scriptscriptstyle 1-t}^{\scriptscriptstyle-1})^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n_{\scriptscriptstyle 1-t}}}}(1-s_{\scriptscriptstyle n}n_{\scriptscriptstyle 1-t}^{\scriptscriptstyle-1})^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle n_{\scriptscriptstyle 1-t}}}}. By (9), $\eta_{\scriptscriptstyle a,b}\rightarrow-{n_{\scriptscriptstyle 1-t}}^{\scriptscriptstyle-1}rs\exp\{-(r+s)\}$. Since \|u\|^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}\|v\|^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}rs\exp\{-(r+s)\} is integrable on (u,v)\in\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}\times\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}}, similar to Abadie and Imbens (2006), we have E\{u^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes{k}}}}(v^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes l}}})^{\prime}\}-E(u^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes k}}})E^{\prime}(v^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes l}}})=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1}}}). Next, suppose that $i$ and $j$ are control and treated subjects, respectively. Conditioning on (X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}},X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(j)}}})=(a,b), let u=n_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 1/p}}}(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(a))}}}-a) and v=n_{\scriptscriptstyle 0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 1/p}}}(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(b))}}}-b), we have

[TABLE]

in which $\delta_{\scriptscriptstyle w}(x)$ is the Dirac function of $x$ such that it is zero whenever $x\neq w$ and \textstyle\int_{\scriptscriptstyle\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle dim(x)}}}}\delta_{\scriptscriptstyle w}(x)h(x)dx=h(w) for any function $h$ of $x$. Similar to the above, we can show that \mathrm{cov}(u^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes{k}}}},v^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes l}}})=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1}}}). By the compactness of $\Omega(X)$, this convergence is uniform on $(a,b)$, which means that \mathrm{cov}\{(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(i))}}}-X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}})^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes k}}},(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(J(j))}}}-X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(j)}}})^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle\otimes l}}}\}=O(n^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1-2/p}}}). Hence (4) holds, which completes the proof. ∎

The theorem suggests that, compared to the variance, the bias of the matching estimator is more sensitive to the dimension of the balancing score. In fact, if we assume a higher-order differentiability of the outcome regression function in Assumption 1, then a higher-dimensional balancing score will be allowed for the resulting matching estimator to have variance of order $O(n^{\scriptscriptstyle-1})$. In particular, we speculate that if the outcome regression function is smooth, then the variance of the matching estimator is always $O(n^{\scriptscriptstyle-1})$, regardless of the dimension of the balancing score used. The simulation studies reveal this point.

From this theorem, whenever the dimension of the balancing score is greater than four, it determines the convergence order of the resulting estimator of $ACE$. Intuitively, this is caused by the fact that when the sample size is fixed and more components are added to the covariates, the empirical distribution of the covariates becomes more sparse, so that each subject $i$ is less similar to those in its neighborhood $J_{\scriptscriptstyle m}(i)$. From this point of view, the propensity score is optimal among all, as it always generates a $\sqrt{n}$-consistent estimator. By contrast, using the original covariates is the worst choice.

When the central subspaces for both the treatment and control groups are at most four-dimensional, the reduced covariates from SDR become a comparable choice to the propensity score, as they also provide a $\sqrt{n}$-consistent estimator of $ACE$. This is very common in applications. For example, the central subspace is one-dimensional under the popular single-index model. In this special case, the following theorem suggests that the reduced covariates can further outperform the propensity score, in the sense that they produce an asymptotically stabler matching estimator. For regulation purpose, in addition to Assumption 1, we adopt the following assumption:

Assumption 4**.**

For $t=0,1$, almost surely, E\{Y^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle w}}}(t)|X\} is Lipschitz continuous when $w=2$ and uniformly bounded when $w=4$, and $\mathrm{var}\{Y(t)|X\}>0$.

Theorem 4**.**

Suppose the ignorability assumption (1), Assumptions 1, 3 and 4 hold. If both ${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}$ and ${\mathcal{S}}_{\scriptscriptstyle Y(1)|X}$ are one-dimensional, then

[TABLE]

where $V_{\scriptscriptstyle X}=\mathrm{var}[E\{Y(1)-Y(0)|X\}]$, $h(\pi(X^{\prime}\beta_{\scriptscriptstyle 0}),0)=1-\pi(X^{\prime}\beta_{\scriptscriptstyle 0})$ and $h(\pi(X^{\prime}\beta_{\scriptscriptstyle 1}),1)=\pi(X^{\prime}\beta_{\scriptscriptstyle 1})$. In addition, $\mathrm{var}(\hat{\mu}_{\scriptscriptstyle r})\leq\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})$ for all large $n$ when the latter exists, and the equality holds in one of the following cases:

• (i)

$T\perp\perp X$, and $E\{Y(t)|X\}\equiv E\{Y(t)\}$ for $t=0,1$.

• (ii)

${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}={\mathcal{S}}_{\scriptscriptstyle Y(1)|X}$, $T\perp\perp X|X^{\prime}\beta_{\scriptscriptstyle 0}$, and $E\{Y(t)|X^{\prime}\beta_{\scriptscriptstyle t}\}=E\{Y(t)|\pi(X^{\prime}\beta_{\scriptscriptstyle t})\}$ for $t=0,1$.

Proof.

The form of the variance follows directly from Theorem $5$ of Abadie and Imbens (2006) and the sufficiency of $X^{\prime}\beta_{\scriptscriptstyle t}$ for $Y(t)|X$. Let \sigma_{\scriptscriptstyle t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}(R)=\mathrm{var}\{Y(t)|R\} for any random element $R$. To see that $\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})\geq\mathrm{var}(\hat{\mu}_{\scriptscriptstyle r})$ for all large $n$, following Theorem $5$ of Abadie and Imbens (2006), under Assumption 2, we have

[TABLE]

We write $n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})=V_{\scriptscriptstyle\pi}+I_{\scriptscriptstyle\pi}+I\!I_{\scriptscriptstyle\pi}$ and $n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})=V_{\scriptscriptstyle X}+I_{\scriptscriptstyle r}+I\!I_{\scriptscriptstyle r}$, and additionally introduce $I_{\scriptscriptstyle X}$ and $I\!I_{\scriptscriptstyle X}$, in which

[TABLE]

Then the inequality $n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})\geq n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle r})$ can be implied if

[TABLE]

Let $\mu_{\scriptscriptstyle c,t}(X)=E\{Y(t)|X\}-E\{Y(t)|\pi(X)\}$ for $t=0,1$. By definition,

[TABLE]

Let \pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}(X)=\pi(X)/\{1-\pi(X)\}, the logit function of $\pi(X)$. By the triangle inequality,

[TABLE]

By plugging it back to (13), we have

[TABLE]

On the other hand, we have

[TABLE]

Hence $V_{\scriptscriptstyle X}-V_{\scriptscriptstyle\pi}\leq I_{\scriptscriptstyle\pi}-I_{\scriptscriptstyle X}$, which implies the first inequality in (12). Because E\{\sigma_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}(X)|\pi(X)\}=E[\mathrm{var}\{Y(1)|X\}|\pi(X)]\leq\mathrm{var}\{Y(1)|\pi(X)\}=\sigma_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}(\pi(X)), we have

[TABLE]

Similarly, we can show the corresponding inequality for the part when $t=0$ with $1-\pi(X)$ in place of $\pi(X)$. Thus $I\!I_{\scriptscriptstyle\pi}\geq I\!I_{\scriptscriptstyle X}$. Next, let \phi:\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle+}}}\rightarrow\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle+}}} be that $\phi(x)=x^{\scriptscriptstyle-1}$. Then $\phi$ is a convex function. By Jensen’s inequality, we have $\pi(X^{\prime}\beta_{\scriptscriptstyle 1})^{\scriptscriptstyle-1}=[E\{\pi(X)|X^{\prime}\beta_{\scriptscriptstyle 1}\}]^{\scriptscriptstyle-1}\leq E\{\pi^{\scriptscriptstyle-1}(X)|X^{\prime}\beta_{\scriptscriptstyle 1}\}$, which, together with that \sigma_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}(X)=\sigma_{\scriptscriptstyle 1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}(X^{\prime}\beta_{\scriptscriptstyle 1}), implies that

[TABLE]

Similarly, we can show the corresponding inequality for the part when $t=0$ with $1-\pi(X)$ in place of $\pi(X)$. Hence $I_{\scriptscriptstyle X}\geq I_{\scriptscriptstyle r}$. Finally, let \psi:\mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle+}}}\rightarrow\mathbb{R} be that $\psi(x)=x^{\scriptscriptstyle-1}-x$, then $\psi$ is also a convex function. Thus the fourth inequality in (12) can be shown similarly. From the arguments above, the equality between the limits of $n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle\pi})$ and $n\mathrm{var}(\hat{\mu}_{\scriptscriptstyle r})$ holds if and only if all the four equalities in (12) hold, which means that

• (a)

\mu_{\scriptscriptstyle c,1}(X){\pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}(X)}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/2}}}=-\mu_{\scriptscriptstyle c,0}(X){\pi^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}(X)}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 1/2}}},

• (b)

$E\{Y(t)|X\}=E\{Y(t)|\pi(X)\}$ for $t=0,1$,

• (c)

$\pi(X)=\pi(X^{\prime}\beta_{\scriptscriptstyle 1})=\pi(X^{\prime}\beta_{\scriptscriptstyle 0})$.

Note that (b) is equivalent to that $\mu_{\scriptscriptstyle c,t}(X)=0$ for $t=0,1$, which implies (a). If ${\mathcal{S}}_{\scriptscriptstyle Y(0)|X}={\mathcal{S}}_{\scriptscriptstyle Y(1)|X}$, then (c) is equivalent to that $\pi(X)=\pi(X^{\prime}\beta_{\scriptscriptstyle 1})$. Otherwise, it means that $\pi(X)$ is a constant, which, together with (b), indicates that $E\{Y(t)|X\}\equiv E\{Y(t)\}$ for $t=0,1$. This completes the proof. ∎

Case (i) in this theorem means that the covariates are redundant to both the propensity score and the outcome regressions. That is, the treatment is assigned completely at random, and both outcomes have constant mean over the subjects. However, the covariates may affect the other aspects of the outcomes, such as the variance, so that the central subspaces do not vanish. Case (ii) in the theorem means that both the treatment assignment and the outcomes are affected by the same linear combination of the covariates, and furthermore, that the propensity score is sufficient for the outcome regressions.

Because both cases can be easily violated in practice, the theorem justifies that when the reduced covariates from SDR are univariate and known a priori, they are likely to outperform the propensity score, in terms of producing an asymptotically stabler estimator in the subsequent matching.

The underlying reason for the superiority of the reduced covariates from SDR, is their sufficiency for modeling the outcomes. In this sense, they are comparable to the prognostic score (Hansen, 2008), defined as $E\{Y(0)|X\}$, which is also sufficient subject to the conditional independence between $Y(1)$ and $X$ given $Y(0)$. Following the arguments in Theorem 3 and Theorem 4, this prognostic score can further outperform the proposed reduced covariates in reducing both the bias and the variance of the matching estimator. However, the additional assumption for the prognostic score involves the joint distribution of both potential outcomes, and its estimation usually relies on parametric models, neither of which can be readily justified in practice.

The reduced covariates from SDR are also comparable to the aforementioned active set in de Luna et al. (2011), in the sense that they both pre-assume a sparsity structure in the data. As mentioned before, the cardinality of the active set is always greater than or equal to the dimensions of the central subspaces, which, by Theorem 3, makes the subsequent matching potentially less effective.

To further reduce the size, de Luna et al. (2011) used the above-mentioned active set in place of the original covariates to further conduct variable selection in modeling the treatment assignment, and to remove the redundant covariates. The cardinality of the reduced set is possibly smaller than the dimensions of the central subspaces. However, as the reduced set is no longer sufficient for modeling the outcomes, Theorem 4 suggests that it may enlarge the variance of the subsequent matching estimator. Nonetheless, this procedure motivates us to develop a sequential SDR procedure, which can be useful when there is need to further reduce the working central subspaces. For continuity of context, we leave more details to the Discussion.

5 Sample estimation and implementation

We now present the details in using SIR to estimate the central subspaces, and then using the resulting reduced covariates in matching to estimate the causal effects.

Following the notations in Section 2, for standardized covariates with zero mean and identity covariance matrix, SIR estimates the central subspace through an intermediate parameter $M=E\{E(X|W)E(X^{\prime}|W)\}$, called the candidate matrix. Given a matrix estimate $\hat{M}$, the central subspace is estimated as the linear space spanned by the leading eigenvectors of $\hat{M}$. Because no model is assumed on $E(X|W)$, to construct $\hat{M}$, Li (1991) proposed the slicing strategy, which divides the sample into $H$ slices according to the similarity in $W$. $E(X|W)$ is then estimated as the sample mean of $X$ within each slice, and $\hat{M}$ is given by the usual sample moments.

Given the reduced covariates, we are free to use any existing matching method. As an illustration, here we use the Abadie-Imbens estimator with the Mahalanobis distance (Rosenbaum and Rubin, 1985) to measure the similarity of subjects. For a random vector $Z$ with the sample values \{Z^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}},i=1,\ldots,n\}, let $\hat{\Sigma}_{\scriptscriptstyle Z}$ be its sample covariance matrix. The Mahalanobis distance between Z^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}} and Z^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(j)}}} is

[TABLE]

When $Z$ is used for matching, for each subject $i$, the matched set $J_{\scriptscriptstyle m}(i)$ in (3) contains the $m$ subjects in the other treatment group whose $Z$ values are closest to Z^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}} in terms of Mahalanobis distance. After calculating the imputed values by (3), we estimate $ACE$ by

[TABLE]

In summary, the proposed matching approach can be implemented in the following steps:

1. Step 1.

SDR: estimate the reduced covariates in the control group as follows.

• •

Standardize $X$ to be \hat{\Sigma}_{\scriptscriptstyle X|T=0}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/2}}}\{X-\hat{E}(X|T=0)\}, where $\hat{E}(X|T=0)$ and $\hat{\Sigma}_{\scriptscriptstyle X|T=0}$ are the sample mean and the sample covariance matrix of $X$ in the control group, respectively. For simplicity, we still denote the standardized variable by $X$, if no ambiguity is caused.

• •

For $j=1,\ldots,H$, let $q_{\scriptscriptstyle j}$ be the $(j/H)$th sample quantile of \{Y^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(0):T^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}=0\}. Construct the $j$th slice as $S_{\scriptscriptstyle j}(Y(0))=I(q_{\scriptscriptstyle j-1}<Y(0)\leq q_{\scriptscriptstyle j})$, where $I(\cdot)$ is the indicator function and $q_{\scriptscriptstyle 0}$ is a constant that is less than any observed Y^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(0). Estimate E(X|S_{\scriptscriptstyle j}(Y(0))=1,T^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}=0) by

[TABLE]

• •

Estimate $M_{\scriptscriptstyle 0}=E[E\{X|Y(0),T=0\}E\{X^{\prime}|Y(0),T=0\}]$ by \hat{M}_{\scriptscriptstyle 0}=\textstyle\sum_{\scriptscriptstyle j=1}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle H}}}\hat{\mu}_{\scriptscriptstyle X,j}\hat{\mu}_{\scriptscriptstyle X,j}^{\prime}/H, and use the sequential tests to determine the rank of $M_{\scriptscriptstyle 0}$. Let $\hat{r}(0)$ be the selected rank and $\hat{\beta}_{\scriptscriptstyle 0}$ be the $\hat{r}(0)$ eigenvectors of $\hat{M}_{\scriptscriptstyle 0}$ with the largest eigenvalues. 2. Step 2.

Match each subject in the treatment group with those in the control group: for subject $i$, let {\mathcal{D}}_{\scriptscriptstyle i}=\{D(X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)\prime}}}\hat{\beta}_{\scriptscriptstyle 0},X^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(j)\prime}}}\hat{\beta}_{\scriptscriptstyle 0}):T^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(j)}}}=0\}, where $D$ is the Mahalanobis distance in (14), and $J_{\scriptscriptstyle m}(i)$ be the set that corresponds to the $m$ smallest values in ${\mathcal{D}}_{\scriptscriptstyle i}$. Use (3) to compute \hat{Y}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(0). 3. Step 3.

Conduct Step 1 for the treatment group, and conduct Step 2 to compute \hat{Y}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle(i)}}}(1) for each subject $i$ in the control group. 4. Step 4.

Estimate $ACE$ by (15).

In practice, the number of slices $H$ is usually fixed at a small constant, such as $5$ or $10$. Although it imposes the constraint that the reduced covariates derived from SIR are at most $(H-1)$-dimensional, the issue can be easily addressed: if one finds that the estimated reduced dimension $\hat{r}(0)$ or $\hat{r}(1)$ coincides with $H-1$, then one can raise the value of $H$ to see whether more directions in the central subspace can be recovered.

When the causal effect of interest is $ACET$, the implementation can be adjusted by skipping step 3 and using the results in step 2 to compute

[TABLE]

6 Simulation Studies

In this section, we illustrate the effectiveness of SDR for matching, by applying the reduced covariates to the Abadie-Imbens estimator. For comparison, we also apply the estimated propensity score derived from CBPS (Imai and Ratkovic, 2014), the true propensity score, the original covariates, and the active set of covariates (de Luna et al., 2011) to the same estimator, respectively. The active set is assumed known a priori. To assess the effect of SDR estimation in the proposed method, we also apply the reduced covariates induced by the true central subspaces, which is oracle, just like the active set of covariates. In addition, we include genetic matching on the original covariates in the comparison, for the reason that the method has the potential to detect and use the sparsity structure in the data.

Because the linearity condition (5) required for the consistency of SIR in each treatment group can be violated in practice, we do not treat the condition as granted. Instead, we let the covariates be elliptically distributed conditional on the treatment assignment, elliptically distributed marginally, and composed of a mixture of continuous and discrete components, sequentially in this section.

6.1 Case 1: Elliptical covariates conditional on treatment

In the first scenario, we generate $T$ marginally from a Bernoulli distribution with mean $0.5$, and then generate $X|T=t$ from $N(\mu_{\scriptscriptstyle t},\Sigma_{\scriptscriptstyle t})$, where the $(i,j)$th entry of $\Sigma_{\scriptscriptstyle t}$ is \delta_{\scriptscriptstyle t}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle|i-j|}}}. We set $\mu_{\scriptscriptstyle 0}$ at the origin of \mathbb{R}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle p}}} and $\delta_{\scriptscriptstyle 0}$ at $0.2$, and let $\mu_{1}$ and $\delta_{1}$ vary with the models. When $\delta_{\scriptscriptstyle 1}=0.2$, the distribution $X|T=t$ depends on $t$ only through the conditional mean $\mu_{t}$, which is the case discussed in Rubin and Thomas (1992). The true propensity scores can be easily calculated using Bayes’ Theorem.

We set $n=500$, $p=10$, and study the following four models:

• I.

Y(t)=(t+2)(X_{\scriptscriptstyle 1}+1.5)^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}+t+\varepsilon, $\mu_{\scriptscriptstyle 1,k}=\mu_{\scriptscriptstyle 0}$, $\delta_{\scriptscriptstyle 1}=0.5$;

• I​I.

$Y(t)=2\sin\{0.5(X_{\scriptscriptstyle 1}+X_{\scriptscriptstyle 2}-X_{\scriptscriptstyle 3})\}+t+\varepsilon$, $\mu_{\scriptscriptstyle 1,k}=I(k\leq 3)/3$, $\delta_{\scriptscriptstyle 1}=0.5$;

• I​I​I.

$Y(t)=X_{\scriptscriptstyle 1}+X_{\scriptscriptstyle 2}+X_{\scriptscriptstyle 3}+t(X_{\scriptscriptstyle 4}+X_{\scriptscriptstyle 5})+\varepsilon$, \mu_{\scriptscriptstyle 1,k}\equiv p^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle-1/2}}}, $\delta_{\scriptscriptstyle 1}=0.2$;

• I​V.

Y(t)=3\sin\{X^{\prime}\Sigma_{\scriptscriptstyle 1}^{\scriptscriptstyle-1}(\mu_{\scriptscriptstyle 1}-\mu_{\scriptscriptstyle 0})/3\}+t\,\{X^{\prime}\Sigma_{\scriptscriptstyle 1}^{\scriptscriptstyle-1}(\mu_{\scriptscriptstyle 1}-\mu_{\scriptscriptstyle 0})\}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}/3+\varepsilon, $\mu_{\scriptscriptstyle 1,k}=c_{\scriptscriptstyle 1}k$, $\delta_{\scriptscriptstyle 1}=0.2$.

Here $\epsilon$ is an independent error generated from N(0,0.5^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle 2}}}), $\mu_{\scriptscriptstyle 1,k}$ is the $k$th component of $\mu_{\scriptscriptstyle 1}$ for $k=1,\ldots,p$, and $c_{\scriptscriptstyle 1}$ is the constant that makes $\|\mu_{\scriptscriptstyle 1}\|_{\scriptscriptstyle 2}=0.5$ in Model I​V. Obviously, the ignorability assumption (1) is satisfied in all the models.

In Model I​I, the causal effects are constant among all the subjects. In Models I and I​I, $\delta_{\scriptscriptstyle 1}$ differs from $\delta_{\scriptscriptstyle 0}$, so that the correlation structure of the covariates varies between treatment groups. From the view of SDR, in Model I​I​I, the two central subspaces differ, and in all the models the central subspaces are one-dimensional. Thus, by Theorem 4, the reduced covariate from oracle SDR should outperform the true propensity score in Models I, I​I and I​I​I, in the sense that the induced Abadie-Imbens estimator has a smaller asymptotic variance. In Model I​V, the true propensity score is an invertible function of $X^{\prime}\Sigma_{\scriptscriptstyle 1}^{\scriptscriptstyle-1}(\mu_{\scriptscriptstyle 1}-\mu_{\scriptscriptstyle 0})$, which is also the reduced covariate from oracle SDR. Therefore, the model falls into case (ii) of Theorem 4 indicating that the two balancing scores are equally good for matching asymptotically.

The strong common support condition is satisfied at the population level in all four models. However, because the dimension $p$ is moderately large, the covariates have a relatively sparse support at the sample level, which makes the actual realization of the condition questionable. In contrast, because the reduced covariate from SDR is univariate in all the models, the issue is less severe. In fact, depending on how different the central subspaces are from $\mu_{\scriptscriptstyle 1}-\mu_{\scriptscriptstyle 0}$, $\Sigma_{\scriptscriptstyle 0}$ and $\Sigma_{\scriptscriptstyle 1}$, the weak common support condition is realized almost perfectly in Models I and I​I and to a reasonable extent in the other models.

We set $m=1$, and estimate both $ACE$ and $ACET$ in each model. The results for $ACE$ are summarized in Table 1 and the results for $ACET$ are displayed in Table 1 in the Appendix B (Supplementary Material), respectively, which include the evaluation of the bias and the standard deviation of each estimator, based on $1000$ independent runs. For readers’ convenience, we also record the root mean square error (RMSE) of each estimator, which can be calculated from the bias and the standard deviation, subject to rounding errors.

From these tables, the reduced covariates from SDR estimation substantially outperform the original covariates for matching, in all the aspects in Models I to I​I​I and in terms of reduced bias in Model I​V. They are comparable with the active set of covariates in Models I and I​I, and superior to the latter on bias reduction in Models I​I​I and I​V. The improvement is due to the relaxation of the “curse of dimensionality”, for the reason that all these three balancing scores are sufficient for the outcome regressions, and that the cardinality of the active set is small in Models I and I​I, and is relatively large in Models I​I​I and I​V. The genetic matching is generally superior to the Abadie-Imbens estimator based on the original covariates, but is inferior to the same estimator based on the reduced covariates from SDR estimation. This is a sign that at the current sample size, the sparsity structure of the data can be detected partially, but not completely, by genetic matching.

Interestingly, the true propensity score is slightly inferior to its estimator in terms of a larger variance in the resulting matching. The phenomenon has also been noticed by multiple authors, see Lunceford and Davidian (2004). Compared to the estimated propensity score, the reduced covariates from SDR estimation makes matching more stable in Models I, I​I, and I​I​I, and equally consistent in Model I​V, which complies with the theoretical anticipation from Theorem 4. Compared to the oracle SDR, the SDR estimation does cost a larger variation in matching, although the price is nearly negligible.

6.2 Case 2: Elliptical covariates in the merged sample

When applying SIR to data with continuous but non-elliptically distributed covariates, it has been a common practice to transform each covariate to be univariate normally distributed, and assume that the transformed covariates have a jointly normal distribution. In causal inference, if the researcher believes that the covariates affect the outcome in similar patterns before and after the treatment assignment, it will be desired to transform the covariates uniformly across the treatment groups, which corresponds to an elliptically distributed $X$ but not necessarily elliptically distributed $X|T$. Consequently, the linearity condition (5) on $X|T$ can be violated.

To study the performance of the proposed approach in this case, we still set $n=500$, $p=10$, and generate $X$ from $N({\mathbf{0}},I_{\scriptscriptstyle p})$. Because the outcome regressions are similar in different treatment groups in Models I and I​I, given $X$, we generate the outcomes from these two models, and use the corresponding propensity score to generate $T$ from an independent Bernoulli distribution. The new models satisfy the ignorability assumption (1), and differ from their counterparts in the previous subsection only in the marginal distribution $X$. We label them as Model I{}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}} and Model I​I{}^{\raisebox{1.2pt}{\mbox{\scriptscriptstyle*}}}, respectively. The results from $1000$ independent runs are summarized in Table 2.

The results are similar to the previous subsection. The reduced covariates from SDR estimation still outperform the original covariates and the propensity score, no matter whether the latter is estimated or known a priori. They are generally comparable to the reduced covariates from oracle SDR and the active set of covariates, indicating that the price of estimating the central subspaces is nearly negligible, and that the dimensionality of the balancing score does not harm the accuracy of matching when it is less than four, as expected in Theorem 3.

6.3 Case 3: Covariates that contain discrete components

As mentioned in Section $2$, when $X$ is non-elliptically distributed, the linearity condition can be approximately satisfied, as long as the dimension $p$ is moderately large. Consequently, the proposed matching approach can still be reasonably effective in finite samples. In this subsection, we examine its performance when $X$ is a mixture of elliptically distributed covariates and binary covariates.

The simulation setup follows Lee et al. (2010). First, we generate ten covariates $(X_{\scriptscriptstyle 1},\ldots,X_{\scriptscriptstyle 10})$, in which $(X_{\scriptscriptstyle 1},\ldots,X_{\scriptscriptstyle 3})$ are only associated with $Y$, $(X_{\scriptscriptstyle 4},\ldots,X_{\scriptscriptstyle 6})$ are only associated with $T$, and $(X_{\scriptscriptstyle 7},\ldots,X_{\scriptscriptstyle 10}$) are associated with both $Y$ and $T$. Six covariates, $(X_{\scriptscriptstyle 1},X_{\scriptscriptstyle 3},X_{\scriptscriptstyle 5},X_{\scriptscriptstyle 6},X_{\scriptscriptstyle 8},X_{\scriptscriptstyle 9})$, have Bernoulli distribution with mean 0.5 marginally, and the others follow standardized normal distribution marginally. The correlation within each pair of covariates is zero, except that:

[TABLE]

$T$ is generated from the conditional distribution $T|X$, in which the propensity score satisfies:

[TABLE]

The observed outcome $Y(T)$ is generated from the linear model

[TABLE]

where $\epsilon$ follows $N(0,0.1^{2})$. The function $g$ and the coefficient vectors $\alpha$ and $\omega$ are specified in the each of scenarios listed in the following. These scenarios differ in the degree of linearity and additivity in the propensity score (16). For more details, see Lee et al. (2010).

A:

Additivity and linearity (main effects only);

B:

Mild non-linearity (one quadratic term);

C:

Moderate non-linearity (three quadratic terms);

D:

Mild non-additivity (four two-way interaction terms);

E:

Mild non-additivity and non-linearity (three two way interaction terms and one quadratic term);

F:

Moderate non-additivity (ten two-way interaction terms);

G:

Moderate non-additivity and non-linearity (ten two-way interaction terms and three quadratic terms).

Again, $1000$ independent samples are generated, and the Abadie-Imbens estimators based on different balancing scores are used to estimate the causal effects. The results for $ACE$ are summarized in Table 3 and the results for $ACET$ are displayed in Table 2 in the Appendix B (Supplementary Material). Same as before, the proposed approach yields smaller bias and variance than the existing methods in all the scenarios, and is only slightly worse than that based on the oracle SDR.

7 Data Application

We now illustrate the proposed methodology using a well-known dataset: LaLonde dataset. This dataset has been analyzed by LaLonde (1986) and Dehejia and Wahba (1999) to evaluate the causal effect of a labor training program called National Supported Work (NSW) Demonstration on earnings for job-seekers, who had economic and social problems before the enrollment in the program. In this study, the outcome is the individuals’ earnings in 1978 (Ee78). The treatment variable is the indicator for the enrollment in the labor training program (Treat). The ten potential confounders are: age (Age), years of schooling (Educ), indicator for Blacks (Black), indicator for Hispanics (Hisp), indicator for being married (Married), indicator for high school diploma (Nodegr), real earnings in 1974 (Re74), real earnings in 1975 (Re75), indicator variable for earnings in 1974 being zero (U74), and indicator variable for earnings in 1975 being zero (U75). Here we consider the subset of the original dataset referred as “CPS-3” by Dehejia and Wahba (1999), which contains 185 subjects who participated in NSW and 429 controls who did not participate in NSW.

Since the focus is to examine the impact of the labor training program on postintervention earnings for those who are eligible for the program, the parameter of interest is $ACET$. Dehejia and Wahba (1999) suggested a one-to-one matching approach with replacement, based on the estimated propensity score from parametric models. To apply the proposed matching approach to the dataset, we apply SIR in the control group to obtain the reduced covariates. The sequential test for SIR (see the Appendix A) indicates that the dimension of the reduced covariates is two. A scatterplot between each reduced covariate and the outcome in the control group is shown in Figure 1. From these plots, both reduced covariates clearly affect the outcome marginally, suggesting that two is a reasonable choice for the reduced dimension.

To check the quality of matching, we draw the boxplot for each reduced covariate, grouped by the treatment assignment. For comparison, we also draw the boxplots for the estimated propensity scores from the three parametric models (See Appendix C in the Supplementary Material) used in Dehejia and Wahba (1999). From the upper panel of Figure 2, when any of the estimated propensity scores is used, the majority of the treated subjects will be matched with a small number of extreme outliers in the control group, which result in loss of power in the subsequent causal effect estimation. Besides, the sparsity of the empirical support of the outliers will also cause non-negligible difference within the matched pairs, and lead to large bias in the estimated casual effect. By contrast, from the lower panel of Figure 2, the issue is much less severe when the reduced covariates from SIR are used for matching.

We further plot the histogram of the two reduced covariates (Figure 3) and the estimated propensity scores (Figure 4) between the treatment and the control group. For propensity score-based approaches, there are a few bins in which the control group is much smaller than the treatment group, which means it will be hard to find matches for the treated subjects in those bins. However, for the proposed approach, the control group is almost always larger than the treatment group, which indicates there is enough overlapping for the estimation of the causal effect to be reliable.

The proposed method estimates $ACET$ to be 205 suggesting that the job-seeking program increases the individual earnings at 1978. This is consistent with the conclusion from an experimental study in Dehejia and Wahba (1999). In contrast, if we apply the Abadie-Imbens estimator based on the original covariates or the estimated propensity scores, the Abadie-Imbens estimator estimates $ACET$ to be -361 (original covariates), -835 (PS1), -528 (PS2), -304 (PS3), respectively, which leads to the opposite conclusion.

8 Discussion

When the dimension of the proposed reduced covariates in either treatment group is greater than four, as Theorem 3 shows, it is desirable to further reduce the covariates to enhance the estimation accuracy of matching. Following the discussion in Section $4$, we can develop sequential SDR by using the outcome and the treatment assignment alternatively as the response variable in dimension reduction. The procedure gives lower-dimensional covariates compared to the proposed approach, when certain sparsity structure exists between the treatment assignment and the covariates. Alternatively, as the proposed reduced covariates are linear combinations of the original covariates, they do not need to be minimal sufficient statistics. Hence, instead of assuming additional sparsity structure in the data, one can also search for finer sufficient statistics in subsequent matching. This can be done in a data-adaptive manner, for example, by using genetic matching based on the proposed reduced covariates.

When the dimension of the covariates is relatively large compared to the sample size, it has been commonly recognized that the conventional SDR methods such as SIR lose effectiveness. In this case, we recommend to use sparse SDR, which constructs the central subspaces on the active set of covariates rather than the original covariates. For more detail, see Chen et al. (2010).