Nonparametric Estimation of the Continuous Treatment Effect with Measurement Error

Identifying and estimating the causal effect of a treatment or policy from observational studies is of great interest to economics, social science, and public health researchers. There, confounding issues usually exist (i.e. individual characteristics are related to both the treatment selection and the potential outcome), making the causal effect not directly identifiable from the data. Early studies focused on whether an individual receives the treatment or not (e.g. Rosenbaum and Rubin, 1983, 1984; Hahn, 1998; Hirano et al., 2003). More recently, as opposed to such binary treatments, researchers have been investigating the causal effect of a continuously valued treatment, where the effect depends not only on the introduction of the treatment but also on the intensity; see Hirano and Imbens (2004); Galvao and Wang (2015); Kennedy et al. (2017); Fong et al. (2018); Huber et al. (2020); Dong et al. (2021); Ai et al. (2021, 2022), and Huang et al. (2021), among others. However, all these methods require the treatment data to be measured without errors.

In many empirical applications, the treatment and confounding variables may be inaccurately observed. For example, Mahajan (2006) studied the binary treatment effect when the observed treatment indicators are subject to misclassification; see also Lewbel (2007) and Molinari (2008), among others. Battistin and Chesher (2014) investigated the extent to which confounder measurement error affects the analysis of the treatment effect. Continuous treatment variables are also likely to be measured with error in practice. For example, in the Epidemiologic Study Cohort data from the first National Health and Nutrition Examination (NHANES-I) (see Carroll et al., 2006), over 75% of the variance in the fat intake data is made up of measurement error. However, to the best of our knowledge, no existing work has considered the identification and estimation of the causal effect from such error-contaminated continuous treatment data. To bridge this gap in the body of knowledge, we focus on continuous treatment data measured with classical error; that is, instead of observing the received treatment, researchers only observe the sum of the treatment and a random error.

In particular, we study a causal parameter of primary interest to scholars in the literature, the average dose-response function (ADRF). It is defined as the population mean of the individual potential outcome corresponding to certain levels of the treatment. To resolve the confounding problem, we adopt the most widely imposed assumption in the literature on treatment effects, the unconfoundedness condition (e.g. Rosenbaum and Rubin, 1983, 1984; Hirano and Imbens, 2004), which assumes that the assignment of the treatment is independent of the potential outcome of interest given a set of observable covariates.

Perhaps, the most straightforward approach to this problem is based on some parametric model specifying how the outcome relates to confounders and the treatment. Then one may apply the parametric methods to measurement error data in the literature (see e.g. Carroll et al., 2006). However, the parametric approach suffers from the model misspecification problem and does not incorporate available information on the treatment mechanism. Thus, this paper focuses on robust nonparametric estimation of ADRF.

In the literature on error-free continuous treatment data, nonparametric estimation of ADRF under unconfoundedness assumption has been extensively studied. For example, Galvao and Wang (2015) identified the ADRF via an unconditional weighted expectation, with the weighting function being the ratio of two conditional densities of the treatment variable. They estimated these two conditional densities separately and then constructed the estimator for the ADRF. However, it is well known that such a ratio estimator can be unstable owing to its sensitivity to the denominator estimation (Kang and Schafer, 2007). To improve the robustness of estimation, Kennedy et al. (2017) developed a doubly robust estimator for ADRF by regressing a doubly robust mapping on the treatment; see more detailed discussion in Remark 2 of Section 3.1. Ai et al. (2021) identified ADRF via a weighted conditional expectation on the treatment variable, where the weighting function is the ratio of the treatment variable’s marginal density and its conditional density given the confounders. They estimated the weighting function (but not the two densities in the ratio separately) by maximising entropy subject to some moment restrictions that identify the weighting function, then obtained the nonparametric weighted estimator of the ADRF. The idea of estimating the density ratio directly has also been exploited in the literature in (bio)statistics based on parametric modeling (Qin, 1998), semiparametric modeling (Cheng and Chu, 2004), and machine learning with an augmented dataset (Díaz et al., 2021).

