Testing for superiority between two variance functions

Statistics and Probability Letters 156 (2020) 108612

Contents lists available at ScienceDirect

Statistics and Probability Letters journal homepage: www.elsevier.com/locate/stapro

Testing for superiority between two variance functions Graciela Boente a , Juan Carlos Pardo–Fernández b ,

∗

a

Departamento de Matemáticas and Instituto de Cálculo, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires and CONICET, Ciudad Universitaria, Pabellón 1, 1428, Buenos Aires, Argentina b Departamento de Estatística e Investigación Operativa, Universidade de Vigo, Campus Universitario As Lagoas-Marcosende, Vigo, 36310, Spain

article

info

a b s t r a c t This paper focuses on the problem of testing the null hypothesis that the variance functions of two populations are equal versus one-sided alternatives under a general nonparametric heteroscedastic regression model. The asymptotic distribution of the test statistic is studied. © 2019 Elsevier B.V. All rights reserved.

Article history: Received 27 February 2019 Received in revised form 19 August 2019 Accepted 1 September 2019 Available online 16 September 2019 Keywords: Conditional variance Hypothesis testing Nonparametric regression models Smoothing techniques

1. Introduction Let us assume that the random vectors (Xj , Yj )t in R2 , j = 1, 2, follow the heteroscedastic nonparametric regression models given by Yj = mj (Xj ) + σj (Xj )ϵj ,

(1)

where mj : R → R and σj : R → R are nonparametric smooth functions and the error ϵj is such that E[ϵj |Xj ] = 0 and E[ϵj2 |Xj ] = 1, to identify the regression and variance functions, respectively, so that mj (x) = E[Yj |Xj = x] and σj2 (x) = Var[Yj |Xj = x], for j = 1, 2. The covariates Xj , j = 1, 2, are assumed to be continuous with density fj and with common support R. The nonparametric model (1) offers an alternative to the standard linear model allowing for more complex relations between the response variable and the covariate. It also includes in particular the situation in which the error ϵj is independent of the covariate Xj , often called location–scale model. When considering two independent populations following model (1), the interest has mainly focused on comparing the regression curves mj . The problem of one-sided regression alternatives has been considered, among others, by Koul and Schick (2003), Neumeyer and Dette (2005), Neumeyer and Pardo-Fernández (2009) and Boente and Pardo-Fernández (2016). However, in many situations, it is of interest to compare the conditional variance functions σ12 and σ22 to check if the same dispersion structure appears in both populations. Furthermore, from some previous knowledge the practitioner may know that σ12 (x) ≤ σ22 (x) and then the null hypothesis of equality of the variance curves versus a one-sided alternative is in this case the target. One of such examples may arise in quality control, where it is important to verify the stability and the homogeneity of the production process under different experimental and practical conditions. Under a controlled framework, the process may have larger variability under one setting than under the other one. Without controlling for the effect of covariates, the classical F -test based on the ratio of the sample variances is a well-known procedure. +

∗ Corresponding author. E-mail addresses: [email protected] (G. Boente), [email protected] (J.C. Pardo–Fernández). https://doi.org/10.1016/j.spl.2019.108612 0167-7152/© 2019 Elsevier B.V. All rights reserved.

2

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

In this paper, we focus on comparing conditional variances, that is, we assume that model (1) holds and our goal is to test the null hypothesis of equal variance curves versus a one-sided alternative. The comparison will be performed in the common support of the covariates X1 and X2 , which is denoted by R. Let A ⊂ R be such that P(Xj ∈ A) > 0, for j = 1, 2. The null hypothesis is then H0 : σ1 (x) = σ2 (x) for all x ∈ R, and the alternative hypothesis is of the following one-sided type H1 : σ1 (x) ≤ σ2 (x) for all x ∈ R and σ1 (x) < σ2 (x) for x ∈ A .

(2)

