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.
80,082 characters · 13 sections · 80 citation commands
A Unified Framework for Specification Tests of Continuous Treatment Effect Models
\sloppy
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Consistent tests; Continuous treatment effect; Series estimation; Bootstrap.
\spacingset{1.45}
Causal inference is a central topic in economics, statistics, and machine learning. Although a randomized trial is the gold standard for identifying causal effects, such trials are often unavailable or even unethical in practice. Observational data, which are collected when the participation of an intervention is only observed rather than manipulated by scientists, are predominantly the type of data that are available. A major challenge for inferring causality in observational studies is confoundedness, whereby individual characteristics are correlated with both the treatment variable and the outcome of interest. To identify causality, the unconfounded treatment assignment condition is frequently imposed in the literature; see rosenbaum1983central,rosenbaum1984reducing. For a comprehensive review of causal inference and its applications, see imbens2009recent and abadie2018econometric.
Treatment effect models are used extensively in economics and statistics to evaluate the causal effect of a treatment or policy. Most of the existing literature focuses on binary treatment, whereby an individual either does or does not receive the treatment (e.g., hahn1998role,Hirano03,donald2014testing,imai2014covariate,abrevaya2015estimating,chan2016globally,Athey2018Approximate,hsu2020counterfactual,chen2020stochastic,fan2020estimation,sant2018covariate,aisimple). Some studies focus on multivalued treatment (see, e.g., cattaneo2010efficient,lee2018efficient,ai2020mann,ao2021multivalued). However, in many applications, the treatment variable is continuously valued, and its causal effect is of great interest to decision makers. For example, when evaluating how non-labor income affects the labor supply, the causal effect may depend on not only the introduction of the non-labor income but also the total non-labor income. Similarly, when evaluating how advertising affects the campaign contributions for political analysis, the causal effect may depend not only on whether any advertisements are released but also on how many of them are distributed.
Estimation of continuous treatment effects has received considerable attention from researchers (see hirano2004propensity,galvao2015uniformly,kennedy2017non,Fong_Hazlett_Imai_2018,dong2019regression,huber2020direct,colangelo2020double,ai2021estimation, among others). hirano2004propensity, galvao2015uniformly, and Fong_Hazlett_Imai_2018 applied fully parametric methods by modeling either the conditional distribution of the treatment given the confounders or that of the observed outcome given the treatment and the confounders. The shortcoming of these parametric methods is that modeling and testing the relationship between the treatment and the observed outcome regarding the confounders are difficult, especially when multiple confounding variables are involved. If the model is mis-specified, the conclusion can be biased and completely misleading. kennedy2017non and huber2020direct estimated the continuous treatment effects using the nonparametric kernel method. Although nonparametric approaches are much more flexible than parametric ones, the former require smoothing of the data rather than estimating finite dimensional parameters, which leads to less precise fits and slower convergence rates (slower than $N^{-1/2}$). Furthermore, it is usually difficult to interpret nonparametric results.
In a recent article, Ai_Linton_Motegi_Zhang_cts_treat studied continuous treatment effects by imposing a univariate generalized parametric model for the functionals of the potential outcome over the treatment variable. The general framework includes many important causal parameters as special cases, for example, average and quantile treatment effects. They proposed a generalized weighting estimator for the causal effect with the weights modeled nonparametrically and estimated by solving an expanding set of equations. They further derived the semiparametric efficiency bound for the causal effect of treatment under the unconfounded treatment assignment condition and showed that their estimator is $\sqrt{N} $-asymptotically normal and attains the semiparametric efficiency bound. Although Ai_Linton_Motegi_Zhang_cts_treat's estimator enjoys superior asymptotic properties and satisfactory finite sample performance, they did not detail the specifications of the parametric models for the functionals of the potential outcomes. If the parametric model is mis-specified, the results developed in Ai_Linton_Motegi_Zhang_cts_treat do not hold.
We study the question of model specification. In particular, we propose a consistent specification test for the most generalized continuous treatment effect model. That is, we consider the generalized parametric model in Ai_Linton_Motegi_Zhang_cts_treat as the null model while testing our hypothesis. The potential outcome variable in the model is not observable. However, under the unconfounded treatment assignment condition, the model can be identified by a semiparametric weighted conditional model. There is abundant literature on the specification tests for conditional models ( e.g., ait2001goodness,bierens1982consistent,bierens1990consistent,fan1996consistent,zheng1996consistent,bierens1997asymptotic,stute1997nonparametric,li1999consistent,chen1999consistent,fan2000consistent,li2003consistent,crump2008nonparametric ). Most authors have considered the problems of testing a parametric/semiparametric null model using an integrated type test statistic. ait2001goodness and chen1999consistent considered testing nonparametric/semiparametric null models using nonparametric kernel methods. li2003consistent considered testing the nonparametric/semiparametric using series methods. crump2008nonparametric derived a nonparametric Wald test statistic for testing the conditional average treatment effects under the unconfoundedness condition. For the binary treatment effect model, shaikh2009specification parametrically modeled the propensity score as a conditional expectation of the treatment given the confounders and proposed an associated specification test. This differs from our problem in the sense that we consider the specification test for the function of the potential outcome with a continuous treatment variable.
Specifically, we estimate our semiparametric weighted null model using the framework developed in Ai_Linton_Motegi_Zhang_cts_treat and construct a Cramér–von Mises test statistic and a Kolmogorov–Smirnov test one to test the null model. Although the weights in our null model are estimated nonparametrically, we show that our proposed test statistic is more efficient than that constructed from the true weights. Moreover, our proposed test statistic can detect local alternatives that deviate from the null model at the rate of $O(N^{-1/2})$.
Under the null hypothesis our test statistic is shown to converge in distribution to a weighted sum of independent chi-squared random variables. It is known that obtaining the exact critical values of such a distribution is extremely difficult in practice. Most of the literature suggests using a residual wild bootstrap procedure to approximate the critical values. This is not applicable in our case because our null model does not imply any explicit form of relationship among the observed outcome, the treatment, and the confounders for residual sampling. To resolve this problem, we adopt a special case of the exchangeable bootstrap to approximate the null limiting distribution. Monte-Carlo simulations and real data analysis were conducted to demonstrate the numerical properties of our test method and limiting distribution approximation.
The remainder of the paper is organized as follows. We introduce the problem formulation and notations in Section (ref). Section (ref) constructs the test statistic, followed by the study of the asymptotic properties under null hypothesis, the fixed and the local alternatives in Section (ref). In Section (ref), we discuss how to approximate the limiting distribution under the null hypothesis. Finally, Section (ref) discusses the choice of the tuning parameters in the estimation and investigates the finite sample performance through simulations and U.S. campaign advertisement data. All proofs are detailed in the supplementary file.
Let $T$ denote a continuous treatment variable with support $\mathcal{T}\subset\mathbb{R}$, where $\mathcal{T}$ is a continuum subset, and $T$ has a marginal density function $f_{T}(t)$. Let $Y^{\ast}(t)$ denote the potential response when treatment $T=t$ is assigned. We are interested in testing the null hypothesis:
against the alternative hypothesis
where $\Theta$ is a compact set in $\mathbb{R}^{p}$ for some integer $p\geq1$, $m(\cdot)$ is some generalized residual function which could possibly be non-differentiable, and $g(t;\boldsymbol{\theta})$ is a parametric working model which is differentiable with respect to $\boldsymbol{\theta}$. If $H_{0}$ holds, for each $t$, the dose-response function (DRF) is defined as the value $g(t;\boldsymbol{\theta}^{*})$ that solves the moment condition in (ref). The following examples show that the average dose-response function (ADRF) and the quantile dose-response function (QDRF) are special cases of $g(t;\boldsymbol{\theta}^{*})$, which result from choosing specific forms of $m(\cdot)$.
We consider an observational study in which the potential outcome $Y^{\ast }(t)$ is not observed for all $t$. Let $Y:=Y^{\ast }(T)$ denote the observed response. Under the null hypothesis, one may attempt to solve the following equation to find $\boldsymbol{\theta}^*$:
However, if there is a selection into treatment, even under the null hypothesis, the true value $\boldsymbol{\theta }^*$ does not solve the above equation. Indeed, in this case, the observed response and the treatment assignment data alone cannot identify ${\boldsymbol{\theta}}^*$. To address this identification issue, most studies in the literature impose a selection on the observable condition Hirano03, imai2004causal, Fong_Hazlett_Imai_2018,Ai_Linton_Motegi_Zhang_cts_treat. Specifically, let $\boldsymbol{X}\in\mathbb{R}^r$, for some integer $r\geq 1$, denote a vector of observable covariates. The following condition shall be maintained throughout the paper.
Let $\{T_i,\boldsymbol{X}_i,Y_i\}_{i=1}^N$ be an independent and identically distributed ($i.i.d.$) sample drawn from the joint distribution of $(T,\boldsymbol{X},Y)$. Let $f_{T|X}$ denote the conditional density of $T$ given the observed covariates $\boldsymbol{X}$. Under Assumption (ref), Ai_Linton_Motegi_Zhang_cts_treat showed that $\mathbb{E}[m\{Y^*(t);g(t;\boldsymbol{\theta})\}]$ can be identified as follows:
where
The function $\pi_{0}(T,\boldsymbol{X})$ is called the stabilized weights in robins2000marginal.
The null and alternative hypothesis in (ref) can then be re-written as
against the alternative hypothesis
This converts the test for (ref) to a specification test for a univariate regression model, if both $\pi_0(T,\boldsymbol{X})$ and $\boldsymbol{\theta}^*$ were given. Specially, letting
the null hypothesis $H_0$ is equivalent to $\mathbb{P}\{\mathbb{E}(U_i|T_i)=0\}=1$. A popular technique for testing such a conditional moment model is to convert it to an unconditional one.
Note that $\mathbb{P}\{\mathbb{E}(U_i|T_i)=0\}=1$ if and only if $\mathbb{E}\{U_iM(T_i)\}=0$ for all bounded and measurable functions $M(\cdot)$. Following bierens1997asymptotic, stinchcombe1998consistent, stute1997nonparametric, and li2003consistent, by choosing a proper weight function $\mathscr{H}(\cdot,\cdot)$, $\mathbb{E}(U_i|T_i)=0$ is a.s. equivalent to
Popular choices of such a weight function are the logistic function $\mathscr{H}(T_i,t)=1/\{1+\exp(c-t\cdot T_i)\}$ with $c\neq 0$, cosine-sine function $\mathscr{H}(T_i,t)=\cos(t\cdot T_i)+\sin(t\cdot T_i)$ and the indicator function $\mathscr{H}(T_i,t)=\mathds{1}(T_i\leq t)$ (see stinchcombe1998consistent and stute1997nonparametric for more detailed discussion). Now, letting
the sample analogue of $\mathbb{E}\left\{U_i\mathscr{H}(T_i,t)\right\}$ multiplied by $\sqrt{N}$, one can test $H_0$ using the Cramér–von Mises (CM)-type statistic
or the Kolmogorov-Smirnov (KS)-type statistic
where $\widehat{F}_T(\cdot)$ is the empirical distribution of $T_1,...,T_N$. However, both $\pi_0(T,\boldsymbol{X})$ and $\boldsymbol{\theta}^*$ are unknown in practice so that the $U_i$'s are unavailable. We must replace the $U_i$'s with some estimates, which is studied in the following section.
Remark 1. Note that in our model, the stabilized weights $\pi_{0}(T,\boldsymbol{X})$ are nonparametric. The conditional distribution of the treatment, given the confounders $f_{T|X}(T|\boldsymbol{X})$, which is known as the generalized propensity score hirano2004propensity, is also nonparametric. Under the unconfoundedness assumption, an alternative identification of the null model is through the conditional distribution of the outcome given the treatment and the confounders and the marginal distribution of the confounders, that is, $\mathbb{E}[m\{Y^*(t);g(t;\boldsymbol{\theta})\}]=\mathbb{E}\Big(\mathbb{E}[m\{Y;g(T;\boldsymbol{\theta})\}|\boldsymbol{X},T=t]\Big)$ holds under Assumption (ref).
One obvious approach for estimating the $U_i$'s is to estimate $f_{T}(T_i)$ and $f_{T|X}(T_i|\boldsymbol{X}_i)$, then construct the estimators of $\pi_0(T_i,\boldsymbol{X}_i)$ and $\boldsymbol{\theta}^*$. However, it is well-known that this ratio estimator of $\pi_0(T,\boldsymbol{X})$ is very sensitive to small values of $ f_{T|X}(T|\boldsymbol{X})$ because small estimation errors in estimating $ f_{T|X}(T|\boldsymbol{X})$ result in large estimation errors of the estimator of $\pi_0(T,\boldsymbol{X})$. To avoid or mitigate this problem, Ai_Linton_Motegi_Zhang_cts_treat directly estimated $\pi_{0}(T,\boldsymbol{X})$ as a whole using the generalized empirical likelihood (GEL). We adopt their estimator and elaborate its construction as follows. Note that the weighting function satisfies
for any suitable functions $u(t)$ and $v(\boldsymbol{x})$. Ai_Linton_Motegi_Zhang_cts_treat showed that the restriction (ref) identifies the weighting function $\pi_{0}(T,\boldsymbol{X})$. This result suggests that one may estimate the $\pi_0(T_i,\boldsymbol{X}_i)$'s by solving the sample analogue of (ref). The challenge is that ((ref)) implies an infinite number of equations, which is impossible to solve with a finite sample of observations. To overcome this difficulty, Ai_Linton_Motegi_Zhang_cts_treat suggested approximating the infinite-dimensional function space by a sequence of finite-dimensional sieve spaces. Specifically, let $u_{K_{1}}(T)=(u_{K_{1},1}(T), \ldots,u_{K_{1},K_{1}}(T))^{\top }$ and $v_{K_{2}}(\boldsymbol{X} )=\big( v_{K_{2},1}(\boldsymbol{X}),\ldots,$ $v_{K_{2},K_{2}}(\boldsymbol{X} )\big) ^{\top }$ denote some known basis functions with dimensions $ K_{1}\in \mathbb{\ N}$ and $K_{2}\in \mathbb{N}$ respectively, and let $ K:=K_{1}\cdot K_{2}$. The functions $u_{K_{1}}(t)$ and $v_{K_{2}}( \boldsymbol{x})$ are called the approximation sieves, such as B-splines or power series Newey97, chen2007large. Because the sieve approximating space is a subspace of the original function space, $ \pi_{0}(T,\boldsymbol{X})$ also satisfies
Following Ai_Linton_Motegi_Zhang_cts_treat, we estimate the $\pi_0(T_i,\boldsymbol{X}_i)$'s consistently by the $\widehat{\pi}_i$'s that maximize the generalized empirical likelihood (GEL) function, subject to the sample analog of (ref): { \
}Two observations are immediately clear. First, by including a constant of one in the sieve base functions, (ref) guarantees that $N^{-1}\sum_{i=1} ^{N}\widehat{\pi}_{i}=1$. Second, we notice that
The entropy maximization problem minimizes the Kullback-Leibler divergence between the weights $\{N^{-1}\pi_{i}\}_{i=1}^{N}$ and the empirical frequencies $\{N^{-1}\}$, subject to the sample analogue of (ref). Further, Ai_Linton_Motegi_Zhang_cts_treat showed that the dual solution of the primal problem (ref) is
where $\rho^\prime$ is the first derivative of $\rho$ with $\rho(u)=-\exp(-u-1)$, and $\widehat{\Lambda}_{K_{1}\times K_{2}}$ is the maximizer of the strictly concave function $\widehat{G}_{K_{1}\times K_{2}}$ defined by
The first order condition of (ref) implies that $\{\widehat{\pi}_{K}(T_{i},\boldsymbol{X}_{i})\}_{i=1}^N$ satisfies the sample analog of (ref); such restrictions reduce the chance of obtaining extreme weights. The concavity of (ref) enables us to obtain the solution quickly via the Gauss-Newton algorithm. To ensure a consistent estimate of $\pi_0(T,\boldsymbol{X})$, the dimensions of the bases, $K_1$ and $K_2$, shall increase as the sample size increases. The choice of $K_1$ and $K_2$ in practice will be discussed in Section (ref).
Having estimated the weights, we now propose an extremum estimator for $\boldsymbol{\theta}^*$ (e.g., pakes1989simulation, chen2003estimation, de2019smoothed). Note that under $H_0$, the true value $\boldsymbol{\theta}^*$ solves the following equation:
where $w(T;\boldsymbol{\theta})$ (which may possibly not involve $\boldsymbol{\theta}$) is a prespecified $q$-dimensional vector with $q\geq p$ such that, under $H_0$, $\boldsymbol{\theta}^*$ is identified or over-identified. Examples of such vectors include $w(T;\boldsymbol{\theta})=(1,T,...,T^{q-1})^\top$ or $w(T;\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}g(T;\boldsymbol{\theta})$, where $``\nabla_{\boldsymbol{\theta}}"$ denotes the derivative with respect to $\boldsymbol{\theta}$. We then estimate $\boldsymbol{\theta }^*$ by
where $\|\cdot\|$ is the Euclidean norm, and
With the estimators $\{\widehat{\pi}_K(T_i,\boldsymbol{X}_i)\}_{i=1}^N$ of $\{\pi_0(T_i,\boldsymbol{X}_i)\}_{i=1}^N$ and $\widehat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$, we estimate $U_i$ by $ \widehat{U}_i=\widehat{\pi}_{K}(T_i,\boldsymbol{X}_i)m\{Y_i;g(T_i;\widehat{\boldsymbol{\theta }})\}$, for $i=1,\ldots,N$. Replacing the $U_i$'s in (ref) by the $\widehat{U}_i$'s, we have a feasible test statistic for $H_0$ based on
the corresponding estimators of the Cramér–von Mises (CM)-type statistic in (ref) and the Kolmogorov-Smirnov (KS)-type statistic in (ref) are, respectively,
where the supremum is calculated as the maximum value over a discretization of $\mathcal{T}$ in practice.
Remark 2. An alternative estimator of $\boldsymbol{\theta}^{\ast}$ can be constructed under $H_{0}$. Suppose that, under $H_{0}$, $\boldsymbol{\theta}^{\ast}$ is identified by the unique solution to the following optimization problem: \[ \boldsymbol{\theta}^{\ast}=\arg\min_{\boldsymbol{\theta}\in\Theta }CM(\boldsymbol{\theta}):=N\times\int_{\mathcal{T}}\left\{ \mathbb{E}\left[ U_{i}(\boldsymbol{\theta})\mathcal{H}(T_{i},t)\right] \right\} ^{2} f_{T}(t)dt, \] where $U_{i}(\boldsymbol{\theta}):=\pi_{0}(T_{i},\boldsymbol{X}_{i} )m\{Y_{i};g(T_{i};\boldsymbol{\theta})\}$. Let $\widehat{U}_{i} (\boldsymbol{\theta}):=\widehat{\pi}_{K}(T_{i},\boldsymbol{X}_{i} )m\{Y_{i};g(T_{i};\boldsymbol{\theta})\}$ and $\widehat{J}_{N} (t;\boldsymbol{\theta}):=N^{-1/2}\sum_{i=1}^{N}\widehat{U}_{i} (\boldsymbol{\theta})\mathscr{H}(T_{i},t)$. Under $H_{0}$, the estimator of $\boldsymbol{\theta}^{\ast}$ can be defined by
Therefore, the alternative test statistic is $\widehat{CM}_{N} (\widehat{\boldsymbol{\theta}}_{opt})$. However, seeking the global minimizer of $\widehat{CM}_{N}(\boldsymbol{\theta})$ is difficult as $\widehat{CM} _{N}(\boldsymbol{\theta})$ may not be differentiable, convex, and even continuous. For example, taking $m\{Y_{i};g(T_{i};\boldsymbol{\theta} )\}=\tau-\mathds{1}\{Y_{i}\leq g(T_{i};\boldsymbol{\theta})\}$ for QDRF, a unique solution to the problem does not exist. Under a stronger condition that $m(y;g)$ is differentiable in $g$, we establish the asymptotic results for both $\widehat{J}_{N}(t;\widehat{\boldsymbol{\theta}}_{opt})$ and $\widehat{CM}_{N}(\widehat{\boldsymbol{\theta}}_{opt})$ in section E in the supplementary file.\\ \noindentRemark 3. In order to estimate $\pi_0(T,\boldsymbol{X})$, Fong_Hazlett_Imai_2018 noted the moment conditions
which are special cases of our moment condition (ref). They then proposed estimating $\pi_0(T,\boldsymbol{X})$ by maximizing the empirical likelihood of $T$ and $\boldsymbol{X}$ under the constraints of the sample analogue of (ref) and estimating $\mathbb{E}\{Y^\ast(t)\}$ by a simple linear model. This can be considered as fixing $u_{K_1}(T)=(1,T)^\top$ and $v_{K_2}(\boldsymbol{X})=(1,\boldsymbol{X}^\top)^\top$, taking $m\{Y^\ast(t),g(t,\boldsymbol{\theta}^\ast)\}=Y^\ast(t) - g(t,\boldsymbol{\theta}^\ast)$ and $g$ as a simple linear model in the estimation method of Ai_Linton_Motegi_Zhang_cts_treat. However, the equation (ref) is of finite dimension and cannot nonparametrically identify $\pi_0(T,\boldsymbol{X})$. Hence, Fong_Hazlett_Imai_2018 imposed a parametric model for the stabilized weights to achieve consistent estimation. We adopt the estimator proposed by Ai_Linton_Motegi_Zhang_cts_treat that does not impose any parametric structure on the stabilized weights.\\ \noindentRemark 4. Once our specification test rejects the null model and no better parametric model can be proposed, several solutions are available. For example, one may consider the double robustness estimator. colangelo2020double estimate the average dose-response function $\mathbb{E}[Y^*(t)]$ based on the following double robustness representation:
where $K_h\left(T-t\right)$ is a kernel weighting observation $T$ with treatment value of approximately $t$ in a distance of $h$. They estimate both the general propensity score $f_{T|X}(t|\boldsymbol{X})$ and the outcome regression function $\mathbb{E}\left[Y|T=t,\boldsymbol{X}\right]$ using nonparametric techniques with cross-fitting. \\ An alternative solution to the misspecification of the dose-response function $g(t;\boldsymbol{\theta})$ is to consider a fully nonparametric specification $g(t)$, that is, to estimate $g(t)$ from the moment $\mathbb{E}[m(Y^*(t);g(t))]=0$ for all $t\in\mathcal{T}$. Under Assumption (ref), $g(t)$ can be identified through the conditional moment $\mathbb{E}[\pi_0(T,\boldsymbol{X})m(Y;g(T))|T]=0$. We can define the sieve minimum distance (SMD) estimator ai2003efficient of $g(T)$ by
where $\Sigma(T_i)$ is a user-specified weighting function, and
and $\mathcal{H}_{K_3}:=\{\lambda^{\top} u_{K_3}(T):\lambda\in\mathbb{R}^{K_3}\}$ is a linear sieve space. Ai_Linton_Motegi_Zhang_cts_treat established the large sample property for the nonparametric estimator of the average dose-response function $\widehat{E}[\widehat{\pi}_K(T,\boldsymbol{X})Y|T=t]$. \\ The extension of these methods to the general dose-response function including the quantile dose-response and the development of the corresponding large sample property is beyond the scope of this paper. \\ \noindentRemark 5. If the residual function $m(y;g(t;\boldsymbol{\theta}))$ is smooth in $(t,y)$, the sieve minimum distance (SMD) estimator of $\boldsymbol{\theta}^*$ developed by ai2003efficient is semiparametrically efficient with respect to the conditional model $\mathbb{E}[\pi_0(T,\boldsymbol{X})m\{T;g(T;\boldsymbol{\theta}^*)\}|T]=0$. This efficient estimation result can also be achieved based on our unconditional moment (ref); indeed, by replacing $\pi_0(T,\boldsymbol{X})$ with its estimate $\widehat{\pi}_K(T,\boldsymbol{X})$ and setting $w(T;\boldsymbol{\theta})$ to be a $q$-dimensional sieve basis, for example, $w(T;\boldsymbol{\theta})=(1,T,...,T^{q-1})^\top$, with $q\to \infty$ at an appropriate rate, it can be shown that the generalized method of moments (GMM) hansen1982large estimator of $\boldsymbol{\theta}^*$ constructed from (ref) is asymptotically equivalent to the SMD estimator (see ai2020simple for an analogous finding). However, if the residual function $m(y;g(t;\boldsymbol{\theta}))$ is non-smooth, it remains an open problem regarding whether the efficient estimation of $\boldsymbol{\theta}^*$ from the conditional model $\mathbb{E}[\pi_0(T,\boldsymbol{X})m\{T;g(T;\boldsymbol{\theta}^*)\}|T]=0$ can be established. For this reason, and to avoid introducing an extra tuning parameter $q$, we estimate $\boldsymbol{\theta}^*$ through (ref) with $w(T;\boldsymbol{\theta})$ as a fixed vector. }
This section studies the asymptotic properties of $\widehat{J}_N(\cdot)$, the test statistics $\widehat{CM}_N$ and $\widehat{KS}_N$.
To establish the asymptotic properties of $\widehat{J}_N(\cdot)$, $\widehat{CM}_N$ and $\widehat{KS}_N$, the following additional assumptions are imposed.
Assumption (ref) is essentially stating that the estimating equation is a.s. approximately satisfied; see pakes1989simulation and chen2003estimation. Assumption (ref) is needed to bound the asymptotic variance of the test statistic. Assumption (ref) (i) and (ii) impose sufficient regularity conditions on both the link function $g$ and residual function $m$. Assumption (ref) (iii) ensures that the variance of the test statistic is finite. Assumption (ref) is a stochastic equicontinuity condition, which is needed to establish the weak convergence of our test statistic; see andrews1994empirical. Again, this is satisfied by widely used residual functions such as $m\{y,g(t;\boldsymbol{\theta})\}=y-g(t;\boldsymbol{\theta})$ and $m\{y,g(t;\boldsymbol{\theta})\}=\tau-\mathds{1}\{y<g(t;\boldsymbol{\theta})\}$ discussed in Section (ref).
To aid presentation of the asymptotic properties of the test statistic, define the following quantities:
and
and
The next theorem establishes the weak convergence of $\widehat{J}_N(\cdot)$ and $\widehat{CM}_N$ under $H_0$.
The proof of Theorem (ref) is relegated to section B in the supplementary file. Similar to bierens1997asymptotic, chen1999consistent, it can be shown that $\int \{J_{\infty}(t)\}^2dF_T(t)$ can be written as an infinite sum of weighted (independent) $\chi_1^2$ random variables with weights depending on the unknown distribution of $(T_i,\boldsymbol{X}_i,Y_i)$. Hence, it is difficult to obtain the exact critical values. We suggest a simulation method to approximate the critical values for the null limiting distribution of $\widehat{CM}_N$; see Section (ref).
The effect of the vector $w(T;\boldsymbol{\theta})$ on the asymptotic property of our test statistic is reflected in the term $\psi(T_i,\boldsymbol{X}_i,Y_i,t)$. It is unclear which choice of $w(T;\boldsymbol{\theta})$ would minimize the variance of $\Sigma(t,t)$. A common choice is $w(T;\boldsymbol{\theta}) = \nabla_{\boldsymbol{\theta}}g(T;\boldsymbol{\theta})$. Then, the second and fourth terms of $\psi(T_i,\boldsymbol{X}_i,Y_i,t)$ in (ref) are canceled out, which also simplifies the calculation approximating the null limiting distribution in practice.
The next theorem shows that the proposed test statistic is more efficient than the infeasible test statistic constructed by using the true $\pi_0(T,\boldsymbol{X})$. Suppose that $\pi_0(T,\boldsymbol{X})$ was known, let $\widehat{\boldsymbol{\theta}}_0$ be the estimator of $\boldsymbol{\theta}^*$ constructed by using the true ratio function $\pi_0(T,\boldsymbol{X})$, which is defined by minimizing the following criterion function:
The infeasible test statistic for $H_0$ is then based on
Let
and
The following theorem establishes the weak convergence of $\widehat{J}_0(\cdot)$ under $H_0$ and shows that the asymptotic variance of the proposed test statistic $\widehat{J}_N(t)$ is smaller than that of $\widehat{J}_0(t)$ for any $t\in \mathcal{T}$.
The proof of Theorem (ref) is presented in section C in the supplementary file. In the estimation of the average treatment effects with binary and multiple treatments, it is a well-known paradox that using a nonparametric estimated propensity score is more efficient than using the true one; see Hirano03, chan2016globally, and lee2018efficient among others. Theorem (ref) shows that this is also the case for continuous treatments.
This section discusses two important special continuous treatment effect models, the average and quantile continuous treatment models. In the case of testing for the average dose-response model, that is,
against the alternative hypothesis
$m\left\{Y^*(t);g(t;\boldsymbol{\theta}^*)\right\} = Y^*(t)- g(t;\boldsymbol{\theta}^*)$, $U_i^{ADRF}=\pi_0(T_i,\boldsymbol{X}_i)\{Y_i-g(T_i;\boldsymbol{\theta}^*)\}$ and the test statistics for $H_0$ are
where
In this special case, the notations $\phi(T_i,\boldsymbol{X}_i;t)$, $\psi(T_i,\boldsymbol{X}_i,Y_i;t)$, and $\eta(T_i,\boldsymbol{X}_i,Y_i;t)$ in Theorem (ref) become
and
and
Then Theorem (ref) implies the following result.
In the case of testing for the quantile dose-response model, that is,
against the alternative hypothesis
$m\left\{Y^*(t);g(t;\boldsymbol{\theta}^*)\right\} = \tau-\mathbbm{1}\{Y^*(t)<g(t;\boldsymbol{\theta}^*)\}$, $U_i^{QDRF}=\pi_0(T_i,\boldsymbol{X}_i)\big[\tau-\mathbbm{1}\{Y_i<g(T_i;\boldsymbol{\theta}^*)\}\big]$, and the test statistics for $H_0$ are
where
Again, in this special case, the notations $\phi(T_i,\boldsymbol{X}_i;t)$, $\psi(T_i,\boldsymbol{X}_i,Y_i;t)$, and $\eta(T_i,\boldsymbol{X}_i,Y_i;t)$ in Theorem (ref) become
and
and
Then, Theorem (ref) implies the following result.
This section studies the asymptotic distribution of $\widehat{J}_N(\cdot)$ under the fixed and Pitman local alternatives. The Pitman local alternative is given by
where $\int \{\delta (t)\}^2dF_T(t)<\infty$. With Assumption (ref), $H_L$ can be represented by
which deviates from the null model at the rate of $O(N^{-1/2})$. Let $\boldsymbol{\theta}^*$ be the limit of $\boldsymbol{\theta}^*_N$ as $N\rightarrow\infty$, hence it solves the following equation:
Define
The following theorem gives the asymptotic distribution of $\widehat{J}_N(\cdot)$ under the local alternative $H_L$ and the fixed alternative $H_1$.
Comparing Theorem (ref) (ii) to Theorem (ref) (ii), we see that our test statistic is able to detect the local alternatives deviated from the null model at the rate of $O(N^{-1/2})$.
We know from Theorem (ref) that $\widehat{CM}_{N}$ converges in distribution to $\int\{J_{\infty} (t)\}^{2}\,dF_{T}(t)$. Using techniques similar to those in bierens1997asymptotic and chen1999consistent, one can show that $\int\{J_{\infty}(t)\}^{2}\,dF_{T}(t)$ is an infinite sum of weighted (independent) $\chi_{1}^{2}$ random variables, where the weights depend on the unknown distribution of the $(\boldsymbol{X}_{i},T_{i},Y_{i})$'s (see also li2003consistent). Obtaining the exact critical values is difficult and we here propose a simulation method to approximate the null limiting distribution. The method is a special case of the exchangeable bootstrap praestgaard1993exchangeably,van1996weak,chernozhukov2013inference,Donald2014Estimation. Specifically, we first generate $B$ sets of $N$ independent standard normal random variables $w_{1,b},\ldots,w_{N,b}$, for $b=1,\ldots,B$ and $B$ a large enough integer. Then we define
where $\widehat{\eta}(T_{i},\boldsymbol{X}_{i},Y_{i};t) = \widehat{U} _{i}\mathscr{H}(T_{i},t) - \widehat{\phi}(T_{i},\boldsymbol{X}_{i};t) - \widehat{\psi}(T_{i},\boldsymbol{X}_{i},Y_{i};t)$, with $\widehat{\phi} (T_{i},\boldsymbol{X}_{i};t)$ and $\widehat{\psi}(T_{i},\boldsymbol{X} _{i},Y_{i};t)$ respectively some consistent nonparametric plug-in estimators of $\phi(T_{i},\boldsymbol{X}_{i};t)$ and $\psi(T_{i},\boldsymbol{X}_{i} ,Y_{i};t)$ defined above in Theorem (ref), for example the additive penalized spline estimator(see ruppert2003 for example) or the series estimator used in Donald2014Estimation.
It is easy to see that $\mathbb{E}^{*}\{w_{i,b}\widehat{\eta}(T_{i} ,\boldsymbol{X}_{i},Y_{i};t)\} = 0$ and $\mathbb{E}^{*}\{w^{2}_{i,b} \widehat{\eta}(T_{i},\boldsymbol{X}_{i},Y_{i};t)\widehat{\eta}(T_{i} ,\boldsymbol{X}_{i},Y_{i};$ $t^{\prime})\}= \widehat{\eta}(T_{i} ,\boldsymbol{X}_{i},Y_{i};t)\widehat{\eta}(T_{i},\boldsymbol{X}_{i} ,Y_{i};t^{\prime})$, for $i=1,\ldots,N$, $b=1,\ldots,B$ and all $t,t^{\prime }\in\mathcal{T}$, where $\mathbb{E}^{*}\{\cdot\}$ is the conditional expectation given the data $(T_{i},\boldsymbol{X}_{i},Y_{i})_{i=1}^{N}$. Because $\widehat{\eta}$ is a consistent estimator of $\eta$, then $\widehat{J}^{*}_{N,b}(\cdot)$ has the same asymptotic behavior as $\widehat {J}_{N}(\cdot)$ for $b=1,\ldots,B$. Then, we can approximate the limiting distributions of $\widehat{CM}_{N}$ and $\widehat{KS}_{N}$ under $H_{0}$, respectively, by \[ \widehat{CM}_{N,b}^{*} = \frac{1}{N}\sum^{N}_{i=1}\big\{\widehat{J}_{N,b} ^{*}(T_{i})\big\}^{2} \quad \text{and} \quad \widehat{KS}_{N,b}^{*} = \sup_{t\in\mathcal{T}}\big|\widehat{J}_{N,b} ^{*}(t)\big|, \] for $b=1,\ldots,B$. That is, we can approximate the $p$-value for the CM-type statistic by $B^{-1} \sum^{B}_{b=1}\mathbbm{1}(\widehat{CM}_{N,b}^{*}\geq\widehat{CM}_{N})$ and that for the KS-type statistic by $B^{-1} \sum^{B}_{b=1}\mathbbm{1}(\widehat{KS}_{N,b}^{*}\geq\widehat{KS}_{N})$.
The large-sample properties of the proposed estimator hold for a range of values of $K_1$ and $K_2$. This presents a dilemma for applied researchers, who have only one finite sample. Too little smoothing yields a large variance and too much smoothing yields a large bias. Therefore, applied researchers would benefit from guidance on the choice of $K_1$ and $K_2$. In this section, we propose a cross-validation method for choosing the smoothing parameters $K_1$ and $K_2$. Specifically, we split the data set into $F$ sets (say $F=5$ or 10), and select $K_1$ and $K_2$ that minimize the following quantity
where $S_j$ denotes the $j$th set of data of $T,\boldsymbol{X}$ and $Y$, $|S_j|$ denotes the number of individuals in the set $S_j$, and for $j=1,\ldots,F$, $$ \widehat{\boldsymbol{\theta }}^{(-j)} = \arg \min_{\boldsymbol{\theta}\in\Theta}\left\|M_N^{(-j)}\left(\boldsymbol{\theta},\widehat{\pi}_{K}^{(-j)}\right)\right\|\,, $$ where $$ M_N^{(-j)}\left(\boldsymbol{\theta},\widehat{\pi}_{K}^{(-j)}\right):= \frac{1}{N}\sum_{i\notin S_j}\widehat{\pi}_{K}^{(-j)}(T_{i},\boldsymbol{X}_{i})m\{Y_i;g(T_i,\boldsymbol{\theta})\}w(T_i;\boldsymbol{\theta})\,, $$ with $\widehat{\pi}_{K}^{(-j)}(T_{i},\boldsymbol{X}_{i})$ obtained in a method identical to that introduced in Section (ref) via (ref) and (ref), but excluding samples in $S_j$.
To assess the performance of our specification test method, we conducted Monte Carlo simulation studies on the following four data generating processes (DGPs):
where $\xi$ and $\epsilon$ are independent standard normal random variables, and $X$ is a uniform random variable supported on $[0,1]$. For all the four scenarios, we considered the two-sided hypothesis testing in (ref), where $m\{Y^*(t);g(t;\boldsymbol{\theta}^*)\}=Y^*(t)-g(t;\boldsymbol{\theta}^*)$ (average) and $m\{Y^*(t);g(t;\boldsymbol{\theta}^*)\}=0.5 - \mathbbm{1}\{Y^*(t)<g(t;\boldsymbol{\theta}^*)\}$ (median), and $$ g\{t;(\theta^*_0,\theta^*_1)\} = \theta^*_0 + \theta^*_1 t\,. $$ We take the vector $w(T;\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}g(T;\boldsymbol{\theta})$ and adopt the algorithm of de2019smoothed for estimating the quantile dose-response function to overcome the computational difficulties with the discontinuity of the indicator function.
Clearly, $H_0$ is true for {\bf DGP0-L} and {\bf DGP0-NL}, but fails for {\bf DGP1-L} and {\bf DGP1-NL}. For each case, we generated 1000 samples of size 100, 200, and 500. The number of samples for the simulation-based approximation of the limiting process is $B=500$ and the number of folds in the cross-validation (ref) was taken to be $F=10$. We compared the three commonly used weight functions $\mathscr{H}$ that are mentioned in Section (ref), namely logistic, cosine-sine, and indicator functions. Specifically, for the logistic weight function, we took the constant $c=5$. We tested all models using both CM-type and KS-type statistics. The results of the two methods are similar; here, we present those of the CM-type statistic. Results of the KS-type one can be found in the supplementary materials of this paper.
Tables (ref) and (ref) summarize the empirical rejection probabilities computed at significance levels $1\%, 5\%$, and $10\%$ for each case, which respectively show the estimated sizes (DGP0-L and DGP0NL) and the estimated powers (DGP1-L and DGP1-NL) of our test method.
We can see from Table (ref) that the estimated sizes of our method with cosine-sine and indicator weight functions are quite close to the nominal sizes from $N=100$ to $500$ for all cases. The estimated sizes when using the logistic weight function are obviously over-sized when the sample size is small, especially for nonlinear $\boldsymbol{X}$ cases, but they also improve as the sample size increases and are close to the nominal sizes when $N=500$.
From Table (ref), we observe that all tests become more and more powerful as $N$ or significance level increases and reach a considerably high power level when $N=200$.
Overall, the simulation studies confirmed our asymptotic theorems and showed that, in practice, the cosine-sine and indicator weight functions might perform better than the logistic one for nonlinear $\boldsymbol{X}$ cases.
In this section, we applied our method to examine the model assumption made on the U.S. presidential campaign data in Ai_Linton_Motegi_Zhang_cts_treat. The data have been analyzed several times in the treatment effect literature Urban_Niebler_2014,Fong_Hazlett_Imai_2018, where the interest was to explore the casual relationship between advertising and campaign contributions. The treatment of interest is the number of political advertisements aired in each zip code from non-competitive states, which ranges from 0 to 22379 across $N=16265$ zip codes.
The data were first analyzed by Urban_Niebler_2014, who used a binary model to compare the campaign contributions of the 5230 zip codes that received more than 1000 advertisements with those of the other 11035 zip codes that received less than 1000 advertisements. Their research suggested that advertising in non-competitive states had a significant casual effect on the level of campaign contributions.
By contrast, Ai_Linton_Motegi_Zhang_cts_treat considered the treatment variable (number of political advertisements) as continuous and assumed that
where the observed outcome $Y^*(T) = \log(\text{Contribution}+1)$ and $T = \log(\#\text{ads}+1)$, where \#ads denotes the number of advertisements. The covariates $\boldsymbol{X}$ considered were $$ \boldsymbol{X} =
\,. $$ The definition of each covariate is almost self-explanatory, and one can refer to \cite{Fong_Hazlett_Imai_2018} for more details. \cite{Ai_Linton_Motegi_Zhang_cts_treat} found that the 95\% confidence intervals for $\beta_2$ and $\beta_3$ were respectively $[-0.025, 0.232]$ and $[-0.025, 0.001]$, indicating that no significant causal link between advertising and campaign contributions was found from the linear model. Similar results were also reported by Fong_Hazlett_Imai_2018. The authors then concluded that such opposing results from binary models and continuous linear models suggested a rather complex relationship between advertising and campaign contributions.
We reached the same conclusion in our data analysis. Indeed, when we applied our method with logistic, cosine-sine, and indicator weight functions with a $B=500$ simulation-based approximation to test the model in (ref), all the methods rejected the model with the $p$-values equal to 0.
We examined the histogram of the original campaign contribution data and the number of advertisements $T$. From the first row of Figure (ref), we can see that both histograms are highly right-skewed. That is, they are not likely to fit any linear models. We then conducted a log-transformation, as in Ai_Linton_Motegi_Zhang_cts_treat. However, the results were similar. To make the data more likely to fit a linear model, we searched across Box-Cox transformations of the response data of the form $\text{BoxCox}(\text{Contribution},\lambda_1,\lambda_2):=\{(\text{Contribution}+\lambda_2)^{\lambda_1}-1\}/\lambda_1$ w.r.t. $\lambda_1, \lambda_2$ to find a transformation of the contribution whose sample quantiles have the largest correlation with those of a standard normal distribution. This yielded $(\tilde{\lambda}_1,\tilde{\lambda}_2) = (0.1397,0.0176)$. We then take
so that the minimum response data is 0. The histogram of $Y$ is shown in the bottom left of Figure (ref). We can see that the transformed data remain highly right-skewed. However, now, it appears to be a truncated normal distribution.
Now, it seems more reasonable to assume a Tobit model for the data. Specifically, let $$ Y(t) = \boldsymbol{\beta}^\top \boldsymbol{t} + \epsilon\,, $$ for some unknown parameter $\boldsymbol{\beta}$ in a compact set in $\mathbb{R}^p$, where $\boldsymbol{t} = (1,t,t^2,\ldots,t^{p-1})^\top$ and $\epsilon$ is a normal random variable with mean 0 and unknown variance $\sigma^2$. We test a Tobit linear model on the potential outcome: $$ Y^\ast(t) =
$$ The details of the estimation of this model and the test statistics can be found in the supplementary material of this paper. We then tested the Tobit model with several different transformations of the treatment data \#ads and the polynomial order $p$. We found the Tobit model with $T = \log(\log(\log(\#ads+1)+1)+2)$ and $p=5$ gives the most reasonable results. The corresponding $p$-values are shown in Table (ref), and the estimated model is depicted in Figure (ref). It indicates that the campaign contribution increases rapidly with a relatively small increase in the number of advertisements from 0, and then the improvement gradually becomes marginal.
Supplementary materials are only for online publication. The supplementary file contains the simulation results of the KS-type statistic, the details of the estimation and the test statistics for the Tobit model used in section (ref), the proofs of Theorems (ref), (ref), (ref), and the asymptotic properties of $\widehat{J}_N(t;\widehat{\boldsymbol{\theta}}_{opt})$ and $\widehat{CM}_N(\widehat{\boldsymbol{\theta}}_{opt})$. \if00 {
We thank two anonymous referees, an associate editor, and the editor Christian Hansen for constructive comments which lead to substantial improvement of the paper. The first author, Wei Huang's research was supported by the Professor Maurice H. Belz Fund of the University of Melbourne. The second author, Oliver Linton, acknowledges Cambridge INET for financial support. The final author, Zheng Zhang, acknowledges financial support from the National Natural Science Foundation of China through project 12001535, and the fund for building world-class universities (disciplines) of the Renmin University of China. } \fi
\if10 \fi