However, no existing methods for the error-free treatment data can be easily extended to the error-contaminated treatment data. Indeed, under the unconfoundedness assumption, all those methods depend on estimating a weighting function with the treatment’s conditional density given the confounders in the denominator. Its estimation is a critical step and full of challenges in the presence of measurement error. Although the nonparametric method of estimating the (conditional) density of a variable measured with the classical error has been developed (see e.g. Fan, 1991a, b; Meister, 2006), estimating the density and plugging it into the denominator suffers from the instability issue mentioned above. Such an issue can be even worse when the treatment is measured with errors. Unfortunately, existing techniques of improving the robustness cannot apply to the measurement error data; for example, Ai et al.’s 2021 approach heavily relies on the moment information of a general function of the treatment variable, which is hard to obtain in the presence of measurement error. To our best knowledge, estimating the moment of a general function of a variable measured with errors and its associated theoretical behaviour are challenging problems and remain open in the literature (see Hall and Lahiri, 2008 and our discussion in section 3.1 for more details).

We propose a broad class of novel and robust nonparametric estimators for the ADRF when the treatment data are measured with error. We first represent the ADRF as in Ai et al. (2021), i.e. a weighted expectation of the observed outcome on the treatment variable, where the weighting function is a ratio of the treatment variable’s marginal density and its conditional density given the confounders. Then we propose a novel approach to identify the weighting function with no need to estimate the moment of a general function of the treatment. Specifically, we propose a local empirical likelihood technique that only requires estimating the conditional moment of a general function of the error-free confounders given each fixed value of the treatment. A consistent nonparametric estimator of such a conditional moment can be obtained by using a local constant deconvolution kernel estimator in the literature on errors-in-variables regression; see Fan and Truong (1993); Carroll and Hall (2004); Hall et al. (2007); Delaigle et al. (2009), and Meister (2009), among others. Moreover, by choosing different criterion functions for the local generalised empirical likelihood, the proposed method produces a wide class of weights, including exponential tilting, empirical likelihood and generalised regression as important special cases. Once the weights are consistently estimated, we construct our estimator of the ADRF by estimating the weighted conditional expectation of the observed outcome using another local constant deconvolution kernel estimator.

Based on our identification theorem and the literature on the local constant deconvolution kernel estimator, it is not hard to see the consistency of our estimator. However, the asymptotic behaviour of a local constant deconvolution kernel estimator incorporated with such local generalised empirical likelihood weight has never been investigated. We show the asymptotic bias of our estimator consists of two parts: the one from estimating the weights and the one from estimating the weighted conditional expectation, but it can achieve the same convergence rate as that in the error-free case. The asymptotic variance depends on the type of measurement error. We study the most commonly used two types of measurement errors: the ordinary smooth and the supersmooth ones. The convergence rate of the variance is slower than that in the error-free case but can achieve the optimal rate for the nonparametric estimation in measurement error literature.

Moreover, without imposing additional conditions, it is well known that the explicit expression of the asymptotic variance of the deconvolution kernel for supersmooth error (e.g. Gaussian error) is not derivable; even the exact convergence rate remains unknown. Therefore, it is common to use residual bootstrap methods to make the inference (see e.g. Fan, 1991a; Delaigle et al., 2009, 2015). However, residual bootstrap methods do not work well in our treatment effect framework. We thus provide the asymptotic linear expansion of our ADRF estimator. Based on this, we propose an undersmoothing (pointwise) confidence interval for the ADRF.

A data-driven method to select the bandwidth for a local polynomial deconvolution kernel estimator is challenging. The only method in the literature to our best knowledge is the simulation-extrapolation (SIMEX) method proposed by Delaigle and Hall (2008). However, the linear back-extrapolation method suggested by the authors performed unstably in our numerical studies. We, thus, propose a new local constant extrapolant function. Monte Carlo simulations show that our estimator is stable and outperforms the naive estimator that does not adjust to the measurement error. We also demonstrate the practical value of our estimator by studying the causal effect of fat intake on breast cancer using the Epidemiologic Study Cohort data from NHANES-I (Carroll et al., 2006).