The problem of testing equality of several variance curves has been considered in Pardo-Fernández et al. (2015). Their approach is based on empirical distribution functions and empirical characteristic functions of residuals. In the particular case of two populations, Pardo-Fernández et al. (2015) focus on the alternative hypothesis H1 : σ1 (x) ̸ = σ2 (x), while in our setting we assume an ordered alternative. This will allow us to give a simpler test statistic, as we will see in the next section. As in Pardo-Fernández et al. (2015), the regression and variance functions as well as the distribution of the errors and covariates are unknown. The proposed methodology is therefore fully nonparametric. The aim of this paper is to propose a class of tests for H0 versus H1 in (2), detecting root−n alternatives, based on smooth estimators of the variance functions combined with the ideas given in Neumeyer and Pardo-Fernández (2009) to compare regression functions. The paper is organized as follows. In Section 2, we introduce the tests statistics and their asymptotic behaviour under the null and contiguous alternatives is derived in Section 3. The results of a numerical experiment designed to evaluate the finite-sample performance of the proposed procedure are included in Section 4. Section 5 presents the analysis of a real dataset. All proofs are deferred to the Supplementary material available online. 2. The test statistic In this paper, we consider independent and identically distributed observations (Xij , Yij )t , 1 ≤ i ≤ nj , with the same distribution as (Xj , Yj ), j = 1, 2. The regression functions mj (x) in (1) equal E[Yj |Xj = x] and can be estimated using a Nadaraya–Watson estimator or a local-linear estimator (see, for example, Fan and Gijbels, 1996). For the sake of simplicity, we will use a Nadaraya–Watson estimator to describe the proposed methodology. Let K be a kernel function (usually a symmetric density) and h = hn a sequence of strictly positive real numbers. Denote as Kh (u) = h−1 K (u/h) and let Wij (x) =

{∑n

j

(

Kh x − Xij

i=1

{ nj ∑

ˆj (x) = m

)}−1

(

Kh x − Xij

(

)

}−1 nj ∑ )

(

Kh x − Xij . Then, the Nadaraya–Watson regression estimators of mj are defined as

i=1

)

Kh x − Xij Yij =

i=1

nj ∑

Wij (x)Yij .

(3)

i=1

On the other hand, under model (1), it holds σj2 (x) = E[(Yj − mj (x))2 |Xj = x] and then, as in Pardo-Fernández et al. (2015), the following estimator of the variance function can be considered:

ˆ σj2 (x) =

nj ∑

(

ˆj (x) Wij (x) Yij − m

)2

.

(4)

i=1

As mentioned in the Introduction, our aim is to develop a class of consistent tests for H0 versus H1 in (2) allowing to detect root−n local alternatives. We will use similar ideas to the ones in Neumeyer and Pardo-Fernández (2009) for regression functions. Let σ be any function satisfying σ1 (x) ≤ σ (x) ≤ σ2 (x) for all x ∈ R, and define the random variables εj0 = (Yj − mj (Xj ))/σ (Xj ) = ϵj σj (Xj )/σ (Xj ), for j = 1, 2. Moreover, let wj : R → R be non-negative weight ◦

functions with compact support Sj ⊂ R such that A ∩ S1 ∩ S2 ̸ = ∅. Then, taking conditional expectation and using that σ1 (x) ≤ σ (x) ≤ σ2 (x) and that E[ϵj |Xj ] = 0 and E[ϵj2 |Xj ] = 1, we have that

)] σ12 (X1 ) a1 := E ε − 1 w1 (X1 ) = E w1 (X1 ) −1 ≤0, σ 2 (X1 )

(5)

)] [( 2 ) ] σ22 (X2 ) − 1 = E ε20 − 1 w2 (X2 ) =: a2 . 2 σ (X2 )

(6)

