A non-parametric ANOVA-type test for regression curves based on characteristic functions

A nonparametric ANOVA-type test for regression curves based on characteristic functions J.C. Pardo-Fernández1∗, M.D. Jiménez-Gamero2†, A. El Ghouch3‡ 1Departamento de Estatística e I.O., Universidade de Vigo, Spain 2Departamento de Estadística e I.O., Universidad de Sevilla, Spain 3ISBA, Université catholique de Louvain, Belgium Abstract This article studies a new procedure to test for the equality of kregression curves in a fully nonparametric context. The test is based on the comparison of empirical estimators of the characteristic functions of the regression residuals in each population. The asymptotic behaviour of the test statistic is studied in detail. It is shown that under the null hypothesis the distribution of the test statistic converges to a finite combination of independent χ2 1random variables. The coefficients in this linear combination can be consistently estimated. The proposed test is able to detect contiguous alternatives converging to the null at the rate n−1/2. The practical performance of the test based on the asymptotic null distribution is investigated by means of simulations. Key Words: Comparison of regression curves; empirical characteristic function; regression residuals. ∗Corresponding author. Departamento de Estatística e I.O., Facultade de Ciencias Económicas e Empresariais, Campus Universitario As Lagoas-Marcosende, 36310 Vigo, Spain. e-mail: [email protected]. Telephone: +34 986 813505. Fax: +34 986 812 401. †Departamento de Estadística e I.O., Facultad de Matemáticas, Calle Tarfia s.n., 41012 Sevilla, Spain. e-mail: [email protected]. ‡Institut de statistique, biostatistique et sciences actuarielles (ISBA), Université catholique de Louvain, Voie du Roman Pays 20, B-1348 Louvain-la-Neuve, Belgium. e-mail: [email protected]. 1 1 Introduction Testing for the equality of the means of kpopulations (k≥2) is a classical problem in Statistics. When the populations are assumed to follow a normal distribution with equal variance, then the ANOVA F-test is the classical way to perform the test. In this paper we consider a more general setting. We assume that in each population along with the response variable, Y, we observe another variable, X, the covariate, so that the mean and the variance of the response variable depend on the values of the covariate. More specifically, let (Xj, Yj),1≤j≤k, be kindependent random vectors satisfying general nonparametric regression models Yj=mj(Xj) + σj(Xj)εj,(1) where mj(x) = E(Yj|Xj=x)is the regression function, σ2 j(x) = V ar(Yj|Xj=x)is the conditional variance function and εjis the regression error, which is assumed to be independent of Xj. Note that, by construction, E(εj) = 0 and V ar(εj)=1. The covariate Xjis continuous with density function fj. The regression functions, the variance functions, the distribution of the errors and the distribution of the covariates are unknown and no parametric models are assumed for them. Under this framework our approach is fully nonparametric. In this conditional setting, the hypothesis of equality of means is stated in terms of the conditional means or regression functions H0:m1=m2=. . . =mk, or, in other words, the mean effect of the covariates over the responses is equal in the k populations. Since the objective is to compare the regression curves, it is reasonable to assume that the covariates have common support. The alternative hypothesis is H1:H0is not true. Note that this testing problem contains the simpler case described in the first paragraph as a particular case by only eliminating the covariates in the models. The problem of testing for the equality of regression curves in nonparametric settings has been previously treated in the statistical literature. The majority of the available papers are either devoted to the comparison of only two curves or impose some restrictive assumptions like fixed design, equal sample sizes, identical design points or homoscedasticity of the residuals. Examples of such works include Delgado (1993), Kulasekera (1995), 2 Munk and Dette (1998), Neumeyer and Pardo-Fernández (2009) and Srihera and Stute (2010), among many others. To the best of our knowledge, the most related works to ours are the papers of Neumeyer and Dette (2003) and of Pardo-Fernández et al. (2007). These papers proposed and studied a procedure to test the hypothesis of equality of k (k≥2) regression functions based on the comparison of marked empirical processes of the residuals for the former and the comparison of distribution functions of the residuals for the latter. In this work, we investigate a new test procedure based on a weighted process generated by the characteristic functions of the residuals. Compared to the competitors, the main advantages of the proposed method are the following. First, the test is universally consistent for any fixed alternative without any restrictions on the weight function involved in the definition of the test statistic (except for the fact that it should be positive) and without any restrictions on the distribution of the regression errors. Thus, the method can be applied when the distributions of the errors are arbitrary: continuous, discrete or mixed. This is illustrated in the simulation study where an example with errors having a mixed distribution is given. Second, the test can be used to check the equality of any number of regression curves and can detect local alternatives converging to the null at the rate n−1/2. Third, the critical values or the p-values can be obtained from the asymptotic null distribution of the proposed test statistic. This is not the case for most of the existing methods, which typically rely on bootstrap. In our case, although we study a bootstrap version of our test statistic in order to perform a comparison with other methods, the bootstrap is not needed because we are able to get the desired level and a good power even for moderate sample sizes by using the asymptotic null distribution. This is clearly shown for many scenarios considered in the simulation study. Our approach is based on comparing the characteristic functions of the regression errors. More specifically, let εj={Yj−mj(Xj)}/σj(Xj)be the regression error in population j. Let m0be the common regression curve under the null hypothesis, and define ε0j={Yj−m0(Xj)}/σj(Xj) = εj+{mj(Xj)−m0(Xj)}/σj(Xj),(2) 1≤j≤k. It turns out that the null hypothesis H0is true if and only if, for all 1≤j≤k, the random variables εjand ε0jhave the same distribution (see Theorem 1 in Pardo- Fernández et al., 2007). This assessment can be interpreted in terms of the cumulative distribution function (cdf) or in terms of any other function characterizing the probability law of the errors. Pardo-Fernández et al. (2007) restricted their attention to the cdf. The probability law of any random variable Xis also characterized by its characteristic function (cf), ϕ(t) = E{exp(itX)}. Recent years have witnessed an increasing number 3 of proposals for hypothesis testing whose test statistics measure deviations between the empirical characteristic function (ecf) of the available data and an estimator of the cf under the null hypothesis. In the line of the setting considered in this paper, that is, by assuming that the data are generated by regression models, are the papers by Jiménez- Gamero et al. (2005) and Hušková and Meintanis (2007, 2010), for testing goodness of fit (gof) for the errors, and Hušková and Meintanis (2009) for testing gof of the regression function to a parametric function. An advantage of the cf approach over the one based on the cdf, as observed in Hušková and Meintanis (2009), is that the former usually requires less stringent assumptions for its validity. In addition, from the simulation results that can be found in the literature (see, for example, Hušková and Meintanis, 2010) the tests based on the ecf compete very satisfactorily with those based on the empirical cdf (ecdf). Having in mind the reasons above, the purpose of the present paper is to test H0by comparing consistent estimators of the cfs of the random variables εjand ε0j, say ˆϕj(t) and ˆϕ0j(t), respectively, 1≤j≤k. To measure deviations between these estimators, we consider a Cramér-von Mises type test statistic. In order to derive the asymptotic null distribution of the test statistic, we first give a stochastic expansion for the differences ˆϕj(t)−ˆϕ0j(t). As a consequence of this expansion, it is shown that the test statistic converges in law to a finite linear combination of independent chi-square variables. Moreover, under certain weak conditions on the distributions of the errors and the covariates, the asymptotic null distribution is proportional to a χ2 k−1distribution. We provide consistent estimators of the coefficients appearing in this linear combination, which allow us to propose a consistent estimator of the asymptotic null distribution. The behaviour of the test under fixed and local alternatives is also studied. Specifically, it is shown that the proposed test is able to detect any fixed alternative and contiguous alternatives converging to the null at the rate n−1/2, where ndenotes the total sample size. In contrast to the procedure in Pardo-Fernández et al. (2007), to derive these properties we do not require further assumptions on the distribution on the errors, such as that they have a probability density. For finite sample, the quality of the proposed approximation of the null distribution of the test statistic as well as the power are investigated numerically. From this numerical study, we conclude that the proposed approximation of the null distribution works adequately for moderate sample sizes, and in terms of power the proposed test competes very satisfactorily with those based on the ecdf. The paper is organized as follows. Section 2 introduces the test statistic and also provides an alternative expression, which is useful from a computational point of view. The asymptotic null distribution of the test statistic and the behaviour of the test under fixed and local alternatives are studied in Section 3. Section 4 reports a summary of a 4 numerical study conducted to study the practical performance of the test and to compare it with other existing methods. Section 5 concludes the paper. All proofs of the theoretical results are deferred to the Appendix. The following notation will be used along the paper: P0denotes probability assuming that H0is true; E0denotes expectation assuming that H0is true; P∗denotes the conditional probability law, given the data; all limits in this paper are taken when n→ ∞; L →denotes convergence in distribution; P →denotes convergence in probability; a.s. →denotes the almost sure convergence; if x∈Rk, with x0= (x1, . . . , xk), then diag(x)is the k×k diagonal matrix whose (i, i)entry is xi,1≤i≤k; for any complex number z=a+ib, Re(z) = ais its real part, ¯z=a−ibis its conjugate and |z|is its modulus; Nk(µ, Σ) denotes the multivariate normal distribution with mean vector µand variance-covariance matrix Σ; an unspecified integral denotes integration over the whole real line R; for a given non-negative real-valued function wand for any complex-valued measurable function g, we denote kgkw=R|g(t)|2w(t)dt1/2to the norm in the Hilbert space L2(R, w). 2 The test statistic Let (Xjl, Yjl),1≤l≤nj, be independent and identically distributed (iid) observations from (Xj, Yj),1≤j≤k. Let fj(x)be the probability density function (pdf) of Xj,n= Pk j=1 nj, and let fmix(x) = Pk j=1 pjfj(x)be the pdf of the mixture of covariates according to the weights p1, . . . , pk, where pj= lim nj/n. In order to estimate the errors, we first need to estimate the regression functions, mj(x) = E(Yj|Xj=x), the variance functions, σ2 j(x) = E[{Yj−mj(x)}2|Xj=x], and the common regression function under H0,m0(x) = Pk j=1 pj{fj(x)/fmix(x)}mj(x). With this aim we use nonparametric estimators based on kernel smoothing techniques. Let Kdenote a nonnegative kernel function defined on R, let 0< hn≡h→0be the bandwidth or smoothing parameter and Kh(x) = h−1K(x/h). We use the following estimators for the functions mj,σ2 jand m0: ˆmj(x) = nj X l=1 wjl(x)Yjl,ˆσ2 j(x) = nj X l=1 wjl(x)Y2 jl −ˆm2 j(x),ˆm0(x) = k X j=1 nj n ˆ fj(x) ˆ fmix(x)ˆmj(x), where ˆ fj(x) = n−1 j nj X l=1 Kh(x−Xjl),ˆ fmix(x) = k X j=1 nj nˆ fj(x), 1≤j≤k. The quantities wjl are, either the local-linear weights given by wjl(x) = Kh(Xjl −x)S2,nj(x)−(Xjl −x)S1,nj S0,nj(x)S2,nj(x)−S2 1,nj(x), 5 with Sk,nj(x) = Pnj l=1(Xjl −x)kKh(Xjl −x),k= 0,1,2, or the Nadaraya-Watson weights wjl(x) = Kh(Xjl −x) Pnj v=1 Kh(Xjv −x), Both are particular cases of local-polynomial weighting (see Fan and Gijbles, 1996). Under the model assumptions that will be stated in the next section, the results in this article are valid for local-linear and for Nadaraya-Watson (local-constant) estimators. Based on these estimators, for each population j,1≤j≤k, we construct two samples of residuals: ˆεjl =Yjl −ˆmj(Xjl) ˆσj(Xjl)and ˆε0jl =Yjl −ˆm0(Xjl) ˆσj(Xjl),(3) 1≤l≤nj, whose ecfs are ˆϕj(t) = 1 nj nj X l=1 exp(itˆεjl)and ˆϕ0j(t) = 1 nj nj X l=1 exp(itˆε0jl),(4) respectively. These ecfs are nothing but (consistent) kernel based nonparametric estimators of the population cfs ϕj(t) = E{exp(itεj)}and ϕ0j(t) = E{exp(itε0j)}, respectively, where ε0jis as defined in (2). The testing procedure consists of comparing ˆϕj(t)and ˆϕ0j(t),1≤j≤k, using a weighted L2−distance. More precisely, following the work of Hušková and Meintanis (2007, 2009, 2010) we define the test statistic T1n≡T1n(w) = k X j=1 nj nkˆϕj(t)−ˆϕ0j(t)k2 w,(5) where wis any given non-negative weight function with finite integral, Rw(t)dt < ∞. The presence of the weight function win the integrals appearing in the expression of T1n is necessary in order to ensure their finiteness, since kˆϕj(t)−ˆϕ0j(t)k2 w≤4Rw(t)dt, for all j. The motivation behind the test statistic T1nis the following: T1nconverges in probability to (see Theorem 7 below) T1≡T1(w) = k X j=1 pjkϕj(t)−ϕ0j(t)k2 w.(6) Under H0,ϕj(t) = ϕ0j(t)for all tand for 1≤j≤k, and thus T1vanishes. As a consequence, under H0,T1nshould be “very small”. We then conclude that, any value of T1nwhich is “significantly large” should lead to the rejection of H0. In practice, given a 6 significance level, a threshold value above which H0is rejected needs to be established. To this end we need to study the null distribution of T1n. Since this distribution is unknown, as an approximation to it we derive the asymptotic null distribution. This will be done in the next section. Remark 1 From Lemma 1 in Alba-Fernández et al. (2008), an alternative expression for T1n, which is useful from a computational point of view, is given by nT1n= k X j=1 1 nj(nj X l,s=1 Iw(ˆεjl −ˆεjs) + nj X l,s=1 Iw(ˆε0jl −ˆε0js)−2 nj X l,s=1 Iw(ˆεjl −ˆε0js)), where Iw(t) = Rcos(tx)w(x)dx. If wis a pdf with cf ϕwthen Iw(t) = Re{ϕw(t)}, which clearly coincides with ϕwwhen wis a symmetric pdf. 3 Asymptotics In order to study the limit behaviour of the test statistic T1nwe first need to introduce some assumptions on the models (1) and on the available data. Recall that we are assuming that {(Xjl, Yjl),1≤l≤nj}are iid observations from (Xj, Yj)and the sets {(X1l, Y1l),1≤l≤n1},. . .,{(Xjk, Yjk),1≤l≤nk}are independent. Assumption (A): (A.1) For 1≤j≤k: (i) Xjhas a compact support R. (ii) fj,mjand σjare two times continuously differentiable on R. (iii) infx∈Rfj(x)>0and infx∈Rσj(x)>0. (A.2) For 1≤j≤k: the samples sizes satisfy lim nj/n =pj, where 0< pj<1. (A.3) Kis a twice continuously differentiable symmetric pdf with compact support. (A.4) The weight function satisfies w(t)≥0, for all t∈R, and Rt4w(t)dt < ∞. (A.5) nh4 n→0and nh2 n/ln n→ ∞. These assumptions are mainly needed to guarantee the uniform consistency of the kernel estimators ˆ fj,ˆσj,ˆmjand ˆm0. Unlike the methods based on the ecdf, observe that we do not impose any restriction on the distribution of the errors, like the existence of a pdf. So the results in this paper could be used to compare two or more regression functions when the distributions of the errors are arbitrary: continuous, discrete or mixed. An example with errors having a mixed distribution is given in Section 4. 7 3.1 Asymptotic null distribution The following theorem gives an asymptotic approximation for √nj{ϕj(t)−ϕ0j(t)},1≤ j≤nj, that will allow us to derive the asymptotic null distribution of the test statistic T1n, given in the subsequent corollary. Let Σ = (σjv)1≤j,v≤kbe the matrix whose elements are σjj = 1 −2pjEfj(Xj) fmix(Xj)+pj k X r=1 prEσ2 r(Xr) σ2 j(Xr) f2 j(Xr) f2 mix(Xr), σjv =√pjpv k X r=1 prEσ2 r(Xr) σj(Xr)σv(Xr) fj(Xr)fv(Xr) f2 mix(Xr) −√pjpvEσv(Xv) σj(Xv) fj(Xv) fmix(Xv)+σj(Xj) σv(Xj) fv(Xj) fmix(Xj), j 6=v. (7) Theorem 2 Under assumptions (A.1)-(A.3) and (A.5), if H0is true, then √nj{ˆϕj(t)−ˆϕ0j(t)}=itϕj(t)Zj+tR1j(t) + t2R2j(t), where sup t|Rsj(t)|=op(1),s= 1,2, and Z:= (Z1, . . . , Zk)0∼Nk(0,Σ). Define the diagonal matrix A=diag(a1, . . . , ak), where aj=ktϕj(t)k2 w,1≤j≤k. The results in the theorems and the corollaries below will hold whenever trace(AΣ) >0. Before stating the results, we briefly discuss this condition. Observe that trace(AΣ) = k X j=1 ajσjj >0if and only if aj>0and σjj >0for some j , 1≤j≤k. The quantities σjj in (7) can be also expressed as σjj =pj k X l=1 plE"σ2 l(Xl) σ2 j(Xl)fj(Xl) fmix(Xl)−I(l=j) pl2#, where I(A)denotes the indicator function of the set A. From Assumptions (A.1)(iii) and (A.2) it follows that σjj >0for all j. Thus, to ensure trace(AΣ) >0we only need to ensure that aj>0for some j. An easy way to get aj>0is by taking w(t)>0, for tin a neighborhood of the origin. The following assumption will appear in the statement of some of the results below. Assumption (B): aj>0for some 1≤j≤k. 8 Corollary 3 Under Assumptions (A) and (B), if H0is true, then nT1n L −→ W1=Z0AZ, where Zis as in Theorem 2. In other words, the limiting distribution of nT1nunder H0is a finite linear combination of independent chi-square variables, Pk j=1 βjχ2 1,j, where χ2 1,1, . . . , χ2 1,k are independent chisquare random variates with one degree of freedom and β1, . . . , βkare the eigenvalues of AΣ. Unfortunately, the quantities βjin this linear combination are unknown. They depend on the distribution of the errors through the a1, . . . , ak, and on the distribution of the covariates through Σ = (σjv)1≤j,v≤k. They also depend on the unknown design densities, fj, and the conditional variance functions, σ2 j. So to use Theorem 2 in practice, one first need to find a consistent estimator, say ˆ βj, for every βj,1≤j≤k. This can be easily done via plug-in method using the estimators defined above instead of the unknown functions ϕj,fj,fmix and σ2 j. In order to perform the test we also need to approximate the distribution of Pk j=1 ˆ βjχ2 1,j which can be done via Monte Carlo method or some numerical method (see for example Kotz et al, 1967, Castaño-Martínez et al, 2005). With such a distribution, we can finally get the critical value and/or the p-value for the test based on T1n. The next result states the validity of this procedure. To estimate ajwe replace ϕj(t)by ˆϕj(t)in its expression obtaining, ˜aj=Zt2|ˆϕj(t)|2w(t)dt =1 n2 j nj X r,s=1 Zt2cos{t(ˆεjr−ˆεjs)}w(t)dt =nj−1 nj ˆaj+1 njZt2w(t)dt, where ˆaj=−1 nj 2!X 1≤r<s≤nj D2Iw(ˆεjr −ˆεjs),1≤j≤k, (8) D2Iw(t) = ∂2 ∂t2Iw(t) = −Rt2cos(tu)w(u)du and Iwis as defined in Remark 1. Since Rt2w(t)dt is a constant term, we estimate ajby ˆaj, which resembles a U-statistic. Let ˆ A=diag(ˆa1,...,ˆak)and ˆ Σ = (ˆσjv)1≤j,v≤k, with ˆσjj = 1 −2ˆpjˆµj+ ˆpjPk r=1 ˆprˆµjjr, ˆσjv =pˆpjˆpvPk r=1 ˆprˆµjvr −pˆpjˆpv(ˆµjv + ˆµvj),j6=v, ˆpj=nj n,ˆµj=1 nj nj X l=1 ˆ fj(Xjl) ˆ fmix(Xjl),ˆµjv =1 nv nv X l=1 ˆσv(Xvl) ˆσj(Xvl) ˆ fj(Xvl) ˆ fmix(Xvl), ˆµjvr =1 nr nr X l=1 ˆσ2 r(Xrl) ˆσj(Xrl)ˆσv(Xrl) ˆ fj(Xrl)ˆ fv(Xrl) ˆ f2 mix(Xrl), 1≤j, v, r ≤k. 9 4.2 The weight function In the present section we discuss the role of the weight function w. For our test to be consistent, it is only needed a positive wsatisfying Rt4w(t)dt < ∞. However, from the results in Section 3, it is clear that the choice of the weight function affects the power. In fact, from Corollary 10, the asymptotic power of the test Ψ1,α is given by P((Z+µ)0A(Z+µ)> t1,α).(11) This probability depends on wthrough A=diag(a1, . . . , ak), where aj=Rt2|ϕj(t)|2w(t)dt, 1≤j≤k. Clearly, the optimal weight (the one that maximizes the power) depends on the characteristic functions ϕ1, . . . , ϕkand on µ1, . . . , µk. Because all these quantities are unknown and involve the unobservable residuals, formula (11) is of little help in practice. An alternative approach is to use different weight functions w1, . . . , wkfor the different populations. This leads to the following test statistic k X j=1 nj nkˆϕj(t)−ˆϕ0j(t)k2 wj. By Theorem 2, this quantity has the same asymptotic properties as T1n. In fact, one only needs to replace aj=ktϕj(t)k2 wby ¯aj=ktϕj(t)k2 wjfor the results given in Section 3 to continue to hold. Following the guidelines in Epps and Pulley (1983) (see also Epps, 2005, Jiménez-Gamero et al., 2009, Hušková and Meintanis, 2010), a reasonable choice for wj(t)is |ϕj(t)|2. This choice attempts to give high weight where the statistic ˆϕj(t) is a relatively precise estimator of ϕj(t). The problem is that ϕj(t)is unknown and needs to be estimated. Taking wj(t)∝ |ˆϕj(t)|2is not possible because R|ˆϕj(t)|2dt =∞, since ˆϕj(t)is a periodic function. To overcome this difficulty we could consider a kernel smoothing estimator, but its application requires to assume rather strong assumptions on the distribution of the errors. A more practical method avoiding the above difficulties is to use a parametric density function as a weight. The density function should put most of the weight near the origin, because the ecf estimates more accurately the population cf around t= 0. This approach is connected with the comparison of kernel density estimators as follows. Let Sbe a pdf symmetric around the origin and let wbe such that w1/2(x) = (2π)−1/2Rexp(itx)S(t)dt. Lemma 2.1 in Fan (1998) (see also Anderson et al, 1994, Henze et al, 2005, Hušková and Meintanis, 2012 and Meintanis, 2013) shows that Z|ˆϕj(t)−ˆϕ0j(t)|2wj(t)dt =Znˆ fj(t)−ˆ f0j(t)o2 dt, 16 where ˆ fj(t) = 1 n nj X l=1 S(t−ˆεjl),ˆ f0j(t) = 1 n nj X l=1 S(t−ˆε0jl), that is to say, for adequate choices of the weight function wj, the quantity kˆϕj−ˆϕ0jk2 w coincides with the integral of the squared of the difference between kernel estimators of the pdf of the errors, both estimated with bandwidth η= 1. This bandwidth can be taken arbitrary by considering Sη(x) = η−1S(x/η)instead of S=S1. The above observation does not narrow down the spectrum of possibilities. So further criteria must taken into account. From a practical point of view, the ease of computation of T1n(also T2n) is closely related to the choice of w. This is specially appealing if instead of using the approximation given in Theorem 4, one wishes to employ a bootstrap approximation, which requires to evaluate the test statistic in a high number of artificial samples. In this sense, a good choice is w(t)∝exp{−t2/2σ2 w},(12) which is tantamount to estimate the pdf of the errors by using the normal kernel and bandwidth σw. Nevertheless, at this point we must say that the choice of σwcannot be guided by standard results on kernel density estimation, because in such a case the asymptotic results in this paper are no longer true. A theoretical study of the optimal choice of σwin terms of Bahadur slopes can be found in Tenreiro (2009) for the problem of testing goodness-of-fit for the normal distribution. The results in the cited paper show that the optimal choice of σwdepends on the alternative, which is unknown in practice. Because of this reason, we did not pursue this line. Neither too large, nor too small values of σware appropriate because, proceeding as in Henze at al (2005), we get lim σw→0σ−2 wT1n= k X j=1 nj n(¯ ˆεj−¯ ˆε0j)2,(13) with ¯ ˆεj=1 njPnj j=1 ˆεjl,¯ ˆε0j=1 njPnj j=1 ˆε0jl,1≤j≤k. The right-side of (13) estimates θ=Pk j=1 pj{E(εj)−E(ε0j)}2. In certain situations θmight be zero, even when H0is not true and therefore the resulting test based on small σwwould not be consistent against all fixed alternatives. A similar phenomenon occurs when σw→ ∞. In the context of testing goodness-of-fit for the distribution of the errors in nonparametric regression models, the simulations in Hušková and Meintanis (2010) reveal that taking σw= 1 gives good results. We have also investigated numerically the effect of changing the parameter σwin models (iv) (level approximation) and (vi) (power) when the critical values are approximated from the asymptotic distribution of T1nas explained 17 in Theorem 4. Table 2 (which is just a part of a larger simulation study) summarizes the obtained results. It can be seen that the results are quite homogeneous with slightly better level with σw= 1, which is the value that we will consider in the rest of the simulations. [ Table 2 to be placed here ] 4.3 Results of the test based on the asymptotic null distribution In this section we present the results of the test based on the approximation of the asymptotic null distribution of T1nas given in Theorem 4. Tables 3 (homoscedastic models) and 4 (heteroscedastic models) display the results for bandwidths of the form h=Cn−0.375, with C= 1,1.5,2, which provide reasonable values for the considered setups. The level –models (i)-(iv)– is slightly overestimated for small sample sizes, but the approximation improves as the sample sizes increase, reaching a good approximation for (n1, n2) = (100,100). The test also reaches good power, both in the homoscedastic case and in the heteroscedastic case. We have also run simulations with error distributions different from the normal distribution. In particular, we have considered two cases: (a) errors with mixed distribution of the form εj=(0with probability 0.5 N(0,√2) with probability 0.5, j= 1,2, which in practice could model a case where the observations come from two devices, one of them with no measurement error; and (b) errors with a recentered Exponential distribution, i.e., εj+ 1 ∼Exponential(1), for j= 1,2. Table 5 displays a brief summary of the obtained results for models (ii), (iv), (vi) and (viii) under heteroscedasticity and sample sizes n1=n2= 100. In the case of the mixed distribution, the approximation of the level is very good. In the case on the exponential distribution, we observe an overestimation of the level, probably cause by a bad approximation of the asymptotic null distribution due to the asymmetry of the error distribution (although not shown in the table, better approximations are achieved for larger sample sizes). Both cases show good power. The same kind of conclusions can be established for the rest of the models, which are not shown here. [ Table 3 to be placed here ] [ Table 4 to be placed here ] [ Table 5 to be placed here ] 18 4.4 Results based on a bootstrap approximation The aim of this subsection is to compare the power of the test proposed in this paper with those in Pardo-Fernández et al. (2007). To approximate the null distribution of their test statistics, Pardo-Fernández et al. (2007) employed a bootstrap procedure based on smoothed residuals (see also Neumeyer, 2009, for a theoretical justification). Of course, the same bootstrap procedure could be used to approximate the null distribution of nT1n. Nevertheless, from a computational point of view, the estimators in Theorem 4 and Corollary 6 are less time consuming. In order to establish a fair comparison, we have also estimated the null distribution of the test proposed in this paper by using the bootstrap algorithm defined in the above-mentioned paper. Besides, we have also incorporated here the test statistic T2n, for which the asymptotic null distribution is difficult to approximate. Table 6 shows the results of the tests based on T1nand T2nand the four tests proposed in Pardo-Fernández et al. (2007), which are denoted by T1 KS,T2 KS,T1 CM and T2 CM . For the sake of brevity of the presentation of the table, we restrict ourselves to the significance level α= 0.05 and bandwidth with C= 1 (similar results have been obtained for other significance levels and other specifications of the bandwidth). In terms of level approximation, we can see that it is good for all test statistics, except for T2 KS. Compared to the asymptotic approximation, the bootstrap approximation improves the behaviour of the test statistic T1nfor small sample sizes. For models (v)-(vii), the highest power is achieved by the test based on the ecf, T1n. For model (viii), the highest power is achieved by T2n, which is also based on the ecf. Note that in this model T1nreaches reasonable power and it is much better than its ecdf-based analog T1 CM . Summarizing, for the models considered, the test based on T1npresents, as a whole, the best behaviour. [ Table 6 to be placed here ] 4.5 The case of three populations We have also investigated the test based on the estimated asymptotic null distribution of T1nin the case of three populations (k= 3). Now the regression models are: (ix) m1(x) = m2(x) = m3(x)=1. (x) m1(x) = m2(x) = m3(x) = x. (xi) m1(x) = x,m2(x) = x+ 0.2,m3(x) = x+ 0.4. (xii) m1(x) = x,m2(x) = x,m3(x) = x+ 0.25. 19 (xiii) m1(x) = 0.5,m2(x) = x,m3(x) = 1 −x. (xiv) m1(x) = 0,m2(x) = sin(2πx),m3(x) = −sin(2πx). Models (ix)-(x) are under the null hypothesis, models (xi)-(xiv) are under the alternative. We only consider homoscedastic models with scale functions σ1(x) = √0.25,σ2(x) = √0.25 and σ3(x) = √0.50. The covariates X1,X2and X3are Beta(1.5,2),Beta(2,1.5) and Beta(2,2), respectively, and all regression errors are N(0,1). As in the previous cases, a bandwidth of the form h=Cn−0.375 is chosen, but now the C= 2,2.5,3are displayed. Other choices for Cwere also tried, but better results where obtained for these values. The results are shown in Table 7. As in Tables 3 and 4, the level is well approximated for large sample sizes and the behaviour in terms of power is correct. [ Table 7 to be placed here ] 5 Conclusions A test for the comparison of kregression functions has been proposed and studied under a totally nonparametric setting. The test statistic compares the ecf of the residuals in each population with the ecf of the residuals under the null hypothesis. For adequate choices of the weight function involved in the definition of the test statistic, the resulting test is consistent against any fixed alternative and is able to detect contiguous alternatives converging to the null at a rate n−1/2. To derive these properties we have assumed certain assumptions, which are weaker than those required by those based on the ecdf. Specifically, no requirement is imposed on the distributions of the errors. An estimation of the asymptotic null distribution has been proposed as an estimator of the null distribution of the test statistic. In the cases tried in our numerical experiments it is observed that this approximation works, in the sense of providing type I errors close to the nominal values, specially when the sample sizes are at least 100. For smaller sample sizes it is recommended to approximate the null distribution through a bootstrap mechanism. 20 6 Appendix See the supplementary material available from the journal website. Acknowledgements The authors thank the anonymous referees for their constructive comments and suggestions which helped to improve the presentation. J.C. Pardo-Fernández has been financially supported by grant MTM2011-23204 (FEDER support included) of the Spanish Ministry of Science and Innovation. M.D. Jiménez-Gamero has been financially supported by grant UJA2013/08/01 (University of Jaén and Caja Rural Provincial of Jaén). A. El Ghouch acknowledges financial support from IAP research network P6/03 of the Belgian Government (Belgian Science Policy), and from the contract ‘Projet d’Actions de Recherche Concertées’ (ARC) 11/16-039 of the ‘Communauté française de Belgique’, granted by the ‘Académie universitaire Louvain’. References Alba Fernández, V., Jiménez-Gamero, M.D., Muñoz-García, J. (2008). A test for the twosample problem based on empirical characteristic functions. Computational Statistics & Data Analysis, 52, 3730–3748. Anderson, N.H., Hall, P., Titterington, D.M. (1994). Two-sample test statistics for measuring discrepancies between two multivariate probability density functions using kernel-based density estimates. Journal of Multivariate Analysis, 50, 41–54. Castaño-Martínez, A., López-Blázquez, F. (2005). Distribution od a sum of weighted central chi-square variables. Communications in Statistics-Theory and Methods, 34, 515–524. de Uña-Álvarez, J. (2013). Comments on: An updated review of Goodness-of-Fit tests for regression models. Test, 22, 414–418. Delgado, M.A. (1993). Testing the equality of nonparametric regression curves. Statistics and Probability Letters, 17, 199–204. Epps, T.W. (2005). Tests for location-scale families based on the empirical characteristic function. Metrika, 62, 99–114. 21 Epps, T.W., Pulley, L.B. (1983). A test for normality based on the empirical characteristic function. Biometrika, 70, 723–726. Fan, Y. (1998). Goodness-of-fit tests based on kernel density estimators with fixed smoothing parameters. Econometric Theory, 14, 604–621. Fan, J., Gijbels. I. (1996). Local Polynomial Modelling and Its Applications. Chapman & Hall. Feller, W. (1971). An Introduction to Probability Theory and its Applications, Vol 2. Wiley. González-Manteiga, W. and Crujeiras, R. (2013). An updated review of Goodness-of-Fit tests for regression models. Test, 22, 361–411. Henze, N., Klar, B., Zhu, L.X. (2005) Checking the adequacy of the multivariate semiparametric location shit model. Journal of Multivariate Analysis, 93, 238–265. Heuchenne, C., Van Keilegom, I. (2010). Goodness-of-fit tests for the error distribution in nonparametric regression. Computational Statistics and Data Analysis, 54, 1942–1951. Hušková, M., Meintanis, S. (2007). Omnibus tests for the error distribution in the linear regression model. Statistics, 41, 363–376. Hušková, M., Meintanis, S.G. (2009). Goodness-of-fit tests for parametric regression models based on empirical characteristic functions. Kybernetika, 45, 960–971. Hušková, M., Meintanis, S.G. (2010). Tests for the error distribution in nonparametric possibly heteroscedastic regression models. Test, 19, 92–112. Hušková, M., Meintanis, S.G. (2012). Tests for symmetric error distribution in linear and nonparametric regression models. Communications in Statistics-Simulation and Computation, 41, 833–851. Jiménez-Gamero, M.D., Muñoz-García, J., Pino-Mejías, R. (2005). Testing goodness of fit for the distribution of errors in multivariate linear models. Journal of Multivariate Analysis, 95, 301–322. Jiménez-Gamero, M.D., Alba-Fernández, V., Muñoz-García, J., Chalco-Cano, Y. (2009). Goodness-of-fit tests based on empirical characteristic functions, Computational Statistics and Data Analysis, 53, 3957-3971. Kotz, S., Johnson, N.L., Boyd, D.W. (1967). Series representations of quadratic forms in normal variables. I. Central case. The Annals of Mathematical Statistics, 38 , 823–837. Kulasekera, K.B. (1995). Comparison of regression curves using quasi-residuals. Journal 22 of the American Statistical Association, 90, 1085–1093. Meintanis, S.G. (2013). Comments on: An updated review of Goodness-of-Fit tests for regression models. Test, 22, 432–436. Munk, A., Dette, H. (1998). Nonparametric comparison of several regression functions: Exact and asymptotic theory. Annals of Statistics, 26, 2339–2368. Neumeyer, N. (2009). Smooth residual bootstrap for empirical processes of non-parametric regression residuals. Scandinavian Journal of Statistics, 36, 204–228. Neumeyer, N., Dette, H. (2003). Nonparametric comparison of regression curves: an empirical process approach. The Annals of Statistics, 31, 880-920. Neumeyer, N., Pardo-Fernández, J.C. (2009). A simple test for comparing regression curves versus one-sided alternatives. Journal of Statistical Planning and Inference, 139, 4006-4016. Pardo-Fernández, J.C., Van Keilegom, I., González-Manteiga, W. (2007). Testing for the equality of kregression curves. Statistica Sinica, 17, 1115-1137. Sperlich, S. (2013). Comments on: An updated review of Goodness-of-Fit tests for regression models. Test, 22, 419–427. Srihera, R., Stute, W. (2010). Nonparametric comparison of regression functions. Journal of Multivariate Analysis, 101, 2039–2059. Tenreiro, C. (2009). On the choice of the smoothing parameter for the BHEP goodness- of-fit. Computational Statistics and Data Analysis, 53, 1038–1053. 23 Table 1: Empirical level of the test based on the asymptotic distribution of T1nfor homoscedastic and heteroscedastic models and smoothing parameters chosen by crossvalidation with n1=n2= 100. homoscedastic models heteroscedastic models model α:0.100 0.050 0.010 0.100 0.050 0.010 (i) 0.154 0.064 0.018 0.148 0.062 0.012 (ii) 0.136 0.052 0.012 0.124 0.052 0.010 (iii) 0.152 0.088 0.016 0.152 0.086 0.012 (iv) 0.106 0.048 0.010 0.106 0.046 0.010 Table 2: Observed rejection proportions of the test based on the asymptotic distribution of T1nfor homoscedastic models (iv) and (vi) and for different choices of the parameter σwwith n1=n2= 100. α: 0.100 0.100 0.100 0.050 0.050 0.050 0.010 0.010 0.010 model σwC: 1.0 1.5 2.0 1.0 1.5 2.0 1.0 1.5 2.0 (iv) 0.50 0.100 0.082 0.078 0.046 0.041 0.039 0.012 0.009 0.007 0.75 0.111 0.090 0.083 0.049 0.043 0.039 0.012 0.009 0.008 1.00 0.122 0.097 0.084 0.052 0.046 0.041 0.012 0.009 0.008 1.25 0.133 0.100 0.086 0.056 0.047 0.042 0.012 0.010 0.008 1.50 0.138 0.109 0.093 0.065 0.051 0.044 0.013 0.010 0.008 (vi) 0.50 1.000 1.000 1.000 1.000 1.000 1.000 0.996 0.997 0.996 0.75 1.000 1.000 1.000 1.000 1.000 1.000 0.996 0.997 0.997 1.00 1.000 1.000 1.000 1.000 1.000 1.000 0.996 0.997 0.997 1.25 1.000 1.000 1.000 1.000 1.000 1.000 0.997 0.997 0.997 1.50 1.000 1.000 1.000 1.000 1.000 1.000 0.999 0.997 0.997 Table 3: See the supplementary material available from the journal website. Table 4: See the supplementary material available from the journal website. 24 Table 5: Observed rejection proportions of the test based on the asymptotic distribution of T1nfor different error distribution with n1=n2= 100. The models are for heteroscedastic. α: 0.100 0.100 0.100 0.050 0.050 0.050 0.010 0.010 0.010 error distribution model C: 1.0 1.5 2.0 1.0 1.5 2.0 1.0 1.5 2.0 mixed (ii) 0.104 0.094 0.093 0.051 0.050 0.047 0.010 0.007 0.008 (iv) 0.105 0.090 0.083 0.049 0.046 0.044 0.011 0.008 0.007 (vi) 0.996 1.000 0.999 0.992 0.998 0.998 0.985 0.993 0.990 (viii) 0.989 0.997 0.994 0.987 0.992 0.985 0.950 0.917 0.836 exponential (ii) 0.149 0.125 0.115 0.076 0.063 0.055 0.015 0.014 0.010 (iv) 0.146 0.115 0.105 0.077 0.061 0.046 0.015 0.012 0.008 (vi) 0.997 0.999 0.999 0.997 0.999 0.998 0.992 0.990 0.988 (viii) 0.998 1.000 0.995 0.995 0.990 0.981 0.948 0.904 0.812 Table 6: See the supplementary material available from the journal website. Table 7: See the supplementary material available from the journal website. 25 View publication statsView publication stats