The remainder of the paper is organised as follows. We introduce the basic framework and notations in Section 2. Section 3 presents our estimation procedure, followed by the asymptotic studies of our estimator in Section 4. We discuss the method of selecting our smoothing parameters in Section 5. Finally, we present our simulation results and real data analysis in Section 6.

We consider a continuously valued treatment in which the observed treatment variable is denoted by $T$ with the probability density function $f_{T}(t)$ and support $\mathcal{T}\subset\mathbb{R}$. Let $Y^{\ast}(t)$ denote the potential outcome if one was treated at level $t$ for $t\in\mathcal{T}$. In practice, each individual can only receive one treatment level $T$ and we only observe the corresponding outcome $Y:=Y^{\ast}(T)$. We are also given a vector of covariates $\boldsymbol{X}\in\mathbb{R}^{r}$, with $r$ a positive integer, which is related to both $T$ and $Y^{\ast}(t)$ for $t\in\mathcal{T}$. Moreover, we consider the situation in which the treatment level is measured with classical error, that is, instead of observing $T$, we observe $S$ such that

where $U$ is the measurement error, independent of $T,\boldsymbol{X}$ and $\{Y^{\ast}(t)\}_{t\in\mathcal{T}}$, and its characteristic function $\phi_{U}$ is known; see Remark 1 for the case in which $\phi_{U}$ is unknown.

The goal of this paper is to nonparametrically estimate of the unconditional ADRF $\mu(t):=\mathbb{E}\{Y^{\ast}(t)\}$ from an independent and identically distributed (i.i.d.) sample $\{S_{i},\boldsymbol{X}_{i},Y_{i}\}_{i=1}^{N}$ drawn from the joint distribution of $(S,\boldsymbol{X},Y)$.

$Y^{\ast}(t)$ is never observed simultaneously for all $t\in\mathcal{T}$ or even on a dense subset of $\mathcal{T}$ for any individual, but only at a particular level of treatment, $T$. Thus, to identify $\mu(t)$ from the observed data, the following assumption is imposed in most of the treatment effect literature (e.g. rubin1990comment; Kennedy et al., 2017; Ai et al., 2021, 2022; D’Amour et al., 2021).

Under Assumption 1, for every fixed $t\in\mathcal{T}$, $\mu(t)$ can be identified as follows:

where

with $f_{T}$ and $f_{T|X}$ being the density function of $T$ and the conditional density of $T$ given $\boldsymbol{X}$, respectively. The function $\pi_{0}(t,\boldsymbol{x})$ is called the stabilised weights in Ai et al. (2021).

If $T$ is fully observable and $\pi_{0}(t,\boldsymbol{x})$ is known, estimating $\mu(t)$ in (2) is reduced to a standard regression problem. For example, $\mu(t)$ can be consistently estimated using the Nadaraya–Watson estimator:

where $L(\cdot)$ is a prespecified univariate kernel function such that $\int_{-\infty}^{\infty}L(x)dx=1$, $h$ is the bandwidth. However, we do not observe $T$ but only $S$ in (1) and $\pi_{0}(t,\boldsymbol{x})$ is also unknown in practice. We address these issues in the next section.

To overcome the problem that the $L\{(t-T_{i})/h\}$’s are not empirically accessible, we apply the deconvolution kernel approach (e.g. Stefanski and Carroll, 1990; Fan and Truong, 1993). This method is often used in nonparametric regression in which the covariates are measured with classical error, as in (1), and the idea is introduced as follows. The density of $S$ is the convolution of the densities of $T$ and $U$, meaning that $\phi_{S}(w)=\phi_{T}(w)\phi_{U}(w)$, where $\phi_{S}$ and $\phi_{T}$ are the characteristic functions of $S$ and $T$, respectively. We consider $U$ with $\phi_{U}(w)\neq 0$ for all $w\in\mathbb{R}$. Using the Fourier inversion theorem, if $|\phi_{T}|$ is integrable, we have