[(

2 10

)

[

]

(

while

[

0 ≤ E w2 (X2 )

(

On the one hand, under the null hypothesis H0 , the inequalities in (5) and (6) are equalities. On the other hand, under the alternative hypothesis, at least one of the inequalities in (5) and (6) must be strict if P(Xj ∈ A ∩ S1 ∩ S2 ) > 0. Effectively, assume that a1 = a2 . Using that a1 ≤[ 0 ≤ a2 , we obtain]that a1 = a2 = 0. Hence, taking into account that a1 = − E[V1 ] [ ] and a2 = E[V2 ], where V1 = w1 (X1 ) 1 − σ12 (X1 )/σ 2 (X1 ) and V2 = w2 (X2 ) σ22 (X2 )/σ 2 (X2 ) − 1 are nonnegative random variables, we obtain that

(

[

P wj (Xj ) 1 −

σj2 (Xj ) σ 2 (Xj )

]

) = 0 = 1 for j = 1, 2 .

(7)

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

3

The fact that P(Xj ∈ A ∩ S1 ∩ S2 ) > 0 together with (7), entails that σ1 (x) = σ2 (x) for x ∈ A ∩ S1 ∩ S2 , which contradicts (2). Therefore, we have that D21 ≥ 0, where

)] [ ( ) ( 2 )] ∑ 2 σj2 (Xj ) σ1 (X1 ) σ22 (X2 ) j − 1 − w1 (X1 ) −1 = −1 , = E w2 (X2 ) (−1) E wj (Xj ) σ 2 (X2 ) σ 2 (X1 ) σ 2 (Xj ) [

D21

(

(8)

j=1

and equality holds if and only if H0 is true. To construct the test statistic, we use (8) and define the common variance function under the null hypothesis as σ02 (x) = π1 (x)σ12 (x) + π2 (x)σ22 (x), for x ∈ R, where 0 ≤ π1 (x) ≤ 1 is a given function and π2 (x) = 1 − π1 (x). An estimator of σ02 (x) can easily be defined from estimators of σj2 (x) as ˆ σ02 (x) = 2 2 π1 (x)ˆ σ1 (x) + π2 (x)ˆ σ2 (x). Then, one may consider the empirical version of the left hand side in (8), that is, T =

2 ∑

n

(−1)

j=1

j 1 ∑

j

nj

i=1

(

ˆ σj2 (Xij ) ˆ σ02 (Xij )

) − 1 wj (Xij ) ,

which, under uniform consistency of ˆ σj , provides a consistent estimator of D21 . The test statistic to be considered is then defined as T1 =

( n n )1/2 1 2 n

T,

(9)

where n = n1 + n2 , and the null hypothesis is rejected for large positive values of T1 . [( 2 ) ] test statistic can be defined by noting that from (5) and (6) we have E ε10 − 1 w1 (X1 ) ≤ 0 ≤ ) ] [(Another 2 ˆj (Xij ))/ˆ − 1 w2 (X2 ) , where equality holds if and only if H0 is true. Then, if we define ˆ εij,0 = (Yij − m σ0 (Xij ), an E ε20 alternative test statistic would be T2 =

( n n )1/2 1 2 n

T⋆ ,

(10)

with T⋆ =

n2 1 ∑(

n2

ˆ εi22 ,0 − 1 w2 (Xi2 ) − )

i=1

n1 1 ∑(

n1

ˆ εi12 ,0 − 1 w1 (Xi1 ) = )

i=1

n

2 ∑

(−1)j

j=1

j 1 ∑

nj

i=1

{[

ˆj (Xij ) Yij − m ˆ σ02 (Xij )

]2

} − 1 wj (Xij ) .

The null hypothesis is also rejected for large positive values of T2 . 3. Asymptotic behaviour The following regularity assumptions are needed in order to prove our asymptotic results: A1 For j = 1, 2, the functions mj and σj are twice continuously differentiable in a neighbourhood of the common support, R. Furthermore, i(σj ) = infx∈Sj σj (x) > 0. A2 E[ϵj4 ] < ∞. Furthermore, the function µ4,ℓ (x) = E (ϵℓ2 − 1)2 |Xℓ = x is bounded in a neighbourhood of R.

[

]

◦

A3 For j = 1, 2, wj : R → R are bounded non-negative continuous weight functions with compact support Sj ⊂ R such that A ∩ Sj ̸ = ∅. The function π1 (x) is continuous in a neighbourhood of Sj . A4 For j = 1, 2, the density fj of Xj is twice continuously differentiable in a neighbourhood of S1 ∪ S2 . Moreover, i(fj ) = infx∈S1 ∪S2 fj (x) > 0. ∫ A5 K : R → R is an even and Lipschitz continuous function with bounded support, say [−1, 1], and K (u)du = 1. A6 nj /n → κ√ j as n = n1 + n2 → ∞, where 0 < κj < 1. A7 hn → 0, nh2n /log n → ∞, nh4 → 0 as n → ∞. It is worth mentioning that condition nh4 → 0 in A7 implies that under-smoothed kernel estimators are considered. Assume that (1) holds. Let ˆ σ02 (x)

Theorem 1.

∑2

j=1

wj (x)π3−j (x)fj (x). Then, under A1 to A7,

= π1 (x)ˆ σ12 (x) + π2 (x)ˆ σ22 (x) , with ˆ σj2 defined in (4) and f (x) =

D

(a) If H0 : σ1 = σ2 holds, we have that T1 −→ N(0, υ 2 ) where

υ = 2

2 ∑ j=1

[ κ3−j E

f 2 (Xj ) fj2 (Xj )

] µ4,j (Xj ) .

(11)

4

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

(b) Let ∆ : R → R be a differentiable non-negative function. Then, under H1,n : σ22 (x) = σ22,n (x) = σ12 (x) + n−1/2 ∆(x), we D

have that T1 −→ N(c , υ 2 ) where c=κ

1/2 1

1/2 2 E

κ

⎤ ⎡ 2 ∑ π 3−j (Xj ) ⎣ ∆(Xj )⎦ . σ12 (Xj ) j=1

(12)

Remark 1. In the special case of location–scale models, that is, when the error ϵj is independent of the covariate Xj , the asymptotic variance of the test statistic T1 equals

υ = 2

2 ∑

[ κ3−j E

j=1

f 2 (Xj )

]

fj2 (Xj )

µ4,j ,

(13)

where µ4,ℓ = E[(ϵℓ2 − 1)2 ] is the unconditional centred fourth moment of the regression error. Furthermore, if additionally X1 and X2 have a common density f1 = f2 = f , then

⎡⎛ ⎞2 ⎤ 2 2 ∑ ⎢ ⎥ ∑ wj (X1 )π3−j (X1 )⎠ ⎦ κ3−j µ4,j . υ 2 = E ⎣⎝ j=1

j=1

In particular, when w1 (x) = w2 (x) = 1 for all x, υ 2 does not depend on the covariate distribution. Remark 2. It is worth noting that the alternative hypothesis affects only the asymptotic mean which depends now on the unknown variance function σ12 (x). Remark 3. The results in Theorem 1 still hold if different bandwidths are used for estimating the regression and the variance function and/or if different bandwidth are used for the two populations, as long as the involved smoothing parameters satisfy A7. The main reason for the first assertion is that, if the bandwidth λ is used for the regression function and a bandwidth h is considered for the variance function and we explicitly denote as

ˆj,λ (x) = m

{ nj ∑

(

Kλ x − Xij

}−1 nj ∑ )

i=1

(

)

Kλ x − Xij Yij =

i=1

∑nj

nj ∑

Wij (x, λ)Yij

i=1

)2

ˆj,h (x) − m ˆj,λ (x) ˆj,λ (x) , we have that ˆ σj2,h,h (x) + m σj2,h,λ (x) = ˆ and ˆ σj2,h,λ (x) = i=1 Wij (x, h) Yij − m −1/2 has uniform rate n when h and λ satisfy A7. (

(

)2

and the second term

σ22 (x) and f (x) = σ12 (x) + π2 (x)ˆ Theorem 2. Assume that (1) holds. Let ˆ σj2 be defined in (4) and denote as ˆ σ02 (x) = π1 (x)ˆ

∑2

j=1

wj (x)π3−j (x) fj (x). Then, under A1 to A7,

D

(a) if H0 : σ1 = σ2 holds, we have that T2 −→ N(0, υ 2 ) where υ 2 is given in (11). (b) Let ∆ : R → R be a differentiable non-negative function. Then, under H1,n : σ22 (x) = σ22,n (x) = σ12 (x) + n−1/2 ∆(x), we D

have that T2 −→ N(c , υ 2 ) where c is given in (12). Remark 4. Note that, in sight of the results in Theorems 1 and 2 the test statistics T1 and T2 are asymptotically equivalent. 3.1. Practical implementation of the test Theorems 1-(a) and 2-(a) give the asymptotic null distribution of the test statistics T1 and T2 , respectively. We will focus on the practical implementation of the statistic T1 . An analogous construction is valid for T2 . First of all, in order the compute the test statistic, estimators of the regression curves and the conditional variances are needed, as given in (3) and (4). In our simulation study (Section 4) and illustration to real data (Section 5), we will use least-squares cross-validation bandwidths to estimate each one of the four curves separately. More precisely, to obtain ˆj , we employ the usual cross-validation bandwidth selector as described, for the estimate of the regression function, m example, in Fan and Gijbels (1996), page 149. For the conditional variance estimator, ˆ σj2 , we employ the same procedure 2 ˆj (Xij )) , i = 1 . . . , nj . by replacing the responses by the squared residuals (Yij − m The asymptotic variance υ 2 in (11) can be consistently estimated in general by

ˆ υ2 =

2 ∑ j=1

n

ˆ κ3−j

j 1 ∑ˆ f 2 (Xij )

nj

i=1

ˆ fj2 (Xij )

ˆ µ4,j (Xij ) ,

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

5

∑2

ˆ µ4,j (x) is a local smoother based where ˆ κj = n(j /n, , ˆ fj is the kernel density estimator of fj , ˆ f (x) = j=1 wj (x)π3−j (x) fj (x) and ˆ ) ˆj (Xij ))2 /ˆ on the pairs Xij , {ξij − 1}2 , 1 ≤ i ≤ nj , with ξij = (Yij − m σj2 (Xij ). However, we should notice that, in practice, the estimation of the conditional fourth moment µ4,j (x) might be cumbersome, specially for small or moderate sample sizes. This can be avoided by assuming independence between the error ϵj and the covariate Xj . In that case, the quantities µ4,j in (13) are unconditional and can be consistently estimated by using the sample variance of {ξij , i = 1, . . . , nj }, that is ˆ µ4,j = (nj − 1)−1 ˆ υ2 =

2 ∑

∑nj

1 ξ − ξ¯ )2 , with ξ¯ = n− j

i=1 ( ij

∑nj

i=1

ξij . The resulting estimator of υ 2 is

n

ˆ κ3−j ˆ µ4,j

j=1

j 1 ∑ˆ f 2 (Xij )

nj

i=1

ˆ fj2 (Xij )

.

(14)

In our simulations in Section 4 and illustration with real data in Section 5 we will restrict our attention to location-scale models and therefore we will implement the test with the estimator of υ 2 given in (14). √ υ 2 exceeds the (1 − α )-quantile of Finally, according to Theorem 1-(a), the test that rejects the null hypothesis if T1 / ˆ the standard normal distribution will have asymptotic level α . 4. Monte carlo study This section includes the results of a simulation study designed to analyse the finite-sample performance of the testing procedures proposed in Section 2. We consider a location–scale model, that is, the situation where the errors are independent of the covariates. Taking into account that the tests based on the statistics T1 and T2 given in (9) and (10), respectively, showed a similar performance, we only report here the numerical results for the former one. More specifically, we use the normal approximation for T1 given in Theorem 1 to compute the critical values. The asymptotic variance is estimated by (14). Table 1 reports the observed frequency of rejections with nominal significance levels α = 0.05 and α = 0.10, over 1000 replications. We choose m1 (x) = x and m2 (x) = sin(2π x). The covariates Xj are generated with uniform distribution on (0, 1) and we set w1 = w2 = I(0,1) and π1 (x) = 0.5. The following four scenarios√were considered to simulate the regression errors, j = 1, 2: NOR: ϵj ∼ N(0, 1); UNIF: ϵj ∼ U (−c , c), where c = 0(.5 12; EXP: ) ϵ√j + 1 ∼ Exponential(1); LGA: ϵj has a centred and standardized log–Gamma distribution, that is, ϵj = Wj − E[Wj ] / Var(Wj ), with Wj = log(Zj ) and Zj ∼ Γ (3, 3), where Γ (α, µ) denotes the mean parametrization of the Gamma distribution, i.e., if Zα,µ ∼ Γ (α, µ), then E[Zα,µ ] = µ and Var(Zα,µ ) = µ2 /α . It is worth noting that if ′ ′ Zα,µ ∼ Γ (α, µ), then E[log(Zα,µ ( )] = − log( ) α/µ √ ) + Ψ (α ) and Var[log(Zα,µ )] = Ψ (α ) with Ψ (α ) = Γ (α )/Γ (α ). Hence, under the scenario LGA, ϵj = Wj − Ψ (3) / Ψ ′ (3). The quantities Ψ (3) and Ψ ′ (3) were computed using the function psigamma in R. Different choices for the scale functions are considered to take into account the homoscedastic and heteroscedastic situations. For that reason, we select three scenarios denoted A, B and C followed by a 0 under the null hypothesis and by a 1 or a 2 under the considered alternatives. More precisely, for the null hypothesis the three situations correspond to A0: σ1 (x) = σ2 (x) = 0.25; B0: σ1 (x) = σ2 (x) = 0.25(1 + x); and C0: σ1 (x) = σ2 (x) = 0.5 exp(0.5 x). The associated alternatives are A1: σ1 (x) = 0.25 and σ2 (x) = σ1 (x) + 0.05; A2: σ1 (x) = 0.25 and σ2 (x) = σ1 (x) + 0.5x, B1: σ1 (x) = 0.25(1 + x) and σ2 (x) = σ1 (x) + 0.25; B2: σ1 (x) = 0.25(1 + x) and σ2 (x) = σ1 (x) + 0.25x; C1: σ1 (x) = 0.5 exp(0.5 x) and σ2 (x) = σ1 (x) exp(0.5 x); and C2: σ1 (x) = 0.5 and σ2 (x) = σ1 (x) exp(0.5 x). The testing procedure relies on nonparametric estimation of the regression and variance functions. The smoothing procedure for the regression function uses a local-linear estimator with the Epanechnikov kernel, while the variance estimators are based on the Nadaraya–Watson method, as this guarantees the positiveness of the obtained estimates. In all cases, the bandwidths were selected using least-squares cross-validation as described in Section 3.1. Table 1 reports the results for balanced designs with sample sizes n1 = n2 = 100, 200 and 400. The empirical level is close to the nominal one in all cases, except when n1 = n2 = 100 with exponential errors. This lack of accuracy may be explained by the fact that for strongly skewed distributions the required fourth moment estimates suffer from a large bias. Regarding the power performance, as expected higher powers are obtained for larger sample sizes. Besides, the test ability to detect alternatives is better for symmetric errors distributions (normal and uniform). To illustrate the behaviour when n1 ̸ = n2 , Table 2 reports the obtained frequencies of rejection under the null hypothesis for normal errors and unbalanced designs. The reported results show that the convergence of the empirical level is much slower than when n1 = n2 . Furthermore, the behaviour of the test is quite different for n1 < n2 and n2 < n1 . In the first situation, the test is liberal while in the second one it is conservative. To understand this behaviour, we have analysed the density of the standardized test statistics obtained by the 1000 Monte Carlo replications. Fig. 1 provides the kernel density estimates for model A0 for several choices of unbalanced sample sizes. In both plots, the solid black line is the target standard normal density and the vertical dotted line identifies its 95%-quantile. In Fig. 1(a), the dashed blue and the dotted–dashed red lines correspond to sample sizes (n1 , n2 ) = (400, 200) and (n1 , n2 ) = (200, 400), respectively. The vertical lines indicate their corresponding 95% empirical quantiles. Notice that, when (n1 , n2 ) = (400, 200), the density right tail lies below the desired asymptotic normal density, which implies that the resulting test is conservative. This can also be concluded when comparing the quantiles. In contrast, for (n1 , n2 ) = (200, 400), the density is mainly over the normal density. To illustrate that the convergence is still slow for large sample sizes, when n1 > n2 , Fig. 1(b) gives the

6

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

Table 1 Proportion of rejections in 1000 simulated datasets for different sample sizes n1 = n2 = 100, 200, 400.

α

Error:

NOR

model

100

200

400

UNIF 100

200

400

100

EXP 200

400

100

LGA 200

400

0.050

A0 B0 C0 A1 A2 B1 B2 C1 C2

0.039 0.050 0.048 0.438 0.998 0.992 0.705 0.680 0.666

0.039 0.047 0.046 0.709 1.000 1.000 0.941 0.916 0.904

0.041 0.052 0.052 0.950 1.000 1.000 0.999 0.997 0.996

0.042 0.050 0.055 0.682 1.000 1.000 0.918 0.895 0.893

0.040 0.044 0.048 0.949 1.000 1.000 0.996 0.996 0.995

0.039 0.048 0.049 1.000 1.000 1.000 1.000 1.000 1.000

0.070 0.081 0.079 0.272 0.907 0.838 0.432 0.405 0.390

0.051 0.058 0.059 0.383 0.992 0.969 0.596 0.559 0.544

0.038 0.041 0.040 0.538 1.000 1.000 0.813 0.762 0.758

0.060 0.062 0.068 0.415 0.993 0.983 0.676 0.634 0.613

0.038 0.040 0.045 0.645 1.000 1.000 0.878 0.865 0.857

0.038 0.043 0.046 0.891 1.000 1.000 0.991 0.985 0.982

0.100

A0 B0 C0 A1 A2 B1 B2 C1 C2

0.086 0.103 0.100 0.576 1.000 0.996 0.805 0.773 0.759

0.095 0.102 0.111 0.801 1.000 1.000 0.978 0.961 0.952

0.073 0.087 0.094 0.981 1.000 1.000 0.999 0.999 0.999

0.083 0.098 0.097 0.778 1.000 1.000 0.953 0.940 0.934

0.071 0.082 0.088 0.971 1.000 1.000 0.998 0.997 0.997

0.078 0.093 0.099 1.000 1.000 1.000 1.000 1.000 1.000

0.120 0.128 0.132 0.386 0.945 0.895 0.551 0.525 0.509

0.100 0.111 0.117 0.519 0.996 0.986 0.693 0.662 0.647

0.095 0.100 0.102 0.661 1.000 1.000 0.897 0.864 0.858

0.095 0.104 0.116 0.539 0.995 0.990 0.755 0.742 0.716

0.081 0.106 0.103 0.766 1.000 1.000 0.937 0.916 0.917

0.076 0.083 0.084 0.932 1.000 1.000 0.995 0.992 0.991

Table 2 Proportion of rejections in 1000 simulated datasets when n1 ̸ = n2 for normal errors.

α = 0.05 (n1 , n2 ) (200, 100) (100, 200) (400, 200) (200, 400) (800, 400) (400, 800) (1600, 800) (800, 1600)

Model:

α = 0.10

A0

B0

C0

A0

B0

C0

0.020 0.087 0.019 0.082 0.023 0.054 0.025 0.058

0.018 0.102 0.023 0.085 0.028 0.054 0.030 0.060

0.020 0.106 0.026 0.087 0.030 0.057 0.032 0.058

0.046 0.150 0.048 0.139 0.048 0.109 0.049 0.109

0.053 0.169 0.062 0.148 0.061 0.115 0.058 0.111

0.054 0.172 0.063 0.155 0.057 0.114 0.060 0.112

plots corresponding to sample sizes (n1 , n2 ) = (1600, 800) and (n1 , n2 ) = (800, 1600). In the first case, the density of the test statistic is still shifted to the left, whereas for n1 < n2 the normal approximation seems to be accurate (see also Table 2). 5. Data analysis To illustrate the performance of the proposed test, we consider the wine quality dataset, which consists of different measures regarding red and white variants of the Portuguese Vinho Verde wine. The data were collected from May 2004 to February 2007 using only protected designation of origin samples that were tested at the official certification entity and are available at https://archive.ics.uci.edu/ml/datasets/Wine+Quality. It contains 1599 red wines and 4898 white ones. The measured variables include: fixed acidity, volatile acidity, citric acid, residual sugar, chlorides, free sulphur dioxide, total sulphur dioxide, density, pH, sulphates and alcohol, which are based on physicochemical tests, and the quality of the wine, that is a score between 0 and 10. Details and a classification analysis based on support vector machine can be found in Cortez et al. (2009). In this paper, we focus on the relation between the pH and the fixed acidity (g(tartaric acid)/dm3 ), more precisely on the variance differences between red and white wines. It is worth noticing that, for the Vinho Verde wines, the pH varies between 2.7 and 4 for the red wines and between 2.7 and 3.8 for the white ones (see Table 1 in Cortez et al., 2009). Fig. 2(a) shows the scatter plot of fixed acidity versus pH, the red filled triangles correspond to red wines, while the grey circles to white wines. As an illustration, we have selected samples of size n1 = n2 = 200 from the red and white wines with pH between 3 and 3.6. The sample of white wines was labelled as population 1, while that of red wines as population 2, since it is expected that for white wines the fixed acidity will have smaller conditional variability than the red ones. Fig. 2(b) shows the plot of fixed acidity versus pH, for the selected samples together with the estimated regression curves. The solid grey ˆ1 (x), while the dashed red one to m ˆ2 (x). As it has been widely discussed in the literature, see, for line corresponds to m instance, Cortez et al. (2009), the fixed acidity increases with the wine acidity, that is, as the pH decreases. Furthermore, red wines tend to have higher values of the conditional mean fixed acidity than white wines. In this paper, we go beyond and we want to compare their variability, that is, the goal is to test equality between variances against the ordered

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

7

√

Fig. 1. Kernel density estimates obtained from the 1000 values of T1 / ˆ υ 2 , under model A0. The solid black line is the target standard normal density and the vertical dotted line identifies its 95%-quantile. (a) The dashed blue and the dotted-dashed red lines correspond to sample sizes (n1 , n2 ) = (400, 200) and (n1 , n2 ) = (200, 400), respectively. (b) The dashed blue and the dotted-dashed red lines correspond to sample sizes (n1 , n2 ) = (1600, 800) and (n1 , n2 ) = (800, 1600), respectively. The vertical lines indicate the corresponding 95% empirical quantiles.

Fig. 2. Scatter plot of fixed acidity versus pH, the red filled triangles correspond to red wines, while the grey circles to white wines. (a) Complete dataset. (b) Selected samples of size n1 = n2 = 200 together with the estimates of the regression functions; the solid grey and the dashed red lines correspond to the white and red wines, respectively. (c) Estimates of the conditional variances σj2 (x); the solid grey and the dashed red lines correspond to the white and red wines, respectively.

alternative. More precisely, if we denote as Xj the pH of the white (j = 1) and red wines (j = 2), our goal is to test H0 : σ1 (x) = σ2 (x) for all x ∈ [3, 3.6] against the alternative hypothesis H1 : σ1 (x) ≤ σ2 (x) for all x ∈ [3, 3.6] and σ1 (x) < σ2 (x) for x ∈ A. The estimates of σj2 (x) obtained as described in Section 2 are plotted in Fig. 2(c), which shows that the fixed acidity variance of white wines seems to be homoscedastic. On the other hand, the variability for red wines decreases as the pH increases. We performed the test to compare their variances and the resulting p-value is smaller than 0.0001. 6. Conclusion In this piece of research, we have developed a new procedure to test equality of two conditional variance functions against an ordered alternative. The asymptotic behaviour of the test statistics is derived. In practice, the testing procedure is simple to implement, as critical values and/or p-values can be derived from the asymptotic null distribution. Simulations show a correct behaviour in terms of level approximation and power when the sample sizes are balanced. For unbalanced designs, the asymptotic normal approximation lacks accuracy in terms of level (see the detailed discussion in Section 4). Related to this topic, one can mention the extension to testing for ordered alternatives, when k > 2. In the unconditional case, this problem has been considered, for instance, in Singh et al. (2009), Noguchi and Gel (2010) and Pallmann et al. (2014). Several alternatives for the variance functions may be considered, including an ordered trend

8

G. Boente and J.C. Pardo–Fernández / Statistics and Probability Letters 156 (2020) 108612

(σ1 ≤ σ2 ≤ . . . ≤ σk ), umbrella type (σ1 ≤ σ2 ≤ . . . ≤ σℓ and σℓ ≥ σℓ+1 ≥ · · · ≥ σk ) or tree type (σ1 ≤ σ2 , σ1 ≤ σ3 , . . . , σ1 ≤ σk ). The extension of our method to each of these possibilities is not immediate as a careful analysis for the construction of the test statistic would be required. We leave this interesting problem for future research. Acknowledgements The authors wish to thank two anonymous referees for valuable comments which led to an improved version of the original paper. This research was partially supported by Grants 20020170100022ba from the Universidad de Buenos Aires and pict 2014-0351 from anpcyt, Argentina (G. Boente), and by Grants MTM2016-76969P (G. Boente) and MTM201789422-P (J. C. Pardo-Fernández) from the Ministry of Economy, Industry and Competitiveness, Spain (MINECO/AEI/FEDER, UE). The research began while G. Boente was visiting the University of Vigo. Appendix A. Supplementary data Supplementary material related to this article can be found online at https://doi.org/10.1016/j.spl.2019.108612. References Boente, G., Pardo-Fernández, J.C., 2016. Robust testing for superiority between two regression curves. Comput. Statist. Data Anal. 97, 151–168. Cortez, P., Cerdeira, A., Almeida, F., Matos, T., Reis, J., 2009. Modelling wine preferences by data mining from physicochemical properties. Decis. Support Syst. 47, 547–553. Fan, J., Gijbels, I., 1996. Local Polynomial Modelling and its Applications. Chapman & Hall. Koul, H.L., Schick, A., 2003. Testing for superiority among two regression curves. J. Statist. Plann. Inference 117, 15–33. Neumeyer, N., Dette, H., 2005. A note on one–sided nonparametric analysis of covariance by ranking residuals. Math. Methods Statist. 14, 80–104. Neumeyer, N., Pardo-Fernández, J.C., 2009. A simple test for comparing regression curves versus one-sided alternatives. J. Statist. Plann. Inference 139, 4006–4016. Noguchi, K., Gel, Y.G., 2010. Combination of levene-type tests and a finite-intersection method for testing equality of variances against ordered alternatives. J. Nonparametr. Stat. 22, 897–913. Pallmann, P., Hothorn, L.A., Djira, G.D., 2014. A levene-type test of homogeneity of variances against ordered alternatives. Comput. Statist. 29, 1593–1608. Pardo-Fernández, J.C., Jiménez-Gamero, M.D., El Ghouch, A., 2015. Tests for the equality of conditional variance functions in nonparametric regression. Electron. J. Stat. 9, 1826–1851. Singh, P., Gill, A.N., Kumar, N., 2009. A test of homogeneity of several variances against tree ordered alternative. Statist. Probab. Lett. 79, 2315–2320.