{'fraction_non_alphanumeric': 0.061906946422966416, 'fraction_numerical': 0.0238790834681411, 'mean_word_length': 3.9096728307254622, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonparametric covariate-adjusted regression

We consider nonparametric estimation of a regression curve when the data are observed with multiplicative distortion which depends on an observed confounding variable. We suggest several estimators, ranging from a relatively simple one that relies on restrictive assumptions usually made in the literature, to a sophisticated piecewise approach that involves reconstructing a smooth curve from an estimator of a constant multiple of its absolute value, and which can be applied in much more general scenarios. We show that, although our nonparametric estimators are constructed from predictors of the unobserved undistorted data, they have the same first order asymptotic properties as the standard estimators that could be computed if the undistorted data were available. We illustrate the good numerical performance of our methods on both simulated and real datasets.

1. Introduction. We consider nonparametric estimation of a regression curve m(x) = E(Y |X = x) when X and Y are observed with multiplicative distortion induced by an observed confounder U . Specifically, we observeX,Ỹ and U , whereỸ = ψ(U ) Y ,X = ϕ(U ) X, ψ and ϕ are unknown functions and U is independent of X and Y . This model is known as a covariate-adjusted regression model. It was introduced by Şentürk and Müller (2005a) to generalize an approach commonly employed in medical studies, where the effect of a confounder U , for example body mass index, is often removed by dividing by U . Motivated by the fibrinogen data on haemodialysis patients, whereỸ was fibrogen level,X was serum transferrin level, and U was body mass index, Şentürk and Müller (2005a) pointed that although it is often reasonable to assume that the effect of U is multiplicative, it does not need to be proportional to U , and a more flexible model is obtained by allowing for distortions represented by the functions ϕ and ψ. More generally, this model is useful to describe the relationship between variables that are influenced by a confounding variable, and see if this relationship still exists once the effect of the confounder has been removed.
A number of authors have suggested estimators of the curve m in various parametric settings. Linear regression models were considered by Şentürk andMüller (2005a, 2006) and Şentürk and Nguyen (2006), who generalized them to varying coefficient models (Şentürk, 2006) and generalized linear models (Şentürk and Müller, 2009). A more general nonlinear regression model was suggested by Cui et al. (2009) and Zhang et al. (2012), and in Zhang et al. (2013), the authors considered a partially linear model, where the linear part is observed with multiplicative distortions.
In this work, we propose more flexible nonparametric estimators of the regression function m, which not only relax the parametric assumptions imposed in the existing literature, but also significantly weaken some of the strong assumptions on the curves ϕ and ψ and on the distribution of the data made by previous authors. In particular, we propose estimators which, unlike in the previous studies, can be applied if EX and EY vanish, and even if the functions ψ and ϕ are not strictly positive. Our procedures involve estimating the functions ϕ and ψ, deduce from there predictors of X and Y , and construct nonparametric estimators of m using those predictors. We show that, under the restrictive assumptions made in the existing literature, this is relatively straightforward to do, whereas under the much weaker assumptions we also consider, we need to use a sophisticated approach.
This paper is organized as follows. We describe the covariate-adjusted model and discuss the model restrictions in the existing literature in Section 2. We propose several nonparametric estimators in Section 3, ranging from the most basic estimators which can be applied under similar restrictions as those imposed in the existing literature, to the most sophisticated ones which rely on much milder assumptions. We derive theoretical properties of our estimators in Section 4, where we show that they have the same first order asymptotic properties as the nonparametric estimators that could be computed if X and Y were observed directly. More surprisingly, in some particular cases, our new estimators can even achieve faster convergence rates than the standard estimators based on direct observations from (X, Y ). We discuss practical implementation of our methods in Section 5, where we also investigate their performance on simulated data, and apply them to analyze two real datasets studied in Şentürk and Müller (2005b) and Şentürk and Nguyen (2006). We discuss multivariate extensions in Section 6. Our proofs are provided in Section 7 and in a supplementary file.
2. Model and data. We observe independent and identically distributed (i.i.d.) triplets {(X i ,Ỹ i , U i )} n i=1 generated by the covariate-adjusted model of Şentürk and Müller (2005a), where for any pairs of random variables (Q, R) and (Q i , R i ), i = 1, . . . , n, let S Q,n (x; with w(s) = (1, s) T and where K is a kernel function, h = h n > 0 is a bandwidth and, for every t ∈ R, K h (t) = h −1 K(t/h).
If the (X i , Y i )'s were available, we could estimate m(x) nonparametrically by a standard local polynomial estimator constructed from the (X i , Y i )'s, the two most popular versions of which are the Nadaraya-Watson and the local linear estimators, defined bỹ ,m LL (x) = e T 1 S −1 X,n (x; K, h)T X,Y,n (x; K, h), (3.1) respectively. In our case, the (X i , Y i )'s are not observed and these standard estimators cannot be computed. We develop new nonparametric estimators that can be computed from the (X i ,Ỹ i , U i )'s, and whose complexity depends on whether (2.3)(a) and (b) are satisfied or not. The simplest situation is the one where (2.3)(a) holds. There, we can estimate m by standard nonparametric estimators based on predictors of the (X i , Y i )'s that are similar to, but less restrictive than, those used by Cui et al. (2009); see Section 3.2. The case where we do not assume (2.3)(a) requires more elaborate techniques: in Section 3.3, we suggest a method that can be used when (2.3)(b) is satisfied; we handle the most general case in Section 3.4, where we develop a sophisticated method which is valid regardless of whether (2.3)(a) and (b) hold or not. It involves computing estimators of unknown constant multiples of |ϕ| and |ψ|, estimate the zeros of those functions, construct piecewise estimators of unknown constant multiples of ϕ and ψ, estimate these constants and finally deduce estimators of ϕ and ψ.
3.2. Basic method. We start by deriving simple nonparametric estimators of m that can be computed when (2.3)(a) holds, and which form the basis of the more sophisticated methods we introduce in the subsequent sections. The idea is similar to the one used in the parametric context by i , whereφ andψ denote Nadaraya-Watson estimators of ϕ 0 and ψ 0 , divided by, respectively, EX = n −1 n i=1X i and EY = n −1 n i=1Ỹ i . It is because of this division that Cui at al. (2009) assume (2.3)(b), but the latter can be avoided and replaced by E|X|, E|Y | = 0 (which holds for all non-degenerate random variables), by better exploiting (2.3)(a). Specifically, under (2.3)(a), |ψ| = ψ, |ϕ| = ϕ, and 2) Motivated by this, we propose to estimate ψ and ϕ bŷ U,n (u; L, g 1 )T U,|X|,n (u; L, g 1 ) andψ + 0,LL (u) = e T 1 S −1 U,n (u; L, g 2 )T U,|Ỹ |,n (u; L, g 2 ) are local linear estimators of ϕ + 0 and ψ + 0 computed with a kernel function L and bandwidths g 1 and g 2 .
Then, we predict Y i and X i by takinĝ , we obtain the following estimators of m(x): (3.5) X,n (x; K, h)TX ,Ỹ ,n (x; K, h). Since they require predicting only the X i 's, these estimators seem more attractive than those in (3.5). However, it can be proved that their asymptotic "variance" is larger than that of the estimators in (3.5). Moreover, they cannot be adapted simply to the case where ϕ does not satisfy (2.3)(a); see Remark 3.3 in Section 3.3.
3.3. Refined procedure. As their parametric counterparts developed in the covariate-adjusted literature, the methods introduced in Section 3.2 can only be computed if (2.3)(a) holds. However, in practice, there is no reason why ϕ and ψ would always be positive, and even if they are, their estimators may vanish or get close to zero, which can cause numerical problems. In this section, we suggest a refined approach which can overcome these difficulties when (2.3)(b) holds. The more complex case where (2.3)(b) is violated will be dealt with in Section 3.4.
As in Section 3.2, to estimate m, the first step is to construct predictorŝ X i andŶ i , and thus estimators of ϕ and ψ. Recall the notation in (2.4).
of ϕ 0 and ψ 0 computed with a kernel function L and bandwidths g 1 and g 2 .
To derive consistent estimators of m without imposing (2.3)(a), recall that, for each i, X i and Y i are independent of U i . As a consequence, for any with ρ 1 , ρ 2 > 0 denoting two small numbers, the following modification ofm NW (x) at (3.1) would be consistent: and a similar consistent versionm LL (x; ρ 1 , ρ 2 ) ofm LL (x) at (3.1) could be constructed by replacing, in the definition ofm LL (x), sums over all i by sums over i ∈ C n (ρ 1 , ρ 2 ) as above. The advantage of this approach is that it enables us to exclude the data for which ψ(U i ) or ϕ(U i ) are small, and thus it can be applied even if (2.3)(a) does not hold.
Remark 3.2. While we shall prove in Section 4 that these estimators are consistent and have the same first order asymptotic properties as their counterparts at (3.1) based on undistorted data, in practice performance can be further improved by excluding a small fraction (say 5%) of the observations corresponding to the U i 's such that a kernel density estimator is the smallest. (Indeed, we know from standard properties of kernel regression estimators that, at points u where f U (u) is small,φ(u) andψ(u) are more variable.) Doing this corresponds to enlarging the set S slightly, which does not affect consistency and convergence rates, again due to the fact that the U i 's are independent of the (X i , Y i )'s, Remark 3.3. It is not possible to directly use this approach to modify the estimator discussed in Remark 3.1 for the case where ϕ has zeros, becausẽ Y i and U i are dependent. Particularly, we note that in general E(Ỹ i |X i = x, U i ∈ S) and E(Ỹ i |X i = x) are not equal.

3.4.
Elaborate procedure for the most general case. Finally we construct estimators of m that rely on neither part of (2.3). As before, we start by deriving predictors of the (X i , Y i )'s. Constructing predictorsX i (resp.,Ŷ i ) without assuming (2.3) requires to derive an estimator of ϕ (resp., ψ) without this assumption, which, unlike the methods used in the previous sections, turns out to be a challenging task. Our procedure is based on the fact that, , which implies that we can estimate ϕ * (resp., ψ * ) by a standard local linear estimatorφ * LL (resp.,ψ * LL ) with kernel L and bandwidth g 1 (resp., In what follows, we explain how to deduce an estimator of ϕ fromφ * LL . The same procedure can be applied to derive an estimator of ψ fromψ * LL . Since ϕ * is proportional to |ϕ|, to extract an estimator of ϕ fromφ * LL , we need to estimate the zeros of ϕ, say τ 1 , . . . , τ M for some finite M , at which ϕ changes sign. To do this we assume that, for each j, ϕ (τ j ) = 0. Then, it is straightforward to see that the first derivative of ϕ * has jump discontinuities at the τ j 's. Moreover, the zeros of ϕ coincide with those of ϕ * , so that, at the τ j 's, ϕ * reaches its minimum value, 0. Therefore, the τ j 's can be estimated using procedures for detecting discontinuities in derivatives of a regression curve, such as those in Gijbels et al. (1999) and Gijbels and Goderniaux (2005), combined with the fact that theτ j 's need to correspond to local minima ofφ * LL ; see Section 5.2 for details of implementation. For j = 1, . . . , M , letτ j denote the resulting estimator of τ j , and let I 0 = (−∞,τ 1 ), Our next target is to construct an estimator of ϕ. Recall the notation ϕ + 0 = ϕ · E|X| in (3.2). Recalling that ϕ changes sign at each τ j , we can obtain a consistent estimator of either ϕ + 0 or −ϕ + 0 (we'll see below how to distinguish these two cases) by takingφ + ±,0 (x) = M j=0 (−1) jφ * j,LL (x) · I(x ∈ I j ) , where, for each j,φ * j,LL denotes the local linear estimator of ϕ * constructed using only the (U i , |X i |)'s for which U i ∈ I j . Here we use a different local estimator in each I j because, under our assumptions, the first derivative of ϕ * = |ϕ| · E|X| is discontinuous at the τ j 's. It can be shown using standard kernel smoothing arguments that in this case the bias near the τ j 's is reduced by using this piecewise approach. Our next step is to extract fromφ + ±,0 an estimator of ϕ + 0 (recall that ϕ + ±,0 is an estimator of ϕ + 0 or −ϕ + 0 , but we can't know of which one). To do this, recall that E{ϕ(U )} = 1, which implies that E{ϕ + 0 (U )} > 0. This fact motivates us to estimate ϕ + 0 (x) byφ + 0 (x) =φ + ±,0 (x) sign n i=1φ + ±,0 (U i ) . Since ϕ + 0 (x) = ϕ(x) E|X|, once we have done this, to estimate ϕ it remains to construct an estimator of E|X|.
Noting that i . We can proceed similarly to construct predictorsŶ i of the Y i 's. As in Section 3.3, since, to obtain these predictors, we divide byφ(U i ) and ψ(U i ), when constructing our estimator of m we cannot use the (X i ,Ŷ i )'s for which |φ(U i )| or |ψ(U i )| is too small. Therefore, to estimate m we use the estimatorsm NW (x; ρ 1 , ρ 2 ) andm LL (x; ρ 1 , ρ 2 ) defined in Section 3.3, but with the predictorsX i andŶ i constructed above.
4. Theoretical properties. We start by establishing theoretical properties of the estimatorsm NW andm LL from Section 3.2. While these estimators seem intuitively natural, because they are computed using variables obtained through nonparametric prediction, checking whether they are consistent, and deriving detailed asymptotic properties, are quite difficult. Recently, Mammen et al. (2012) gave a deep account of nonparametric estimators computed from nonparametrically generated covariates, but our estimators do not fall into the class of settings they consider, not least because in our case, not only the covariate X, but also the dependent variable Y , are nonparametrically generated, which makes the problem even more complex than theirs. In addition to the basic model assumptions introduced in the first paragraph of Section 2, we make the following regularity assumptions: The bandwidths (h, g 1 , g 2 ) = (h n , g 1n , g 2n ) are such that h n −α 0 and g 1 n −β 1 and g 2 n −β 2 for some 0 < α 0 , β 1 , β 2 < 1/3.
Condition (B1) is a relaxed version of assumption (2.3) often assumed in the covariate-adjusted regression literature. See, for example, Şentürk andMüller (2005a, 2006) and Cui et al. (2009). Condition (B2) includes standard regularity and smoothness assumptions for the asymptotic results of kerneltype nonparametric regression estimation. In (B3), we relax the conventional boundedness condition on the covariates used by Şentürk andMüller (2005a, 2006) and Mammen et al. (2012), and assume instead that X has a finite exponential moment (for example this is satisfied if the distribution of X comes from the exponential family or is compactly supported). Condition (B4), which requires exponentially light tails of ε, is similar in spirit to Assumption 1. (iv) in Mammen et al. (2012). Like them, we need this technical assumption to employ an argument based on empirical processes. Condition (B5) is standard in the context of kernel regression, and is easy to satisfy since we can choose the kernels. Condition (B6) states the required range of magnitude of the bandwidths, and is easy to satisfy in practice.
The next two theorems establish uniform consistency and asymptotic normality of our estimatorsm NW andm LL defined in Section 3.2. Their proof can be found in Section 7 and in Section D in the supplementary file.
We deduce from the theorems that, although they are constructed from distorted data, when computed with appropriate bandwidths, our estimatorŝ m NW andm LL defined in Section 3.2 have the same uniform convergence rates as the standard estimators in (3.1) used when the (X i , Y i )'s are available. This contrasts with the errors-in-variables models studied by Fan and Truong (1993) and Delaigle et al. (2009), where convergence rates are significantly degraded by the measurement errors. The conclusions arising from the asymptotic distribution of our estimators are also interesting. Abusing terminology, we refer to V and B 0 +B (resp., B 1 +B) as the asymptotic variance and bias and of our estimatorm NW (resp.,m LL ), and we call asymptotic mean squared error (AMSE) the sum of the asymptotic variance and squared bias. We use similar terminology for the standard estimators of m.
We learn from part (ii) of both theorems that, if we choose g 1 and g 2 of order o(h), the asymptotic bias and variance of our estimators are identical to those of standard estimators, and there, as in the standard case, it is optimal to take h n −1/5 , so that AMSE n −4/5 . Perhaps more surprisingly, in cases where B 0 (resp., B 1 form LL ), B ϕ and B ψ do not all have the same sign, it is possible to choose h and g 1 or g 2 an order of magnitude slightly larger than n −1/5 such that the asymptotic bias B 0 +B (resp., B 1 +B) vanishes and the AMSE our estimator is of order o(n −4/5 ), thus smaller than the AMSE of the standard estimator (similar results can be established for the integrated AMSE). However, while it is theoretically interesting, we were not able to exploit this result in practice to make our estimator outperform the standard one, despite several attempts. In part this is because to benefit from this result we need to choose the bandwidths in a very specialized way that requires estimating too many unknowns, and we found that the simpler bandwidths choice suggested in Section 5.2 almost always worked better.
Next, we develop theoretical properties of our estimator defined in Section 3.3. We start by rewritingĈ n (ρ 1 , ρ 2 ) asĈ n (ρ 1 , We can rewrite the estimator at (3.7) aŝ To emphasize the main idea while avoiding repetitive arguments, here we present the theoretical result only for this estimator, assuming that only ϕ may have zeros, and therefore we take ρ 2 = 0 throughout this section. A straightforward adaptation of the arguments used to prove Theorem 4.3 below leads to similar results in the more general case where ϕ has zeros and ρ 2 > 0, and for the local linear estimatorm LL (x; ρ 1 , ρ 2 ).
When ρ 2 = 0,L n (ρ 1 , ρ 2 ) depends only on ρ 1 ; to simplify notation we rewrite it asL n (ρ 1 We will need the following assumptions: consists of finitely many points located in the interior of I U and min u∈∂L(t) |ϕ (u)| > 0 is non-empty.
The next theorem establishes uniform consistency and asymptotic normality ofm NW (x; ρ). See Section E in the supplementary file for its proof.
We deduce from the theorem that our estimator defined in Section 3.3 has the same uniform convergence rate as the standard Nadaraya-Watson estimator in (3.1), used when the data (X i , Y i ) are available. Moreover, as long as we choose g 1 and g 2 of order o(h), the asymptotic "bias" and "variance" of our estimator from Section 3.3 are equal to those of the standard Nadaraya-Watson estimator, where i ∈ {1 ≤ j ≤ n : U j ∈ L n (ρ)}. As we already indicated below Theorems 4.1 and 4.2, in theory in some cases it is possible to choose the bandwidths in such a way that the AMSE of our estimator tends to zero faster than that of the standard estimator, but it seems very hard to find a way to exploit this in practice. Similar results can be established for the local linear estimatorm LL (x; ρ 1 , ρ 2 ).
Establishing theoretical results for the more general procedure described in Section 3.4 is particularly challenging. Recall that this method combines a change point detection algorithm and the ridge-parameter based method introduced in Section 3.3. The complex nature of this approach implies that deriving its theoretical properties rigorously requires long and tedious arguments. Since our paper is already very long, and even the proofs for our simpler methods are fairly tedious, we leave such rigorous derivations for future work. However, our preliminary calculations already indicate that the procedure from Section 3.4 should have asymptotic properties similar to those described in Theorem 4.3. In particular, these calculations indicate that estimating the τ j 's and the sign of ϕ and/or ψ has no first order asymptotic effect on the properties of our estimators of m.

Numerical results.
5.1. Which method to use. The approach in Section 3.4 can be applied in essentially all cases, but since the methods from Sections 3.2 and 3.3 are simpler, the user might prefer to use these if all parts of (2.3) hold. While (2.3) can be verified by standard tests of hypothesis applied to the observed data (see Remark 5.1 below), when these conditions are needed, it is because the techniques employed involve dividing by estimators of ψ, ϕ, EX or EY . Therefore, in practice, to avoid numerical issues, we suggest using the method from Section 3.3, and to use instead the method from Section 3.4 if the absolute values of estimators of EX or EY are small, the extent of which depends on the magnitude of other quantities involved and the precision of the software employed. This is generally rather easy to determine by examining the data, but if unsure the user can just apply the method of Section 3.4, which is valid in the most general case.
We note too that one does not necessarily need to predict the X i 's and the Y i 's with the same method. For example, if one is confident that EX is far from zero, but is not sure about EY , then the X i 's could be predicted using the approach from Section 3.3, and the predictors of the Y i 's could be imsart-aos ver. 2013/03/06 file: Cov_Adj_Main.tex date: January 13,2016 obtained from the approach suggested in Section 3.4.
Remark 5.1. The assumption at (2.3) can be tested in several ways. For example, since EX = EX, we can first test the sign of EX by a standard test of hypothesis for the mean applied to the dataX 1 , . . . ,X n , and then test the sign of the function ϕ 0 = ϕ · EX at (2.4), using for example tests such as those in Dümbgen and Spokoiny (2001), Chetverikov (2012) and Lee et al. (2013), applied to the observed data. 5.2. Details of implementation. As in the case where the (X i , Y i )'s are available, in practice we recommend using the local linear versions of our estimators, and in this section we suggest ways of choosing the parameters required to compute them. Similar ideas can be used for the Nadaraya-Watson estimators. We know from Section 4 that, while we have to choose h with care, we have more flexibility for the bandwidths g 1 and g 2 , which can take a large range of values. If we take h to be of the standard size for nonparametric regression, and g 1 = o(h) and g 2 = o(h), then our estimators have the same first order asymptotic properties as the estimators at (3.1).
Motivated by this, for the estimators in Section 3.2, we take g 1 = n −0.1 g 1,PI , g 2 = n −0.1 g 2,PI and h = h PI , where the subscript PI means that we use a standard plug-in bandwidth for local linear estimators (Ruppert et al., 1995) constructed based on, respectively, the data (U i , |X i |), (U i , |Ỹ i |) and (X i ,Ŷ i ). For the estimators in Sections 3.3 and 3.4, we take g 1 = n −1/10 g 1,PI and g 2 = n −1/10 g 2,PI , where g 1,PI and g 2,PI denote standard plug-in bandwidths for local linear estimators constructed based on, respectively, the data (U i ,X i ) and (U i ,Ỹ i ). Then, in Section 3.3, we choose ρ 1 = max(0.1, ρ * 1 ) and ρ 2 = max(0.1, ρ * 2 ), where ρ * 1 (resp., ρ * 2 ) denotes the square root of an estimator of the asymptotic "mean squared error" ofφ LL (resp.,ψ LL ), integrated over the set of x-values where |φ LL (x)| (resp., |ψ LL (x)|) take its smallest values; see Appendix A in the supplementary file for details. We do the same for the method from Section 3.4, except that we use the estimatorsφ andψ of ϕ and ψ derived there. Finally, we take h = h PI , a standard plug-in bandwidth for local linear estimators computed from the data (X i ,Ŷ i ), i ∈Ĉ n (ρ 1 , ρ 2 ).
The estimators from Section 3.4 also require to estimate the zeros τ 1 , . . . , τ M at which ϕ changes sign, and the same is required for ψ if the method in that section is used to compute predictors of the Y i 's. We proceed as follows. First, since the τ j 's all correspond to a local minimum of ϕ * , we find all the points at whichφ * LL has local minima. Then, among those points we keep only those which are close to the discontinuity points of the derivative ϕ * detected by the method of Gijbels and Goderniaux (2005). Here we define "close" by less than 2h away, where h is the bandwidth in Section 2.2.1 of Gijbels and Goderniaux (2005). Finally, to slightly improve numerical performance, we implement Remark 3.2 and remove the data corresponding to the 5% smallestf U (U i )'s.

5.3.
Simulations. We applied our methods to a variety of simulated examples, ranging from the simplest ones in which ψ > 0 and ϕ > 0, where we can use the method from Section 3.2, to more complex ones in which EX = 0 and both ψ and ϕ oscillate between positive and negative values, where we need to use the sophisticated approach suggested in Section 3.4.
We generated data (X i ,Ỹ i , U i ), i = 1, . . . , n, from model (2.1) for n = 100, 200, 500 and 1000, and considered various combinations of m, ϕ, ψ and σ, and various distributions of X i and U i . We took ε i ∼ N (0, 1), and considered shifted versions of three regression curves m, denoted by m 1 , m 2 and m 3 and defined as m where φ µ,θ denotes the density of a N (µ, θ 2 ). In all cases below, the generic constant const. was chosen so that E{ϕ(U )} = E{ψ(U )} = 1.
First, we considered models where the local linear estimators from Sections 3.2 to 3.4 could all be applied:  Next, we considered models (i.c)-(iii.c) and (i.d)-(iii.d), where we took m, X i and σ as in models (i.a)-(iii.a) and (i.b)-(iii.b), respectively, but took U i ∼ β(3, 5) and ϕ(·) = ψ(·) = const. m 1 (5 · −2). Here ϕ and ψ have zeros and change signs, so that the method from Section 3.2 cannot be applied. Finally, in our last models, ϕ and ψ change signs and have several zeros and E(X i ) = 0, so that we can apply only the method from Section 3.4: (iv.a) m(·) = m 1 (· + 1), X i ∼ N (0, 1.5 2 ), σ(x) = 0.3; (v.a) m(·) = m 2 (· + 1),  We compared each of our estimators with the ideal estimatorm LL at (3.1) computed from the (X i , Y i )'s, which are not available in real data applications but are available when we simulate data, and with the inconsistent naive estimatorm LL,naive , which is the standard local linear estimator computed from the contaminated (X i ,Ỹ i )'s. For each n and each model, we generated 1000 samples and constructed each estimator for each sample. Letm denote any one of the estimators considered below. To summarise the performance ofm, we computed, for each sample, the integrated squared error ISE = b a {m(x) − m(x)} 2 dx, where, in each case, a and b were the quantiles 0.025 and 0.975 of the distribution of X.
In Tables 1 to 4 in Appendix B in the supplementary file, for each method we report the first, second and third quartiles of the resulting 1000 ISEs. See Appendix B for a detailed discussion of the simulation results. In summary, we found that, as expected, when ϕ, ψ, EX and EY were different from zero, but EX and/or EY were relatively close to zero, the estimator that worked best was the one from Section 3.2, but the most complex estimator from Section 3.4 worked well. When EX and EY were far from zero, all three estimators worked well, with the simplest one from Section 3.2 giving the best results and the one from Section 3.4 working the worst. When ϕ and/or ψ had zeros, the estimator from Section 3.2 could not be applied, and when EX and EY were close to zero, the best results were obtained with the estimator from Section 3.4, whereas when EX and EY were far from zero, the estimator from Section 3.3 worked best. Finally, we found that our approach also performed well when the errors were heteroscedastic.
In all cases, our estimators performed considerably better than the naive estimator, but were of course outperformed by the oracle estimator. As expected, the performance of our estimators improved as sample size increased. In all our simulation settings, the estimator from Section 3.4 gave reasonable results. However, if ϕ and ψ were far from zero, we got better results by using the simplest estimator from Section 3.2, and if EX and EY were far from zero, we got better results using the estimator from Section 3.3.
To illustrate these results graphically, we present a few figures that are representative of the conclusions of our simulations. For each estimatorm presented in the figures, we show the three estimated curves corresponding to the first three quartiles of the 1000 ISEs defined above. In Figure 1, using example (i.a), we illustrate the fact that, when all three methods can be applied, they often give similar results. Figure 2 shows estimated curves for examples (ii.c) and (ii.d). We can see that, in case (ii.c), where EX is close to zero, the estimatorm LL (·; ρ 1 , ρ 2 ) from Section 3.4 worked better than the one from Section 3.3, but that the reverse is true in case (ii.d), where EX and EY are both far from zero. In that figure, we also depict the naive estimatorm LL,naive , which performed very poorly. Finally, in Figure 3, we use example (vi.a) to demonstrate the improvement that our estimator m LL (·; ρ 1 , ρ 2 ) from Section 3.4 benefits from as the sample size n increases. Here too, the naive estimator performed very poorly, even for n large.

5.4.
Real data illustrations. We applied our new method to the Boston house-price dataset described in Harrison and Rubinfeld (1978), available at http://lib.stat.cmu.edu/datasets, and which contains information about houses and their owners at 506 locations around Boston. As in Şentürk and Müller (2005b), we are interested in the relationship between the median price (in USD 1000's) of houses,Ỹ , and per capita crime rate by town, X, with the confounding effect of the proportion of population of lower educational status, U , removed. Şentürk and Müller's (2005b), whose interest was in the correlation betweenX andỸ , concluded that this correlation alters dramatically after adjusting for the confounding effect of lower educational status. On the left panel of Figure 4, we depict the covariate-adjusted regression curve obtained using the local linear estimatorm LL (·; ρ 1 , ρ 2 ) from Section 3.3, the estimatorm LL (·; ρ 1 , ρ 2 ) from Section 3.4, and the naive regression estimatorm LL,naive obtained by regressingỸ onX after removing a few outliers. In this example, the estimator from Section 3.2 was identical to the one from Section 3.3. We can see thatm LL,naive indicates a pronounced relationship between house price and crime rate (as crime rate increases, house price decreases), but once we adjust for the effect of lower educational status, the regression curve obtained by both versions of our estimator is almost flat, indicating a weak relationship between the adjusted X and Y . Next, we applied our procedure to the diabetes dataset used by Schorling et al. (1997) and Willems et al. (1997), available at http://biostat.mc.van derbilt.edu/DataSets, which represents a subset of 403 individuals taken from a larger cohort of 1046 subjects who participated in a study for African Americans about obesity, diabetes and related factors in central Virginia. As in Şentürk and Nguyen (2006), our goal was to examine the relationship between glycosolated hemoglobin levelỸ , a biomarker for diabetes, and diastolic blood pressureX, adjusting for the effect of body mass index, U , which was found to be a confounder for both variables. As in Şentürk and Nguyen (2006), we removed a few outliers before our analysis. As in the previous example, EX and EY were far from zero, so that we used the estimatorm LL (·; ρ 1 , ρ 2 ) from Section 3.3, which we compared with the naive estimatorm LL,naive . Here too, the estimator from Section 3.2 was identical to the one from Section 3.3. We also computed the estimator fromm LL (·; ρ 1 , ρ 2 ) from Section 3.4. These estimators, depicted on the right panel of Figure 4, show that after adjusting for body mass index, the relationship between glycosolated hemoglobin level and diastolic blood pressure is noticeably less pronounced. We should highlight that, in this example, the data were rather sparse for diastolic blood pressure greater than 100, and the few patients for whichX was greater than 100 had a rather low value ofỸ , whence the decreasing shape on the right hand side of the graph, which may just be an artifact of the sparseness of the data in that area.
Another interesting application of our method is to the baseline data collected from studies A and B of the Modification of Diet in Renal Disease Study (Levey et al., 1994). The nonlinear relationship between the baseline unadjusted glomerular filtration rate (GFR) and serum creatinine (SCr) is of particular interest. Taking body surface area (BSA) as the confounder, Cui et al. (2009) used a parametric nonlinear model of the form m(x) = β 1 exp(−β 2 − β 3 x 2 ) + β 4 to study the relationship between GFR and SCr after correcting for the distorting effect of BSA. Because this dataset is not publicly available, we shall not compare the proposed nonparametric method with that of Cui et al. (2009) in this paper.
6. Generalizations to the multivariate case. Our approach can be generalized to the d-variate case, d ≥ 1, where we observe data distributed like a vector (U,X T ,Ỹ ), withX ∈ R d a distorted version of X ∈ R d . Reflecting the fact that the components ofX may not all be distorted, we write d = d 1 + d 2 , with d 1 ≥ 0 and d 2 ≥ 1, and let X = (X T 1 , X T 2 ) T and . . , X d ) T , and where we use the convention that X = X 2 if d 1 = 0. In this notation, the data we observe are generated by the model is a curve we wish to estimate, the random variables X, U and ε are mutually independent, E(ε) = 0 and var(ε) = 1. As in (2.2), we assume that E{ψ(U )} = 1, E{ϕ r (U )} = 1, for r = 1, . . . , d 2 .
The procedures from Section 3.2 to 3.4 can each be generalized to the multivariate setting, but for space constraint here we show only how to generalize the approach from Section 3.2. The same ideas can be applied for the methods from Sections 3.3 and 3.4. To construct a nonparametric version of the estimator from Section 3.2, we first construct predictorsŶ i andX i,d 1 +1 , . . . ,X id as in equation (3.4), and let X i = (X i1 , . . . , X id 1 ,X i,d 1 +1 , . . . ,X id ) T .
Next, we use a standard multivariate local linear regression estimator applied to the data (X T i ,Ŷ i ). That is, we define (see Fan and Gijbels, 1996)  It is well known that fully nonparametric estimators suffer from the curse of dimensionality, which means that as d increases, such estimators can only work reasonably well if the sample size is very large. To overcome this problem, a common approach is to restrict the regression model so that only univariate curves have to be fitted. A popular example is the additive model (Hastie and Tibshirani, 1990), which assumes that m(X) = m 0 + d j=1 m j (X j ). In our context, the additive covariate-adjusted regression model can be written as where m 1 , . . . , m d are unknown univariate functions satisfying E{m j (X j )} = 0 for j = 1, . . . , d and m 0 is an unknown parameter.
In the standard setting where the (X T i , Y i )'s are directly observed, there are several ways to fit the additive model; see Horowitz (2014) for an overview of estimation and inference for nonparametric additive models. The simplest approach is to adapt to our setting the iterative backfitting algorithm of Buja et al. (1989), as follows. First, letm 0 = n −1 n i=1Ỹ i andm j ≡ 0 for j = 1, . . . , d. For j = 1, . . . , d, updatem j by taking it equal to a local linear regression estimator using the data Alternatively, instead of takingm j = 0 as initial estimators, we could start with a linear approximation of the model in (6.2). See Appendix C in the supplementary file for details. We could also apply similar transformations to other existing methods for fitting additive models, such as the approach suggested by Horowitz and Mammen (2004). The main theoretical challenge is a delicate analysis on how the presence of generated response and predictors affects the first order asymptotic properties of the final estimators. However, deriving such results requires much more work than can possibly done in this paper, and so we leave this problem for future research. The method proposed in this section can be applied to creatinine data, which was analyzed by Şentürk and Müller (2006). In this study, serum creatinine level is taken as the response and the two predictors include cholesterol level and serum albumin level. The confounder variable U is taken to be body mass index defined as weight/height 2 . The readers can find more details about this dataset in Şentürk and Müller (2006). 7. Proof of Theorem 4.1. We start by introducing basic notations. For a kernel function K, we write µ K, = u K(u) du for non-negative integers . For any set S, we denote its complement by S c and its cardinality by #S. Throughout, we let const. denote a finite positive constant independent of n, which may take different values at each occurrence. We also use the following notation: µ 0 = E(X), m 0 = E(Y ), µ + 0 = E|X|, m + 0 = E|Y | and We proceed with the proof of Theorem 4.1. For u ∈ I U = [0, 1], write andψ + 0,LL are local linear estimators of ϕ + 0 and ψ + 0 defined below (3.3). Noting the model at (2.1), and hence by (3.4) and (7.2), Substituting the expressions in (7.3) into (3.5) giveŝ Proof of (i).
Let α ∈ (α 0 , λ 1 ) be a constant, such that under Condition (B6), n −α = o(h) and P {E n (α) c } → 0 as n → ∞. On the event E n (α) ∩ E 1n (C 1 ), we have (7.16) such that for every x ∈ [a, b], |X i − x| ≤ |X i − x| + C 1 n −α log n. Therefore, on the event E n (α) ∩ E 1n (C 1 ) with n sufficiently large, It follows from (7.13) and (7.17) that, on E n (α)∩E 1n (C 1 ) with n large enough, denotes the empirical distribution function. To further bound the right-hand side of (7.18), we let F X be the distribution function of X and then apply the Dvoretzky-Kiefer-Wolfwitz inequality (Massart, 1990) to obtain that P ( Substituting this into (7.18) and taking A n ≡ E n (α) ∩ E 0n (1/2) ∩ E 1n (C 1 ) imply that for all sufficiently large n, and that P (A c n ) → 0 as n → ∞. Together, (7.9) and (7.21) prove (7.11). Next we studyΠ 01 (x). For this, we first writefX (x)Π 01 (x) as leading to (7.11) can be used to prove that max x∈ [a,b] (1) and the same bound holds if the m(X i )'s are replaced by the σ(X i )'s. Moreover, similarly to (F.9) in the proof of Lemma F.1, it can be proved that max 1≤i≤n |ε i | = O P (log n). (7.24) This, together with (7.9) and the two displays before (7.24) yields max x∈[a,b] Here, the last step follows from Condition (B6) and the assumption that α 0 < 2β 1 . Together, (7.12), (7.22), (7.23) and (7.25) Together with (7.7) and (7.12), this implies max x∈ [a,b] Π 02 (x) − Π 02 (x) = O P {δ n,0 (g 1 )}, (7.28) where Eg n,i (x), To bound max x∈[a,b] |R n (x)|, we create a grid using N points of the form x j = a + j with = (b − a)/N for some N ≥ 1 to be determined below (7.32). Since g n,i ( by the mean value theorem we have, for every x, y ∈ R, |g n,i (x) − g n,i (y)| ≤ ( For each x ∈ R fixed, g n,1 (x), . . . , g n,n (x) are independent random variables satisfying |g n,i (x)| ≤ K ∞ m ∞ and E{g n, Hence, by Bernstein's inequality and Boole's inequality, for every y ≥ 0, 32) such that in view of (7.31), P {C n (λ) c } ≤ 2N exp(−τ λ) for some absolute constant τ > 0. By taking N = n and λ = 2τ −1 log n, it follows from (7.30) and (7.32) that max x∈ [a,b] |R n (x)| = O P {h 1/2 (n/ log n) −1/2 + n −1 log n + (nh) −1 }. (7.33) For the second term on the right-hand side of (7.29), standard arguments show that, under Conditions (B3) and (B5), uniformly in x ∈ [a, b]. Consequently, combining (7.12), (7.28), (7.33) and (7.34), we get max For the last termΠ 03 (x) in (7.4), we need to control the stochastic error (7.3). To this end, we shall use a lattice argument by making a finite approximation of the compact interval [a, b] using a sequence {x j } N j=1 of equidistant points x j = a+j for = (b−a)/N , and then discretize ∆ n,∞ to define ∆ n, Here, N is a positive integer that will be determined after (7.44).
(7.44) Combining (7.43) with N = n, λ = 2c −1 log n, (7.44) and the fact that P (A c n ) → 0 proves (7.38) as claimed. Proof of (ii). To prove the asymptotic normality, we need to use a more refined argument. In what follows, x ∈ [a, b] is fixed and we deal with the sum in (7.4) over eachΠ 0j (x) separately.