This inspired Stefanski and Carroll (1990) to estimate $f_{T}$ by $\widehat{f}_{T,h}(t):=(Nh)^{-1}\sum^{N}_{i=1}L_{U}\{(t-S_{i})/h\}$, where

with $\phi_{L}$ the Fourier transform of the kernel $L$, which aims to prevent $\widehat{f}_{T,h}$ from becoming unreliably large in its tails.

Based on this idea, Fan and Truong (1993) proposed a consistent errors-in-variables regression estimator by replacing the $L\{(t-T_{i})/h\}$’s in (4) with the $L_{U}\{(t-S_{i})/h\}$’s. In our context, an errors-in-variables estimator of $\mu(t)$ is

Note that $U\perp(T,\boldsymbol{X},Y)$, we have

where the last equation comes from the Fourier inversion theorem. Using this property, $\widetilde{\mu}(t)$ has the same asymptotic bias as that of $\mu_{NW}(t)$, which shrinks to zero as $h\to 0$. Then, to verify its consistency to $\mu(t)$, it suffices to show that its asymptotic variance decays to zero as $N\rightarrow\infty$, using a straightforward extension of the proof in Fan and Truong (1993).

However, the challenge is that $\pi_{0}$ is unknown in practice. We next show how to estimate $\pi_{0}(t,\boldsymbol{X})$ from the error-contaminated data $(S_{i},\boldsymbol{X}_{i},Y_{i})$, $i=1,\ldots,N$.

In this section, we establish the $L_{\infty}$ and $L_{2}$ convergence rates of $\widehat{\pi}(t,\cdot)$ for every fixed $t\in\mathcal{T}$. We then investigate the asymptotic behaviour of the proposed ADRF estimator $\widehat{\mu}(t)$. Note that $\widehat{\mu}(t)$ (resp. $\widehat{\pi}(t,\cdot)$) is a nonparametric estimator and that its asymptotic behaviour is affected by the asymptotic bias and variance, which are respectively defined as the expectation and variance of the limiting distribution of $\widehat{\mu}(t)-\mu(t)$ (resp. $\widehat{\pi}(t,\cdot)-\pi_{0}(t,\cdot)$). Based on (8), we will show that the asymptotic biases of the two estimators are the same as their counterparts in the error-free case. That is, they depend on the smoothness of $\pi_{0}$, $\mu$, and the density of $T$, and the approximation error based on the sieve basis $u_{K}$. In particular, the following conditions are required:

Assumption 3 (i) restricts the covariates and the treatment to be bounded. This condition is commonly imposed in the nonparametric regression literature. Assumption 3 (i) can be relaxed if we restrict the tail distributions of $\boldsymbol{X}$ and $T$. For example, Chen et al. (2008, Assumption 3) allowed the support of $\boldsymbol{X}$ to be the entire Euclidean space but imposed $\int_{\mathbb{R}^{r}}(1+|\boldsymbol{x}|^{2})^{\omega}f_{X}(\boldsymbol{x})d\boldsymbol{x}<\infty$ for some $\omega>0$.

Assumption 3 (ii) is a strict positivity condition requires every subject having certain chance of receiving every treatment level regardless of covariates. This condition is also imposed in a large body of literature in the absence of measurement error (see e.g. Kennedy et al., 2017, Assumption 2 and D’Amour et al., 2021, Assumption 3), particularly when no restrictions are imposed on the potential outcome distribution. This condition can be relaxed if other smoothness conditions are imposed on the potential outcome distribution (Ma and Wang, 2020), or if different target parameters are considered; for example, Muñoz and Van Der Laan (2012) studied the estimation of a stochastic intervention causal parameter, defined by $\mathbb{E}[\mathbb{E}[Y|T+a(\boldsymbol{X}),\boldsymbol{X}]]$, based on a weaker positivity condition, i.e. $\sup_{t\in\mathcal{T}}\{f_{T|X}(t-a(\boldsymbol{X})|\boldsymbol{X})/f_{T|X}(t|\boldsymbol{X})\}<\infty$ a.e., where $a(\boldsymbol{X})$ is a user specified intervention function. Díaz and van der Laan (2013) studied the estimation of a conditional causal dose-response curve defined by $\mathbb{E}[\mathbb{E}[Y|T=t,\boldsymbol{X}]|\boldsymbol{Z}]$, where $\boldsymbol{Z}\subset\boldsymbol{X}$ is a subset of observed covariates, based on a weaker positivity condition, i.e. $\sup_{t\in\mathcal{T}}\{b(t,\boldsymbol{Z})/f_{T|X}(t|\boldsymbol{X})\}<\infty$ a.e., for a user specified weight function $b(t,\boldsymbol{Z})$. Although Assumption 3 (ii) is not the mildest condition in the literature, we maintain it throughout this paper owing to its technical benefits, especially in the presence of measurement error.

