Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,552 characters · 8 sections · 88 citation commands
\setcounter{page}{1}
$^{\dag}$ { [email removed]}, $^{*}$ { [email removed]}, $^{\ddag}$ { [email removed]}, $^{\S}$ { [email removed]. }
{ Acknowledgments: Yu-Chin Hsu gratefully acknowledges research support from the National Science and Technology Council of Taiwan (NSTC111-2628-H-001-001), the Academia Sinica Investigator Award of the Academia Sinica, Taiwan (AS-IA-110-H01), and the Center for Research in Econometric Theory and Applications (107L9002) from the Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education of Taiwan. Chu-An Liu gratefully acknowledges research support from the Academia Sinica Career Development Award (AS-CDA-110-H02).}
\thispagestyle{empty}
While most treatment evaluations focus on binary interventions, a growing literature also considers continuously distributed treatments. We propose a Cram\'{e}r-von Mises-type test for testing whether the mean potential outcome given a specific treatment has a weakly monotonic relationship with the treatment dose under a weak unconfoundedness assumption. In a nonseparable structural model, applying our method amounts to testing monotonicity of the average structural function in the continuous treatment of interest. To flexibly control for a possibly high-dimensional set of covariates in our testing approach, we propose a double debiased machine learning estimator that accounts for covariates in a data-driven way. We show that the proposed test controls asymptotic size and is consistent against any fixed alternative. These theoretical findings are supported by the Monte-Carlo simulations. As an empirical illustration, we apply our test to the Job Corps study and reject a weakly negative relationship between the treatment (hours in academic and vocational training) and labor market performance among relatively low treatment values.
{\bf JEL classification:} C01, C12, C21 \\
{\bf Keywords:} Average dose response functions, average structural function, continuous treatment models, doubly robust, high dimension, hypothesis testing, machine learning, treatment monotonicity.
\thispagestyle{empty} \setcounter{page}{1}
Even though many studies on treatment or policy evaluation investigate the effects of binary or discrete interventions, a growing literature also considers the assessment of continuously distributed treatments, e.g.\ hours spent in a training program whose effect on labor market performance is of interest. Most contributions like Imbens2000, HiranoImbens2004, Flores2007, Floresetal2012, GalvaoWang2015, Lee2018 and CL focus on the identification and estimation of the average dose-response function (ADF), which corresponds to the mean potential outcome as a function of the treatment dose. This permits assessing the average treatment effect (ATE) as the difference in the ADF assessed at two distinct treatment doses of interest, while HiranoImbens2004, Floresetal2012, and CL also consider the marginal effect of slightly increasing the treatment dose, which is the derivative of the ADF. Rather than considering the total effect of the treatment, Huberetal2020 suggest a causal mediation approach to disentangle the ATE into its direct effect and indirect effect operating through intermediate variables or mediators to assess the causal mechanisms of the treatment.
In this paper, we propose a method for testing whether the ADF has a weakly monotonic relationship with (i.e.\ is weakly increasing or decreasing in) the treatment dose under a weak unconfoundedness assumption, implying that confounder of the treatment-outcome relation can be controlled for by observed covariates. Such a test appears interesting for verifying shape restrictions, e.g.\ whether increasing the treatment dose always has a non-negative effect, no matter what the baseline level of treatment is. Moreover, the treatment effect model is known to be equivalent to a nonseparable structural model of a nonseparable outcome with a general disturbance, as for instance IN09ETA and Lee2018. In this case, the ADF corresponds to the average structural function in BP03. Therefore our test can be applied to testing monotonicity of the average structural function in a nonseparable structural model under a conditional independence assumption.
To construct our test, we first transform the null hypothesis of a monotonic relationship to countably many moment inequalities based on the generalized instrumental function approach of HsuLiuShi2019 and HsuShen2020. We construct a Cram\'{e}r-von Mises-type test statistic based on the estimated moments, which are shown to converge to a Gaussian process at the parametric regular root-$n$ rate. Importantly, by making use of moment inequalities, our method does not rely on the nonparametric estimation of the ADF or the marginal effects, which would converge at slower nonparametric rates. To compute the critical value for our test, we apply a multiplier bootstrap method and the generalized moment selection (GMS) approach of AndrewsShi2013, AndrewsShi2014. We demonstrate that our test controls asymptotic size and is consistent against any fixed alternative.
To employ nonparametric or machine learning estimators in the presence of possibly high-dimensional nuisance parameters, we propose a double debiased machine learning (DML) estimator. Utilizing a doubly robust moment function based on a Neyman-type orthogonal score and cross-fitting, we give high-level conditions under which the nuisance estimators do not affect the first-order large sample distribution of the DML estimators. Specifically, we give the high-level conditions on the mean-squared convergence rates on the first-step estimators, as for the semiparametric models in CCDDHNR. The nuisance estimators for the conditional expectation function and the conditional density can be kernel and series estimators, as well as modern ML methods, such as lasso and deep neural networks. See CCDDHNR and AtheyImbens for potential ML methods, such as ridge, boosted trees, and various ensembles of these methods. As each ML method has its strength and weakness depending on the data generating process and applications, it is desired to flexibly employ various nuisance estimators. High-dimensional control variables are accommodated via the nuisance estimators; for example, lasso allows the dimension of $X$ to grow with the sample size.
Our paper is related to a growing literature on testing monotonicity in regression problems such as Bowmanetal1998, Ghosaletal2000, Gijbelsetal2000, HallHeckman2000, DumbgenSpokoiny2001, Durot2003, Beraudetal2005, WangMeyer2011, Chetverikov2019 and HsuLiuShi2019. The main difference between our GMS method and the previously suggested tests is that we rely on a two-step estimation procedure when computing the moments, with the first step consisting of estimating the generalized propensity score, i.e.\ the conditional density of a treatment dose given the covariates, and/or the conditional mean function. For this reason, it is necessary to take into account the behavior of the first step when we derive the limiting behavior of the estimated moment inequalities underlying our test.
We investigate the finite sample behavior of the proposed test approach in a simulation study and also extend our method to testing conditional monotonicity given observed covariates. As an empirical illustration, we apply our test to data from an experimental study on Job Corps, see ScBuGl01 and ScBuMc2008, a program aimed at increasing the human capital of youths from disadvantaged backgrounds in the U.S. We consider hours in academic and vocational training in the first year of the program as the continuous treatment and investigate its association with several labor market outcomes: weekly earnings in the fourth year, earnings and hours worked per week in quarter 16, and a binary employment indicator four years after assignment. For all outcomes, our test clearly rejects weakly negative monotonicity in the treatment when considering treatment doses between 40 and 3000 hours of training. In contrast, weakly positive monotonicity is not refuted at conventional levels of statistical significance. When, however, splitting the treatment range into 3 brackets of 40 to 1000, 1000 to 2000, and 2000 to 3000 hours, the test points to a violation of weakly negative monotonicity only in the lowest treatment bracket. In the remaining brackets with larger treatment values, we neither reject weakly positive, nor weakly negative monotonicity. Our results are consistent with a concave ADF as for instance found in Floresetal2012, suggesting that the marginal effect of training on labor market performance is positive for relatively low treatment doses but decreases as hours in training increase. A potential explanation could be that participants attending more training in the first year might be induced to attain more education in the following years rather than to participate in the labor market.
The paper is organized as follows. Section (ref) formulates the hypothesis of weak monotonicity to be tested. Section (ref) propose monotonicity tests under DML estimation. Section (ref) presents a Monte-Carlo simulation and discusses how to choose the tuning parameters of the test in practice. Section (ref) provides an empirical application to the Job Corps data. Section (ref) adapts the method to testing monotonicity with conditional (rather than unconditional) mean potential outcomes given observed covariates. Section (ref) concludes. The technical proofs are relegated to the Appendix. An online supplement contains monotonicity tests under nonparametric and parametric estimations of generalized treatment propensity score.
Let $Y(t)$ denote the potential outcome corresponding to the level of treatment intensity $t\in\mathcal{T}$, where $\mathcal{T}=[a,b]$ with $-\infty< a<b <\infty$. $Y(t)$ is called the unit-level dose-response function in HiranoImbens2004. Let $\mu(t)=E[Y(t)]$ for $t\in\mathcal{T}$ denote the average of the potential outcome function, also known as the average dose-response function or the average structural function. In this paper, we are interested in testing if the average dose-response function is weakly increasing in the treatment intensity within a specific range. We define the null hypothesis of our interest as
where $a\leq t_\ell<t_u\leq b$ so that $[t_\ell,t_u]$ is a convex and compact subset of $[a,b]$. Without loss of generality, we assume that $[t_\ell,t_u]=[0,1]$.\footnote{If $[t_\ell,t_u]$ is not $[0,1]$, we can always apply an affine transformation $\phi$ on $t$ so that $\phi(t_\ell)=0$ and $\phi(t_u)=1$. }
Note that the null hypothesis in ((ref)) has a form that is similar to that in the literature on regression monotonicity, see for instance HsuLiuShi2019. However, the identification of $\mu(t)$ in our case is different from theirs. We apply the generalized instrumental function approach of HsuLiuShi2019 and HsuShen2020 to transform $H_0$ in ((ref)) to countably many moment inequalities without loss of information.\footnote{The generalized instrumental function approach is a generalization of the instrumental function approach in AndrewsShi2013, AndrewsShi2014.} To be specific, suppose that $\mu(t)$ is a continuous function on $t=[0,1]$ and $h(t)$ is a positive weighting function such that $\int_{0}^1h(t)dt<\infty$. Then by Lemma 2.1 of HsuShen2020, $H_0$ in ((ref)) is equivalent to
for any $q=2,\cdots,$ and for any $t_1\geq t_2$ such that $q\cdot t_1,q\cdot t_2\in \{0,1,2,\cdots,q-1\}$. Equations ((ref)) and ((ref)) hold by the fact that if a function is non-decreasing, then its weighted average over an interval will be non-decreasing as well when the interval moves to the right. In addition, by HsuLiuShi2019, Equations ((ref)) and ((ref)) contain the same information as the null hypothesis.
In the following, we discuss the identification of $\int_{t}^{t+q^{-1}} \mu(s)\cdot h(s) ds$.
Assumption (ref) is a commonly invoked identifying assumption based on observational data, also known as conditional independence and selection on observables. It assumes that conditional on observables $X$, $T$ is as good as randomly assigned, or conditionally exogenous. The observed outcome $Y$ satisfies that $Y=Y(T)$. We then have the following lemma concerning the identification of $\int_{t}^{t+q^{-1}} \mu(s)\cdot h(s) ds$. Let $p(t,x)=f_{T|X}(t|x)$ be the generalized propensity score, which is the conditional density of the treatment given the covariates and $p(t,x)>0$ for all $t$ and $x$.
We now apply Lemma 2.1 of HsuShen2020 and the identification result in Lemma 2.1 to transform $H_0$ in ((ref)) to countably many moment inequalities based on which we will construct our test. For $\ell=(t_1,t_2,q^{-1})\in[0,1]^2\times (0,1]$, define
For each $\ell$, we define
The proof of Lemma (ref) is a direct implication of ((ref)) and Lemma (ref). To see this, set $h(t)=1$ and note that $\int_t^{t+r} h(s) ds=r$. By ((ref)) and Lemma (ref), for any $\ell\in \mathcal{L}$,
In Lemma (ref), we pick $h(t)=1$ for simplicity, but the result also holds for any other known valid wight function $h(t)$.
To deliver a reliable distributional approximation in practice, the double debiased ML (DML) method contains two key ingredients: a doubly robust moment function and cross-fitting. The doubly robust moment function reduces sensitivity in estimating $\nu(\ell)$ with respect to nuisance parameters.\footnote{ Our estimator is doubly robust in the sense that it consistently estimates $\nu(\ell)$ if either one of the nuisance functions $E[Y|T,X]$ or $f_{T|X}$ is misspecified. The rapidly growing ML literature has utilized this doubly robust property to reduce regularization and modeling biases in estimating the nuisance parameters by ML or nonparametric methods; for example, BCH14RES, Farrell15, BCFH17, FLM, CEINR, CCDDHNR, RotheFirpo, and references therein. } Cross-fitting removes bias induced by overfitting and achieves stochastic equicontinuity without strong entropy conditions. Our work builds on the results for semiparametric models in IchimuraNewey22QE, CEINR, CCDDHNR, and the nonparametric models for continuous treatments in CL.
We construct the moment function for our DML estimator by the Gateaux derivative limit. Denote as $\nu(t,r) = \int_t^{t+r} \mu(s) ds$ and $\gamma(t,x) = E[Y|T=t, X=x]$. Let $f^0$ be the true pdf of $Z = (Y,T,X)$ and $f_Z^h$ be a pdf approaching a point mass at $Z $ as $h \rightarrow 0$. CL derive the Gateaux derivative of $\mu(t)$ with respect to a deviation from the true distribution $f_Z^h - f^0$ to be
Since $\nu(t,r)$ is a linear functional of $\mu(t)$, the Gateaux derivative limit of $\nu(t,r)$ is
and it follows that
We propose a DML estimator for $\nu(\ell)$ based on ((ref)):
We use $\|\cdot\|_2$ to denote the $L_2$-norm, e.g. $\|\hat\gamma_k-\gamma\|_2 = \left(\int_\mathcal{X}\int_\mathcal{T}\left( \hat\gamma_k(t,x) - \gamma(t,x) \right)^2 f_{TX}(t, x) dtdx\right)^{1/2}$ and $\|\hat p_k-p\|_2 = \left(\int_\mathcal{X}\int_\mathcal{T}\left( \hat p_k(t,x) - p(t,x) \right)^2 f_{TX}(t, x) dtdx\right)^{1/2}$.
Assumptions (ref)(i) and (ii) are the typical conditions on the mean-squared convergence rates, as in CCDDHNR. Assumption (ref)(iii) is to control the approximation error of the numerical integration.
Lemma (ref) establishes the limiting behavior of DML estimators for $\nu$'s. Let $\hat{\sigma}^2_{\nu,DML}(\ell)={K}^{-1}\sum_{k=1}^K n^{-1}_k\sum_{i \in I_k} \hat{\phi}^2_{\ell,DML}(Y_i,T_i,X_i)$ where
and $\hat{\sigma}^2_{\nu,DML}(\ell)$ will be a consistent estimator for the asymptotic variance of $\sqrt{n}(\hat{\nu}_{DML}(\ell)-{\nu}(\ell))$ under the assumptions Lemma (ref). Let $\hat{\sigma}_{\nu,\epsilon,DML}(\ell)=\max\{\hat{\sigma}_{\nu,DML}(\ell), \epsilon\cdot \hat{\sigma}_{\nu,DML}(0,1/2,1/2)\}$, by which we manually bound the variance estimator away from zero. To test the null hypothesis $H_0'$, we make use of a Cram\'{e}r-von Mises test statistic defined as
where $Q$ is a weighting function such that $Q(\ell)>0$ for all $\ell \in \mathcal{L}$ and $\sum_{\ell\in\mathcal{L}}Q(\ell)<\infty$.
We next define the simulated critical value for our test. We first introduce a multiplier bootstrap method that can simulate a process that converges to the same limit as $\sqrt{n}(\hat{\nu}_{DML}(\ell)-{\nu}(\ell))$. Let $\{U_i:~ 1\leq i\leq n\}$ be a sequence of i.i.d.\ random variables that satisfy Assumption (ref). We construct the simulated process as
where $\hat{\phi}_{\ell,np}(Y_i,T_i,X_i)$ is the estimated influence function defined in ((ref)). Under specific regularity conditions, we can show that the simulated process weakly converges to a Gaussian process conditional on the sample path with probability approaching one and that this limiting Gaussian process corresponds to the limiting process of $\sqrt{n}(\hat{\nu}_{np}(\ell)-\nu(\ell))$.
We adopt the GMS method to construct the simulated critical value as
in which $a_n$ and $B_n$ satisfy Assumption (ref).\footnote{The GMS approach is similar to the recentering method of Hansen2005 and DonaldHsu2016, and the contact approach of Lintonetal2010.}
The decision rule is then given by
The high-level conditions in Assumption (ref) are attainable by various estimators, in particular, kernel, series, deep neural networks, and lasso. The theory of the conventional nonparametric kernel and series methods is well established. Recently FLM provide $\|\hat\gamma_k - \gamma\|_2$ of deep neural networks. CL propose GPS estimators that utilize generic estimators of the conditional mean function. Specifically, Lemmas 1 and 2 in CL provide the convergence rates for their GPS estimators using the deep neural networks in FLM, i.e. $\|\hat p_k-p\|_2$.\footnote{ The convergence rates in Lammas 1 and 2 of CL can be shown to hold uniformly over $t\in\mathcal{T}$, so we can obtain $\|\hat p_k-p\|_2$. The additional assumption for the MultiGPS estimator in Lemma 2 CL is $\sup_{t\in\mathcal{T}}\left\| \hat\mu\left(h_1^{d_t}g_{h_1}(T-t); X\right) - {\mathbb{E}}[h_1^{d_t}g_{h_1}(T-t)|X] \right\|_{F_X} = O_p(R_{1n})$. Then Lemma 3 in CL provides the sufficient conditions. } So Assumptions (ref)(i) and (ii) are attainable by the deep neural networks in FLM and CL. In the rest of this section, we provide the sufficient low-level conditions for lasso methods.
We illustrate how to employ lasso methods to estimate the nuisance conditional mean function $\gamma(t,x)$ and the generalized propensity score $f_{T|X}(t|x)$. We provide sufficient conditions to verify the high-level Assumption (ref). We modify the penalized local least squares estimator of $\gamma(t, X)$ in SUZ (SUZ, hereafter). We use the conditional density estimator in SUZ. For completeness, we present the estimators and asymptotic theory in SUZ and refer readers to SUZ for details.
Let $b(T, X)$ be a $p\times 1$ vector of basis functions. We approximate $\gamma(t,x)$ by $b(t,x)'\theta$. The lasso estimator $\hat\gamma_k(t,x) = b(t,x)'\hat\theta_k$ for $k\in\{1,...,K\}$, where
where $n_{k} = \sum_{i=1}^n {\bf 1}\{i \in I_k\}$, $\|\cdot\|_1$ denotes the $L_1$ norm, $\lambda = \ell_n(\log(p\vee n)n)^{1/2}$ for some slowly diverging sequence $\ell_n$, and $\hat\Xi_k = diag(\tilde l_{k1},...,\tilde l_{kp})$ is a generic penalty loading matrix computed by Algorithm (ref) below from the iterative Algorithm 3.1 in SUZ. Denote as $\| f(X) \|_{\mathbb{P}_{nk},2} = \big((n-n_k)^{-1}\sum_{i \notin I_k} f(X_i)^2\big)^{-1/2}$ for a generic function $f(\cdot)$.
Let the final penalty loading matrix $\hat\Xi_k=\hat\Xi^S_k$ from Algorithm (ref). Then the lasso estimator $\hat\gamma_k(t,x) = b(t,x)'\hat\theta_k$ for $k\in\{1,...,K\}$ from ((ref)).
To estimate the conditional density $p(t,x) = f_{T|X}(t|x)$, first estimate the conditional CDF $F_{T|X}$ by the logistic distributional lasso regression and then take the numerical derivative. Let $b(X)$ be a $p\times 1$ vector of basis functions. We approximate $F_{T|X}(t|x)$ by $\Lambda(b(x)'\beta_t)$, where $\Lambda$ is the logistic CDF. For $k\in\{1,...,K\}$, $\hat F_{T|X_k}(t|x) = \Lambda(b(X)'\hat\beta_{tk})$, where
where $n_{k} = \sum_{i=1}^n {\bf 1}\{i \in I_k\}$, $M(y, x; g) = -[y\log(\Lambda(b(x)'g)) + (1-y)\log(1-\Lambda(b(x)'g))]$ is the logistic likelihood, the penalty $\tilde\lambda = 1.1\Phi^{-1}(1-r/\{p\vee nh_1\}) n^{1/2}$, for some $r\rightarrow 0$ and $h_1\rightarrow 0$, and $\Phi$ is the standard normal CDF. A generic penalty loading matrix $\hat\Psi_{tk}$ is computed by Algorithm (ref) below from the iterative Algorithm 3.2 in SUZ.
Let the final penalty loading matrix $\hat \Psi_{tk} = \hat\Psi_{tk}^S$ from Algorithm (ref). Compute $\hat F_{T|X_k}(t|x) = \Lambda(b(X)'\hat\beta_{tk})$ from ((ref)). Then the conditional density estimator \[ \hat p_k(t,x) = \frac{\hat F_{T|X_k}(t+h_1|x) - \hat F_{T|X_k}(t-h_1|x)}{2h_1}. \]
Assumption (ref) collects the conditions in Theorems 3.1 and 3.2 in SUZ. Following SUZ's notations, denote as $\|\cdot\|_{\mathbb{Q},q}$ the $L^q$ norm under measure $\mathcal{Q}$ and $\mathbb{P}$ assigns probability $1/n$ to each observation.
Let Assumption (ref) hold. Then Theorems 3.1 and 3.2 in SUZ imply that $\sup_{(t, x)\in\mathcal{T}\times\mathcal{X}}|\hat\gamma_k(t, x) - \gamma(t, x) |= O_p(A_n)$, where $A_n = \ell_n(\log(p\vee n) s^2\zeta_n^2/n)^{1/2}$ and $\sup_{(t, x)\in\mathcal{T}\times\mathcal{X}}|\hat p_k(t, x) - p(t, x)| = O_p(B_n)$, where $B_n = h_1^{-1}(\log(p\vee n) s^2\zeta_n^2/n)^{1/2}$. Then we can obtain the same rates for the root-mean-squared rates $\|\hat\gamma_k - \gamma\|_2 $ and $\|\hat p_k - p\|_2$ to verify Assumption (ref). Therefore a sufficient condition of Assumption (ref)(i) is $A_n \rightarrow 0$ and $B_n \rightarrow 0$. And a sufficient condition of Assumption (ref)(ii) is $\sqrt{n}A_n B_n \rightarrow 0$.
This section provides a simulation study to examine the finite sample performance of the proposed test. To implement our test in practice, one has to choose several tuning parameters in advance. We make the following propositions concerning the choice of these parameters and present related Monte Carlo simulation results further below.
For all data generating processes (DGPs), the continuous treatment variable $T$, the control variables $X$, and the error term $U_y$ are generated as follows
where the $(i,j)$-entry $\Sigma_{ij}=(0.5)^{|i-j|}$ for $i,j=1,\dots,100$, $U_t\sim\mathcal{N}(0,1)$, and $U_y$, $U_t$, and $X$ are mutually independent. We set $\beta_j=1/j^2$ for mild dependence between $X_j$ and $\beta_j=1/j$ for strong dependence between $X_j$. Three cases of the potential outcomes are studied:
In DGP 1, $\mu(t)=0$, and $H_0$ holds with moment equalities. In this case, we expect that the size of the proposed test will achieve the nominal level since every moment would hold with equality. In DGP 2, $\mu(t)= t^2$, and $H_0$ holds with strict moment inequalities. In this case, we expect the size will converge to zero since every moment would hold with strict inequality. This is because the test statistics will converge to zero and the critical value is bounded away from zero. In DGP 3, $\mu(t)=\sin(\pi T)$, and $H_0$ does not hold. In this case, we expect the power will increase with the sample size.
In these DGPs, $1+d_x=101$. We consider samples of sizes $n = 200$, $400$, $800$, and $1600$. For $q_1$, we set $q_1=4$ for $n=200$, $q_1=8$ for $n=400$, $q_1=16$ for $n=800$, and $q_1=32$ for $n=1600$. The number of subsamples used for cross-fitting is $K\in\{2,5,10\}$. All our Monte Carlo results are based on $1000$ simulations. In each simulation, the critical value is approximated by $1000$ bootstrap replications. The nominal size of the test is set at $10\%$.
To estimate the conditional mean function $\gamma(t,x) = E[Y|T=t, X=x]$, we employ the lasso regression, where the penalization parameter is chosen via grid search utilizing 10-fold cross validation. To estimate the conditional density estimation $p(t,X)$, we first estimate $F_{T|X}(t|x)$ by the logistic distributional lasso regression, and then take the numerical derivative. The penalization parameter of the distributional lasso regression is estimated by Algorithm 3.2 of SUZ. Also, all lasso estimations include an intercept and the covariates. For numerical integration in Step 2, we set $M=[n^{2/3}]$, where $[\cdot]$ is the nearest integer. Our test is based on the trimmed generalized propensity score estimator, defined as $\tilde{p}(T_i, X_i)=\max\{\hat{p}(T_i, X_i),0.025\}$, implying that conditional treatment densities below 2.5% are set to 2.5%.\footnote{In general, one can follow DHL2014 and HLL2020 and trim the estimated generalized propensity scores to prevent them from being too close zero, in order to obtain a more stable IPW estimator whose variance is not affected by extremely low scores.}$^,$\footnote{Based on this trimming rule, around 0.5% of the samples are trimmed.}
Table (ref) shows the rejection probabilities of our test for DGPs 1-3, and the results are consistent with our theoretical findings. For the mild dependence case, the proposed test controls size well in DGP 1 and DGP 2, and the rejection probabilities increase with the sample size and are greater than the nominal size $0.1$ in DGP 3. For the strong dependence case, our test still control size will in both DGP 1 and DGP 2. The power increases with the sample size in DGP 3, but the rejection probabilities are a bit less than the nominal size $0.1$ for $n=200$. Overall, we do not find significant difference for different choices of $K$.
We next investigate the robustness of the performance of our test to the choice of $q_{1}$. Let $N$ denote the expected sample size of the smallest cube. We consider three alternative choices of $q_{1}$, each resulting in $N = 33$, $40$, and $66$, respectively. Table (ref) shows the rejection probabilities of our test for different choices of $N$. The results suggest that the choice of $q_{1}$ does not affect the test performance much. Therefore, the finite sample behavior of our test appears to be reasonably robust to different values of $q_{1}$.
As an empirical illustration, we apply our test to data from the Job Corps study. The latter was conducted between November 1994 and February 1996 to evaluate the publicly funded U.S.\ Job Corps program and used an experimental design that randomly assigned access to the program. Job Corps targets youths from low-income households who are between 16 and 24 years old and legally reside in the U.S. Program participants obtained on average roughly 1200 hours of vocational and/or academic classroom training as well as housing and board over an average duration of 8 months. We refer to ScBuGl01 and ScBuMc2008 for a detailed discussion of the study design and the average effects of program assignment on a range of different outcomes. Their results suggest that Job Corps raises educational attainment, reduces criminal activity, and increases labor market performance measured by employment and earnings, at least for some years after the program.
Particularly relevant for our context is the study by Floresetal2012, who consider the length of exposure to academic and/or vocational training as continuously distributed treatment to assess its effect on earnings based on regression and weighting estimators (using the inverse of the conditional treatment density as weight). As the length of treatment exposure is (in contrast to Job Corps assignment) not random, they impose a selection-on-observables assumption and control for baseline characteristics at Job Corps assignment. While the authors find overall positive average effects of increasing hours in academic and vocational instruction, the marginal effects appear to decrease with length of exposure, pointing to a potential concavity in the association of earnings and time of instruction. Relatedly, Lee2018 and CL assess the effect of hours in training on the proportion of weeks employed in the second year after program assignment based on kernel regression and double machine learning, respectively. Also for this outcome, the plotted regression lines in either study point to a concave association with the treatment dose.\footnote{See also Huberetal2020, who use a causal mediation approach to assess the direct effect of the treatment dose on the number of arrests in the fourth year after program assignment when controlling for employment behavior in the second year based on inverse probability weighting and find a non-linear association.}
However, in the light of estimation uncertainty, mere eye-balling of the outcome-treatment associations in empirical applications does not tell us whether specific shape restrictions can be refuted. For this reason, we use our DML method with lasso regression for nuisance parameter estimation to formally test whether weak positive and negative monotonicity can be rejected in the Job Corps data when considering several labor market outcomes. To this end, we define the treatment variable $T$ as the total hours spent in academic and vocational training in the 12 months following the program assignment. Our outcomes $Y$ include weekly earnings in the fourth year, earnings and hours worked per week in quarter 16, and a binary employment indicator four years after assignment (i.e.\ in week 208).
For invoking weak unconfoundedness (Assumption 2.1), we consider the same set of pre-treatment covariates $X$ as Lee2018, CL, and Huberetal2020, which overlaps with the control variables of Floresetal2012.\footnote{A control variable in Floresetal2012 we do not have access to is the local unemployment rate which was constructed by matching county-level unemployment rates to individual postal codes of residence, which are only available in a restricted-use data set.} We condition on individual characteristics like age, gender, ethnicity, language competency, education, marital status, household size and income, previous receipt of social aid, family background (e.g. parents' education), criminal activity, as well as health and health-related behavior (e.g.\ smoking, alcohol, or drug consumption). Conditioning on such a rich set of socio-economic variables appears important, as the satisfaction of weak unconfoundedness relies on successfully controlling for all factors jointly affecting treatment duration and labor market behavior. Furthermore, we include variables that might be associated with the duration in training, namely expectations about Job Corps and interaction with the recruiters, which might serve as proxies for unobserved personality traits (like motivation) that could also affect the outcomes. Finally, we control for pre-treatment outcomes, namely previous labor market participation and earnings, to tackle any confounders that affect the outcomes of interest through their respective pre-treatment values.
The original Job Corps data set consists of $15,386$ individuals prior to program assignment, but a substantial share never enrolled in the program and dropped out of the study, such that there are only $11,313$ individuals with completed follow-up interviews four years after randomization. Among those, $6,828$ had been randomized into Job Corps and had thus access to academic or vocational training. To define our final evaluation sample, we follow Floresetal2012, Lee2018, CL, and Huberetal2020 and consider observations with at least 40 hours (or one working week) of training for our analysis, all in all $4,166$ individuals. Among these, there are cases of item non-response in various elements of $X$ measured at the baseline survey, for which we account by the inclusion of missing dummies as additional regressors, while observations with missing values in the outcome of interest need to be dropped when running the respective test. Table (ref) provides descriptive statistics for selected covariates $X$ (see Huberetal2020 for a full list of control variabes) as well as for the treatment $T$ and all outcomes $Y$, including the respective number of nonmissing observations (nonmissing).
The choices of nuisance parameters are the same as in the simulations (see the previous section). The number of subsamples used for cross-fitting is $5$, and the expected sample size of the smallest cube is either $40$ or $50$. The lasso estimations include an intercept, the covariates and the squared terms of any non-binary covariates. The $p$-values of the tests for the various outcomes are calculated based on 1000 bootstrap replications.\footnote{In our empirical study, we do not get unstable IPW $\nu(\ell)$ estimates, so we decide not to apply the trimming method. Also, we note that all estimated generalized propensity scores are greater than $0.0001$ in our empirical study.}
In a first step, we apply the test to a treatment interval of $T \in [40,3000]$, where choosing 3000 hours of training as upper bound of the analysis is motivated by the quickly decreasing number of observations beyond that point.
Table (ref) reports the test statistics and p-values for all outcomes under both null hypotheses of weakly increasing mean potential outcomes in the treatment ($\mu(t_1)\geq\mu(t_2)$ for $t_1>t_2$) and weakly decreasing mean potential outcomes ($\mu(t_1)\leq\mu(t_2)$), respectively. Our tests clearly reject the latter hypothesis of weakly negative monotonicity for any labor market outcome at the 1% level of statistical significance. In contrast, weak positive monotonicity is never rejected, as any test yields p-values close to or equal to 1 (or 100%). Our findings therefore suggest that an increase in the treatment does either increase or at least not reduce the outcome over the treatment range $T \in [40,3000]$.
It is worth mentioning that the concavities in the outcome-treatment associations spotted in the previously mentioned empirical applications suggest decreasing marginal effects when increasing the treatment. In our testing context, this implies that weakly negative monotonicity should be more clearly rejected for lower rather than higher ranges of treatment values by our method. To verify this suspicion, we in a second step partition the treatment support into three sets of $[40,1000]$, $[1000,2000]$, and $[2000,3000]$ and run the tests separately within each set.
Table (ref) presents the results for $T \in [40,1000]$. None of the tests rejects weakly positive monotonicity at any conventional level of significance, while all tests strongly reject weakly negative monotonicity. For the intermediate treatment range of $[1000,2000]$ considered in Table (ref), however, neither positive nor negative monotonicity is ever rejected at the 10% level of statistical significance. This implies that marginal treatment effects are generally less positive than for lower values of $T$. The same findings apply to the highest treatment bracket $[2000,3000]$, where all tests yield p-values which are beyond conventional levels of significance. Summing up, our empirical findings are consistent with a concave mean potential outcome-treatment dependence, implying that initially strongly positive marginal treatment effects decrease as the treatment value considered (hours in training) increases. A potential explanation for the concavity could be that individuals attending more training in the first year might be induced to attain more education also in the following years rather than to participate in the labor market.
In this section, we adapt our method to testing monotonicity with conditional (rather than unconditional) mean potential outcomes given observed covariates $X$. In this case, the null hypothesis considered corresponds to
where $\mu(t,x)=E[Y(t)|X=x]$ is the conditional average of the potential outcome function or the average dose-response function. For simplicity and without loss of generality, we henceforth assume that $X$ is a scalar with $\mathcal{X}=[0,1]$. By Lemma 2.1 of HsuShen2020, $H_0$ in ((ref)) is equivalent to
for any $q=2,\cdots,$ and for any $t_1\geq t_2$ such that $q\cdot t_1,q\cdot t_2, q\cdot x \in \{0,1,2,\cdots,q-1\}$. Define $h(t,x)=f(x)$ to be the density function of $X$. Following Lemma (ref) and ((ref)), we have for $r>0$,
For $\ell_x=(t_1,t_2,x,q^{-1})\in[0,1]^3\times (0,1]$, we let
Similar to ((ref)), for each $\ell_x$, we define
This permits establishing the following lemma.
Similar to Section (ref), we estimate $\nu(\ell_x)$ with $\ell_x=(t_1,t_2,x,q^{-1})$ as the following:
Similar to Lemma (ref), we can show that uniformly over $\ell_x\in\mathcal{L}_x$,
where
Let $\hat{\phi}_{\ell_x,DML}(Y,T,X)$ be the estimated influence function similar to ((ref)) and let $\hat{\sigma}^2_{\nu,DML}(\ell_x)=K^{-1}\sum_{k=1}^K{n^{-1}_k}\sum_{i\in I_k} \hat{\phi}^2_{\ell_x,DML}(Y_i,T_i,X_i)$ which will be a consistent estimator for the asymptotic variance of $\sqrt{n}(\hat{\nu}_{DML}(\ell_x)-{\nu}(\ell_x))$ under proper regularity conditions. Furthermore, let $\hat{\sigma}_{\nu,\epsilon,DML}(\ell_x)=\max\{\hat{\sigma}_{\nu,DML}(\ell_x), \epsilon\cdot \hat{\sigma}_{\nu,DML}(0,1/2, 0,1/2)\}$. The Cram\'{e}r-von Mises test statistic is defined as
where $Q$ is a weighting function such that $Q(\ell_x)>0$ for all $\ell_x \in \mathcal{L}_x$ and $\sum_{\ell_x\in\mathcal{L}_x}Q(\ell_x)<\infty$. The simulated process is constructed as
The GMS simulated critical value is given by
Finally, the decision rule is given by
The size and power properties are similar to the unconditional potential outcome cases and the details are omitted for brevity.
In this paper, we propose Cram\'{e}r-von Mises-type tests for testing whether a mean potential outcome is weakly monotonic in a continuously distributed treatment under a weak unconfoundedness assumption. To flexibly employ nonparametric or machine learning estimators in the presence of possibly high-dimensional nuisance parameters, we propose a double debiased machine learning estimator for the moments entering the test. Furthermore, we extend our method to testing monotonicity conditional on observed covariates. We also investigate the test's finite sample behavior in a simulation study and find it to perform decently under our suggested choices of tuning parameters.
As an empirical illustration, we apply our test to the Job Corps study, investigating the associations of several labor market outcomes (earnings, employment, and hours worked) with hours in training as treatment. We find that an increase in the treatment does either increase or at least not reduce the outcome. When splitting the treatment range into subsets, our testing results are consistent with a concave mean potential outcome-treatment dependence, implying that initially stronger marginal treatment effects decrease as the treatment value (i.e.\ hours already spent in training) increases.