Telescope matching
2022-03-21
Source:vignettes/articles/telescope_matching.Rmd
telescope_matching.RmdThis vignette illustrates the telescope matching method developed in Blackwell and Strezhnev (2021) as implemented in the telescope_matching routine. One drawback of sequential g-estimation or other model-based estimators is that they depend on correctly specifying two regression models: one for the effect of the mediator given treatment and pre-/post-treatment covariates and another for the effect of treatment given pre-treatment covariates. Telescope matching provides a more flexible and less model-dependent estimation strategy for controlled direct effects when the treatment and mediator are binary. It combines a two-stage matching procedure to impute the unobserved counterfactuals with a bias-correction to account for biases induced by imperfect matches.
Estimation by telescope matching
Telescope matching provides an alternative, less model-dependent approach to estimating the average controlled direct effect when both treatment and mediator are binary. The average controlled direct effect of interest is the effect of treatment versus control (\(A_i = 1\) vs. \(A_i = 0\)) holding constant the mediator at \(0\) (\(M_i = 0\)). \[ ACDE = E[Y_i(1, 0) - Y_i(0,0)] \]
Identification still depends crucially on the sequential unconfoundedness assumption Assumption 1, which states that there are no unobserved confounders of \(A_i\) and \(Y_i\) given pre-treatment covariates \(X_i\) and that there are no unobserved confounders of \(M_i\) and \(Y_i\) given treatment \(A_i\), pre-treatment covariates \(X_i\) and post-treatment covariates \(Z_i\). We can understand the process of estimating the ACDE as an imputation problem. Our estimator of the ACDE, \(\hat{\tau}\) is simply the average of imputed potential outcomes for each unit in the sample. \[ \hat{\tau} = \frac{1}{N}\sum_{i=1}^N \left(\widehat{Y}_i(1,0) - \widehat{Y}_i(0,0)\right) \]
To obtain the imputations \(\widehat{Y}_i(1,0)\) and \(\widehat{Y}_i(0,0)\) for each unit, the telescope matching approach proceeds in two stages. The first step matches each unit with \(M_i = 1\) to some user-specified number of units with \(M_i = 0\) that share treatment status \(A_i\) and are similar in both pre-treatment covariates \(X_i\) and post-treatment covariates \(Z_i\) as measured by some distance metric (in telescope_match we implement the Mahalanobis distance). Let \(\mathcal{J}_L^m(i)\) denote the set of units matched to unit \(i\). We define the following imputation for each unit’s potential outcome fixing \(M_i\) to 0.
\[ \widehat{Y}_i(A_i, 0) =\begin{cases} Y_i & \text{ if } M_i = 0 \\ \frac{1}{L} \sum_{\ell\in \mathcal{J}^m_L(i)} Y_{\ell} & \text{ if } M_i = 1\end{cases} \]
In the second step, we match each unit to some number of units of the opposite treatment status with similar values of pre-treatment covariates \(X_i\) Let \(\mathcal{J}_L^a(i)\) denote the set of units matched to unit \(i\) such that \(A_j = 1 - A_i\) for all \(j \in \mathcal{J}_L^a(i)\). We then use the first-stage imputations for either unit \(i\) or its second-stage matches to impute the potential outcomes for each unit under treatment and control fixing the mediator to \(0\).
\[ \widehat{Y}_{i}(a,0) = \begin{cases} \widehat{Y}_{i}(A_i,0) & \text{if}\ A_i = a\\ \frac{1}{L} \sum_{j\in \mathcal{J}^a_L(i)} \widehat{Y}_{j}(A_j,0) & \text{if}\ A_i = 1-a \end{cases} \]
Since matching is done with replacement, units may be used multiple times. We define \(K^m_L(i) = \sum_{k=1}^N \mathbb{I}\{i \in \mathcal{J}^m_L(k)\}\) as the number of times that unit \(i\) is used as a match in stage and \(K^a_L(i) = \sum_{j=1}^N \mathbb{I}\{i \in \mathcal{J}^a_L(j)\}\). Moreover, since units with \(M_i = 0\) contribute indirectly to second stage matches we define \(K^{am}_L(i) = \sum_{j=1}^N \mathbb{I}\{i \in \mathcal{J}^m_L(j)\}K^a_L(j)\) to denote the number of times a unit with \(M_i = 0\) matched to a unit with \(M_i = 1\) is implicitly used as a match in the second stage. This allows us to re-write the simple matching estimator as a weighted average \(\widehat{\tau} = N^{-1} \sum_{i=1}^N (2A_i - 1)(1 - M_i)W_iY_i\) where the weight is defined as
\[ W_i = 1 + \frac{K^a_L(i)}{L} + \frac{K^m_L(i)}{L} + \frac{K^{am}_L(i)}{L^2} \]
These weights can be used as a diagnostic for assessing the variance of the estimator and whether particular observations have an extreme influence on the estimate through large weights, resulting in large variances. The telescope_match function returns each of the constituent matching counts along with the combined matching “weight” on each observation for use in diagnostic plots.
Bias correction for matching
Matching estimators with a fixed number of matches exhibit bias even in large samples due to differences in the regression function between units and their matches evaluated at their respective covariate values Abadie and Imbens (2011). While this bias converges in probability to \(0\) as the sample size grows, the rate of convergence is slow enough that the bias terms dominate the distribution To address this, matching estimators for treatment effects typically incorporate a bias correction which adjusts the matched values by the estimated difference in regression functions. We incorporate a similar approach for correcting matching bias in our two-stage procedure.
The bias of the simple matching estimator for \(\hat{\tau}\) consists of two terms, a bias due to matching on the mediator (\(B_L^m\)) and a bias due to matching on treatment (\(B_L^a\))
\[ B^m_L = \frac{1}{N} \sum_{i=1}^N (2A_i - 1)M_i \left(1 + \frac{K^a_L(i)}{L}\right) \left( \frac{1}{L} \sum_{{\ell} \in J^m_L(i)} \mu_{A_i0}(X_{\ell}, Z_{\ell}, A_i) - \mu_{A_i0}(X_i, Z_i, A_i) \right) \]
\[ B^a_L = \frac{1}{N} \sum_{i=1}^N (2A_i - 1) \Bigg[\frac{1}{L} \sum_{j \in J^a_L(i)} \mu_{1-A_i,0}(X_i,1-A_i) - \mu_{1-A_i,0}(X_j,1-A_i) \Bigg] \]
Where \(\mu_{am}(x,z,a) = E[Y(a,m) \mid X_i = x, Z_i = z, A_i = a]\) and \(\mu_{am}(x,a) = E[Y_i(a,m) \mid X_i= x, A_i = a]\) denote the conditional expectations of the potential outcomes given two different conditioning sets (with and without \(Z\)). Under sequential ignorability, \(\mu_{am}(x,z,a) = \mu(x,z,a,m) = E[Y_i \mid X_i= x, Z_i = z, A_i = a, M_i = m]\). The method implemented in telescope_match extends the bias correction strategy of Abadie and Imbens (2011) to the two-stage setting. It estimates the two conditional expectation functions using regression estimators \(\widehat{\mu}(x, z, a, m)\) and \(\widehat{\mu}_{a0}(x, a)\). As shown in Blackwell and Strezhnev (2021), if the regression estimators are consistent for the true regression functions, then the estimated bias correction converges to the true bias. The rate of convergence is fast enough such that the bias can be ignored asymptotically.
Inference
Obtaining valid standard errors in the matching context is difficult as matching with replacement induces dependencies between imputed potential outcomes. We provide two approaches for estimating standard errors. The first implements a version of the Otsu and Rai (2017) wild bootstrap for matching estimators, extended to the two-stage setting. The second (default) approach estimates the components of the asymptotic variance derived in Blackwell and Strezhnev (2021).
Empirical illustration of telescope matching
In this section, we illustrate the implementation of the telescope matching estimator as applied to the Job Corps experiment data used in Huber (2014) to estimate the effect of a randomized job training program on self-assessed health holding constant the intermediate outcome of employment.
The data is supplied along with the package.
data(jobcorps)-
\(Y_i\) is an indicator for whether participants reported “very good” health after 2.5 years out after randomization (
exhealth30) -
\(A_i\) is an indicator for assignment to the job training program (
treat) -
\(M_i\) is an indicator for employment 1 to 1.5 years after assignment (
work2year2q) -
\(Z_i\) are the post-treatment, pre-mediator intermediate confounders (
emplq4,emplq4full) -
\(X_i\) are the pre-treatment characteristics (
chobef,trainyrbef,jobeverbef)
The original Huber (2014) paper looks at separate controlled effects for female and male participants. We start by subsetting the dataset down to the female participants
jobcorps_female <- subset(jobcorps, female == 1)We define the two formula objects used for matching, the first including all pre- and post-treatment covariates \(X_i\) and \(Z_i\) along with all treatment-covariate interactions. The second includes only the pre-treatment covariates and treatment-covariate interactions. Here, we include only a subset of the covariates used by Huber (2014) in their analysis.
## Telescope matching formula - First stage (X and Z)
tm_form <- exhealth30 ~ schobef + trainyrbef + jobeverbef | treat |
emplq4 + emplq4full | work2year2qThe telescope_match() function can handle additional mediators and intermediate covariates by simply adding them to the end of the formula in the same manner as these two groups. We then pass this formula to the function itself:
### Estimate ACDE for women holding employment at 0
tm_out <- telescope_match(tm_form, data = jobcorps_female, L = 3)## Beginning matching...## Matching work2year2q...## Matching treat...## Beginning bias correction...The summary() function will print the output (estimate and standard errors) to the console for each of the possible controlled direct effects. This function also provides some summaries of the matching output. The elements of the output can be accessed directly from the returned tmatch object.
# Prints the summary output
summary(tm_out)## Telescope matching results
## ----------------------------
##
## Call:
## telescope_match(formula = tm_form, data = jobcorps_female, L = 3)
##
## Active treatment: treat
## Controlled treatment(s): work2year2q
##
## Matching summary:
## Term Matching Ratio L:1 N == 1 N == 0 Matched == 1 Matched == 0
## 1 treat 3 2801 1551 2275 1545
## 2 work2year2q 3 2559 1793 1369 1514
##
##
## Summary of units matching contributions:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## treat 0 0.3333333 0.6666667 1 1.333333 4.666667
## treat:work2year2q 0 0.0000000 0.3333333 1 1.000000 124.333333
## work2year2q 0 0.0000000 0.3333333 1 1.333333 67.666667
##
##
## Estimated controlled direct effects of treat:
## work2year2q Estimate Estimate (no BC) Std. Err.
## (1, 0) vs. (0, 0) 0 -0.08189 -0.07207 0.05252
## (1, 1) vs. (0, 1) 1 0.03563 0.03562 0.02255
# The coefficients + SE can be accessed directly as well
tm_out$tau # Point estimate## [1] -0.08189419 0.03562810
tm_out$tau_se #Standard error## [1] 0.05251830 0.02255062Additional diagnostics can be conducted by calling the balance.tmatch() function on the returned tmatch object to assess the change in pre-/post-matching covariate balance in both the first and second stages.
## Assess mediator balance on selected pre- and post- treatment covariates
balance.tmatch(tm_out, vars = work2year2q ~ schobef + emplq4,
data = jobcorps_female)## variable before_0 before_1 after_0 after_1 before_sd before_diff
## schobef schobef 0.6352482 0.6287612 0.7040441 0.6314338 0.4824714 -0.006486953
## emplq4 emplq4 0.1996654 0.7198124 0.4778646 0.5057445 0.5000270 0.520147061
## before_std_diff after_diff after_std_diff
## schobef -0.01344526 -0.07261029 -0.15049658
## emplq4 1.04023787 0.02787990 0.05575679
## Assess treatment balance on selected pre-treatment covariates
balance.tmatch(tm_out, vars = treat ~ trainyrbef + hhsize,
data = jobcorps_female)## variable before_0 before_1 after_0 after_1 before_sd
## trainyrbef trainyrbef 0.01676338 0.01713674 0.01700368 0.01700368 0.1292996
## hhsize hhsize 4.54094133 4.53802213 4.53844975 4.54373468 2.2087850
## before_diff before_std_diff after_diff after_std_diff
## trainyrbef 0.0003733584 0.002887544 0.000000000 0.000000000
## hhsize -0.0029191932 -0.001321629 0.005284926 0.002392685Calling the plotDiag.tmatch() function will return a histogram of the number of times each unit is used as a match (\(K/L\)) in either the first (mediator) stage or the second (treatment) stage.
# Number of times each unit is used as a match in the mediator (first) stage
plotDiag.tmatch(tm_out, stage = "work2year2q")
# Number of times each unit is used as a match in the treatment (second) stage
plotDiag.tmatch(tm_out, stage = "treat")
References
Abadie, Alberto, and Guido W Imbens. 2011. “Bias-Corrected Matching Estimators for Average Treatment Effects.” Journal of Business & Economic Statistics 29 (1): 1–11.
Blackwell, Matthew, and Anton Strezhnev. 2021. “Telescope Matching for Reducing Model Dependence in the Estimation of the Effects of Time-Varying Treatments: An Application to Negative Advertising.” Journal of the Royal Statistical Society (Series A).
Huber, Martin. 2014. “Identifying Causal Mechanisms (Primarily) Based on Inverse Probability Weighting.” Journal of Applied Econometrics 29 (6): 920–43.
Otsu, Taisuke, and Yoshiyasu Rai. 2017. “Bootstrap Inference of Matching Estimators for Average Treatment Effects.” Journal of the American Statistical Association 112 (520): 1720–32.