Assumption 4 includes smoothness conditions required for nonparametric estimation. Under Assumption 1, the parameter of interest $\mu(t)=\mathbb{E}[Y^{*}(t)]$ can be also written as $\mu(t)=\mathbb{E}\left[\pi_{0}(t,\boldsymbol{X})Y|T=t\right]=\mathbb{E}\left[f_{T|Y,X}(t|Y,\boldsymbol{X})Y/f_{T|X}(t|\boldsymbol{X})\right]$. Note that Assumption 4 (i) implies that $t\mapsto f_{T|Y,X}(t|Y,\boldsymbol{X})/f_{T|X}(t|\boldsymbol{X})$ is third-order continuously differentiable almost surely. Furthermore, using Leibniz integral rule and Assumption 4 (ii), we have that the target parameter $t\to\mu(t)$ is third-order continuously differentiable.

Assumption 5 (i) is used to control the complexity (measured by the uniform entropy integral) of the function class $\{\pi_{0}(t,\boldsymbol{x}),\boldsymbol{x}\in\mathcal{X}\}$ such that it forms a Donsker class and the empirical process theory can be applied (Van Der Vaart et al., 1996, Corollary 2.7.2). Despite of its stringency, the smoothness condition of this type is commonly adopted in the literature of nonparametric inference, see Chen et al. (2008, Assumption 4 (i)) and Fan et al. (2021, Condition E.1.7).

Assumption 5 (ii) requires the sieve approximation error of $\rho^{\prime-1}\left\{\pi_{0}(t,\boldsymbol{x})\right\}$ to shrink at a polynomial rate. This condition is satisfied for a variety of sieve basis functions. For example, it can be satisfied with $\alpha=+\infty$ if $\boldsymbol{X}$ is discrete, and with $\alpha=s/r$ if $\boldsymbol{X}$ is continuous and $u_{K}(\boldsymbol{x})$ is a power series or a B-spline, where $s$ is the smoothness of the approximand and $r$ is the dimension of $\boldsymbol{X}$. Assumption 7 imposes a similar sieve approximation error for $m(t,\boldsymbol{x})$.

Assumption 6 (i) rules out near multicollinearity in the approximating basis functions, which is common in the sieve regression literature. Assumption 6 (ii) is satisfied with $\zeta(K)=O(K)$ if $u_{K}(\boldsymbol{x})$ is a power series and with $\zeta(K)=O(\sqrt{K})$ if $u_{K}$ is a B-spline (Newey, 1997). Assumption 8 imposes the boundedness conditions on the moment of the response variable, which are also standard in the errors-in-variables problem (e.g. Fan and Truong, 1993; Delaigle et al., 2009). This condition is needed for deriving the asymptotic distribution of the proposed estimator by applying the Lyapunov central limit theorem.

Depending on the type of the distribution of $U$ and decaying rates of $h_{0}$ and $h$, the asymptotic variance of our estimator differs. This is different from the error-free case. We consider two types of $U$: the ordinary smooth case and supersmooth case, which are standard in the literature of errors-in-variables problem (see e.g. Fan and Truong, 1993, Delaigle et al., 2009, and Meister, 2009, among others, for more details).

An ordinary smooth error of order $\beta\geq 1$ satisfies