Data-Driven Percentile Modified Two-sample Wilcoxon Test

Data-Driven Percentile Modified Two-sample
Wilcoxon Test
O. Thasa ∗and J.C.W. Raynerb
a
Dept. of Applied Mathematics, Biometrics and Process Control,
Ghent University, Belgium
b
School of Mathematical and Physical Sciences, University of
Newcastle, Callaghan, NSW 2308, Australia
∗
Corresponding author: Olivier Thas, Department of Applied Mathematics, Biometrics and
Process Control, Ghent University, Belgium (e-mail: olivier.thas@UGent.be)
1
Summary
The two-sample Wilcoxon rank sum statistic can be derived as the first component
of the Pearson chi-squared statistic in a particularly constructed contingency table.
For this test a “percentile modification” has been proposed, which is equivalent to
splitting the contingency table into two independent subtables, and computing the
Wilcoxon statistic on one of the subtables. Although this procedure does not use
all data in the sample, it often results in a power increase. The splitting position
is determined by an arbitrarily chosen trimming proportion p. We first illustrate
that an inappropriate choice of p may result in substantial power loss. To circumvent this problem, we propose a new test statistic by using a data-dependent
choice for p, say pˆ. We show that its asymptotic null distribution is the supremum
of a time-transformed Brownian Motion. At the rejection of the null hypothesis,
informative conclusions may be obtained by conditioning on pˆ. In a simulation
study it is shown that our solution results in a power advantage for some particular alternatives. Also, instead of using only one subtable, we suggest to compute
the Wilcoxon statistic on both subtables, and to consider their sum as a new test
statistic, which we consider as a recomposition of statistics, rather than a decomposition.
keywords: contingency table, Pearson chi-squared, Brownian Motion
2
1
INTRODUCTION
The Wilcoxon rank sum test (Wilcoxon, 1945) is probably the most common nonparametrical test for testing the two-sample location problem
H0 : F1 (x) = F2 (x) for all x
(1)
H1 : F1 (x) = F2 (x − ∆) for all x,
(2)
versus
where ∆ 6= 0 measures the location shift. This formulation of the two-sample
location problem is generally referred to as the location-shift model. Under the
restriction that F1 and F2 have the same “shape”, the location-shift model implies
that the hypotheses are equivalent to H0 : ∆ = 0 versus H1 : ∆ 6= 0, where ∆ is
now the difference in means (medians).
Although most text books present the location-shift model to introduce the Wilcoxon
rank sum test, this test actually tests the null hypothesis H0 : P [X < Y ] =
1/2, where X and Y have distribution functions F1 and F2 , respectively, without imposing any further restrictions on their shape. Note that F1 = F2 implies
P [X < Y ] = 1/2, but the reverse is not necesarilly true. Some authors have argued (see Hollander & Wolfe, 1999, Chapter 4, for an overview) that a conclusion
as e.g. P [X < Y ] > 1/2 is more informative than a conclusion about differences
in means. For example, suppose X and Y are responses under placebo and active
treatment, respectively. Then, rejecting H0 in favour of P [X < Y ] > 0.5, implies
that it is more likely that giving an active treatment to a patient will result in a
higher response than when the patient would have been given placebo.
Some people prefer conlclusions in terms of stochastic ordenings of F1 and F2 ,
i.e. they want to test F1 = F2 against F1 (x) < F2 (x) (F1 (x) > F2 (x)) for all x.
This is referred to as “X is stochastically larger (smaller) than Y ”. Evidently, this
further implies P [X < Y ] < 1/2 (P [X < Y ] > 1/2), but this is not necessarily
true in the other direction.
As mentioned before, the most common nonparametrical test for testing H0 is the
Wilcoxon rank sum test (Wilcoxon, 1945). Rayner and Best (2001) have shown
how this statistic is obtained as the first component of Pearson’s χ2 statistic computed on a particularly constructed contingency table. This is summarized in Section 2.
1
Gastwirth (1965) has proposed a percentile modified two-sample Wilcoxon test.
Basically, his modification consists in pooling the data of the two samples, and
performing the Wilcoxon test on only a fraction of the data, where the fraction is
determined by percentages r and p of the most extreme small and large data points
in the pooled sample, respectively. In this paper, we set r = 0, i.e. the fraction
only refers to a portion of extreme large observations. This trimming proportion
p must be determined a priori. Although in this way, a part of the data is actually
discarded, Gaswirth has shown that in many occasions a power gain is obtained.
However, if p is chosen badly, then there is a risk of loosing power. In Section 3 we
show that this modification is basically the same as splitting the contingency table
and using only one part for further calculations. We also illustrate the meaning of
p and argue that an appropriate value for p may be estimated from the data, say pˆ,
so that the percentile modified Wilcoxon statistic based on pˆ is maximized. Also
this pˆ has an interesting interpretation, and the maximized Wilcoxon statistic may
be used as a test statistic. In Section 4 we give the asymptotic null distribution.
The results of a small simulation study are presented in Section 5.
2
The Contingency Table Approach
Let X1 , . . . , Xn1 and Y1 , . . . , Yn2 denote the sample observations from F1 and F2 ,
respectively, and let n = n1 + n2 , the total sample size. The pooled sample
observations are denoted by Zi , and the corresponding order statistics by Z(i)
(i = 1, . . . , n). A binary variable S is defined as an indicator for the original
sample, i.e. Si = 1 if Zi is an observation from the first sample, and Si = 2 if Zi
is originally from the second sample. Similarly, S(i) refers to Z(i) . The two-sample
Wilcoxon rank sum statistic is then given by
Wn =
n
X
jI(Z(j) ),
(3)
j=1
where
I(Z(j) ) =
if S(j) = 1
.
if S(j) = 2
1
0
(4)
A 2 × n contingency table, say {Nij }, is constructed as follows (this exposition
follows from Rayner and Best (2001)). Let N1j = I(Z(j) ) and N2j = 1 − I(Z(i) )
2
(j = 1, . . . , n). Note that by construction the n column totals N.j = N1j +
N
P2jn (j = 1, . . . , n) are fixed and all equal to one. The two row totals Ni. =
j=1 Nij = ni are the number of observations in the first (i = 1) and in the
second (i = 2) sample.
Under the two-sample null hypothesis, the two column distributions {N1j }j and
{N2j }j are equal to the marginal distribution, which is the uniform. Therefore,
Pearson’s X 2 test seems an appropriate choice of test statistic. In the contingency
table Pearson’s X 2 statistic for independence (or homogeneity of column distributions) is given by
2 X
n
X
(Nij − Ni. /n)2
2
.
Xn =
N
/n
i.
i=1 j=1
From Rayner and Best (2001), Pearson’s X 2 statistic satisfies
Xn2 =
n−1
X
n−1
X
Uj2 =
j=1
Vjt Vj ,
j=1
where Vjt = (V1j , V2j ) and
Vij =
−1/2
Ni.
n
X
Nik gj (k).
k=1
Here, the {gj } are orthonormal polynomials on the marginal column distribution,
which is, by construction, the uniform discrete distribution {n−1 , n−1 , . . . , n−1 }.
We take gn (k) = 1 for all k = 1, . . . , n, so that for s = 1, . . . , n − 1, gs (.) is a
polynomial of degree s. In particular, the degree one polynomial is
r
12
n+1
k−
.
g1 (k) =
n2 − 1
2
As the first component we find

Wn −
U12 = V1t V1 =  q
2
N1. (n+1)
2

,
(5)
N1. N2. (n+1)
12
which is exactly the squared standardized Wilcoxon rank sum statistic, which
has asymptotically a χ21 distribution under H0 . Moreover, the asymptotic null
3
distribution of U12 follows also directly from the general theory given in Rayner
and Best (2001).
The second component of Xn2 is

2
Pn
n2 −1
n+1 2
√
− 12 N1.
k=1 N1k k − 2
 ,
U22 = V2t V2 =  n q
1
2
2
N N (n − 1)(n − 4)
180 1. 2.
(6)
which is exactly the square of the standardized Mood statistic for testing equality
of dispersion. The third component U32 = V3t V3 of Xn2 , is given by

√
U32 = V3t V3 =  n
Pn
n+1
k=1 N1k k − 2
q
3
−
3n2 −7
20
Pn
k=1 N1k k −
(n2 −1)(n2 −4)(n2 −9)N1. N2.
2800
n+1
2
2
 .
(7)
They all have an asymptotic
χ21
null distribution.
[something about the relation with two-sample tests for dispersion and skewness
with reference to (Eubank, LaRiccia, & Rosenstein, 1987)]
3
3.1
The Percentile Modification
Percentile Modified Statistics
Gastwirth (1965) proposed a percentile modification to increase the power of twosample tests. He claims that differences between distributions or samples are often
masked when all observations are considered. In particular, differences become
sometimes more apparent in the tails of the distribution. Therefore, his suggestion
is to remove a fraction of the sample observations and to apply the Wilcoxon test
on the remaining observations.
The fraction of the observations to be discarded is denoted by p, and let np denote the largest integer smaller than or equal to pn. Given the original order
statistics Z(i) , we define the lower observations as Z(1) , . . . , Z(np ) , and the upper observations as Z(np +1) , . . . , Z(n) . The construction of the contingency table
{Nij }i=1,2;j=1,...,n immediately implies that the lower observations completely de4
termine the left part of the table, {Nij }i=1,2;j=1,...,np , and, similarly, the upper observations determine the right part of the table, {Nij }i=1,2;j=np +1,...,n . The row
totals are denoted by Ni.low = Ni.low (p) and Ni.up = Ni.up (p) (i = 1, 2), respectively.
For notational comfort, the dependence on p will be omitted unless it is needed.
We will continue with the contingency table corresponding to the upper observations and give these upper observations ranks 1, . . . , n − np . In this table the
theory of Section 2 applies. Hence, the first component of Pearson’s X 2 statistic
is the square of the standardized conditional Wilcoxon rank sum statistic, based
on the upper subsample. Conditional refers to the conditioning on the row totals
of the selected subtable. In particular, the conditional Wilcoxon statistic is given
by
n
X
up
up
W = W (p) =
(j − np )I(Z(j) ),
(8)
j=np +1
and the first component is the square of
T
up
W up − 21 N1.up (nq + 1)
,
= T (p) = q
up up
1
N
N
(n
+
1)
q
2.
12 1.
up
(9)
where nq = n − np is the number of upper observations in the sample. Note that
W up (0) is exactly the Wilcoxon statistic. To apply the general distribution theory
of Pearson components (Rayner & Best, 2001), it must be assured that nq → ∞
as n → ∞. Therefore, we will assume that p < 1, which is an obvious restriction
from a ‘data’-analysis point of view. Under this condition it follows immediately
that T up has asymptotically a standard normal null distribution.
It may happen that one believes that the difference between F1 and F2 will occur
in the lower part. Then, the lower contingency table may be used to construct
a statistic similar to T up . Since the construction is basically exactly the same,
we will give no further details. Let T low = T low (p) denote this statistic. Note
that also T low has an asymptotic standard normal null distribution, provided that
p > 0. Since T low and T up are calculated on mutually distinct subsamples, it may
be informative to combine them both into one test statistic,
T tot (p) = T up (p)2 + T low (p)2 .
(10)
Equation 10 looks like a decomposition of the T tot (p) statistic, though we pefer to
look at it as a recomposition of the statistics T up (p) and T low (p).
5
3.2
Profiles
The distribution theory of the percentile modified Wilcoxon test only holds when p
is fixed, i.e. when p is chosen prior to the observation of the data. Of course, it may
often happen that p is chosen inappropriately, resulting in a low power test. In this
section, we propose to compute a profile of the percentile modified Wilcoxon test.
A profile is simply T up (p) considered as a function of the trimming proportion p.
Since each p corresponds to exactly one threshold value Z(np ) , the profile can also
be defined as T up (p(z)) with p(z) the smallest p such that p(Znp ) ≤ p. In a similar
way T low (p(z)) and T tot (p(z)) are the profiles of the T low and T tot statistics.
We argue that profile plots are informative exploratory tools, but we will also use
them here to illustrate that some choices of the trimming proportion p may result
in poor tests. As an example, we consider two samples of 50 observations each.
Figure 1 shows the EDF’s of the two samples and the three profile plots for T up ,
T low and T tot . The EDF’s clearly show that the shapes of the two distributions are
quite different. The profile plots show for each z (or corresponding p) the value of
the percentile modified test statistic. For each fixed p, the level α = 0.05 critical
value of T up (p) and T low (p) is 1.96 (two-sided; asymptotic approximation). From
the profile plot of T up , it is seen that significance would only be obtained with
1 < z < 4, and, similarly, for T low significance only occurs for 2.5 < z < 4. For
each p, T tot (p) is asymptotically χ22 distributed. Thus, significance would result
from T tot (p) > 5.99. This happens when 2 < z < 4. These plots clearly illustrate
that the final results of the percentile modified tests depend heavily on the choice
of p.
The profile plots also suggest another approach: define a test statistic as the maximum absolute value of the profile over p ∈ (0, 1). This is further discussed in the
next section.
4
The Supremum Statistic
4.1
The Test Statistics and Asymptotic Null Distributions
Define
Snup = sup |Tnup (p)|,
p
6
(11)
●●
1
● ●
●
−1
●●
−5
profile
5
10
15
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●●
● ●
●●
● ●●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
0
●●●
●
●●●
5
10
z
lower
total
35
z
●●●
●
●●●
●● ● ●
5
10
15
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
0
z
●●
15
●
●●
15
●● ● ●
●●
●
●●
●
●
●
●
●
●
●
●
●●●
●
● ●
●
●
●● ●
●
●●
●
●
●
●
25
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
● ●
●
●
●
●●
●●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●●
●
●
●
●
●
0
●
●●
0 5
2
1
−1
profile
3
4
5
0
●
●
●●
−3
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●● ●
●
●
● ●
●
●
●
●
●●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
upper
profile
0.4
0.0
Fn(x)
0.8
EDF
●●●
●
●●●
5
10
●● ● ●
●●
15
z
Figure 1: Example data set. Upper left pannel: EDF’s of the two samples; upper
right pannel: profile of T up ; lower left pannel: profile of T low ; lower right pannel:
profile of T tot .
7
and denote pˆ for the value of p for which the supremum is obtained. In a similar
way Snlow = supp |Tnlow (p)| and Sntot = supp Tntot (p) are defined. Of course the
asymptotic null distribution of these new test statistics will no longer be standard
normal. In the following theorem the asymptotic null distribution is discussed. In
Figure 1 the pˆ corresponding to the maxima are indicated by vertical lines. It is
interesting to note that the positions of these lines more or less coincide with the
crossing of the two EDF’s, and that they approximately devide the z-range into
parts with conditional stochastic ordenings.
The asymptotic null distribution of Snup is given in the next theorem. The proof is
given in appendix A.
Theorem 1 Let 0 < θ < 1, Ψ(p) = (1 − p)3 , and let S(t) denote a Brownian
motion on [Ψ(θ), 1]. Then, as n → ∞
(a)
w S ◦Ψ
T up −→ 1/2
Ψ
(b)
S(t) d
up
up
Sn = sup |T (p)| −→ sup 1/2 0≤p≤θ
Ψ(θ)≤t≤1 t
Note that the restriction p ≤ θ < 1 expresses the condition that there should be
observations in the upper sample. A very similar theorem follows immediately
for T low and Snlow . Once the weak limits of T up and T low are known, the following
theorem follows almost immediately (an outline of the proof is given in appendix
B).
Theorem 2 Let 0 < θ < 12 , Ψup (p) = (1 − p)3 , Ψlow (p) = p3 , and let S(t) denote
a Brownian motion on [Ψlow (θ), Ψup (θ)]. Then, as n → ∞
(a)
!
!
T
tot
w
−→
2
S ◦ Ψup
+
1/2
Ψup
2
S ◦ Ψlow
1/2
Ψlow
(b)
Sntot =
d
sup T tot (p) −→
θ≤p≤1−θ
"
sup
θ≤p≤1−θ
8
S(Ψup (p))
Ψup (p)1/2
2
+
S(Ψlow (p))
Ψlow (p)1/2
2 #
Table 1: Simulated critical values for the Snup (Snlow ) and Sntot statistic
Sntot
n1 = n2
Snup (Snlow )
α = 0.01 α = 0.05 α = 0.10 α = 0.01 α = 0.05 α = 0.10
10
2.8284
2.4056
2.2478
9.9709
8.0404
7.1333
25
3.1720
2.7045
2.5311 12.8290
9.8143
8.7879
50
3.3526
2.8999
2.6734 14.2166 11.1685
9.7863
∞
3.7488
3.2341
3.0111
Simulations of Snup under H0 for several sample sizes n1 and n2 have indicated
that the convergence to the limiting null distribution is quite slow. Therefore, we
propose to work with the Monte Carlo approximation to the exact null distribution.
For sample sizes (n1 = n2 ) n1 = 10, n1 = 25 and n1 = 50, and levels α = 0.01,
α = 0.05 and α = 0.10, the simulated and asymptotic critical values are given in
Table 1. All estimates are based on 100,000 simulation runs.
In the example presented in Figure 1 both samples had 50 observations. The
critical value of Snlow and Snup at α = 0.05 is 2.9 (Table 1), and for Sntot we find
11.2. Using these critical values as thresholds in the profile plots of Figure 1,
we find that all three sup-tests result in a rejection of the null hypothesis H0 :
P [X < Y ] = 0.5.
4.2
Interpretation of Alternatives
As mentioned in the introduction, the Wilcoxon rank sum test tests the null hy1 +1
ˆ = n11n2 W − n2n
pothesis H0 : P [X < Y ] = 21 . It is also interesting to note that π
2
is an unbiased estimator of P [X < Y ]. This is a direct consequence of the relationship between the Wilcoxon and the Mann-Whitney statistics. Thus, at the
rejection of H0 , π
ˆ may be directly used to formulate a conclusion.
In many practical situations P [X < Y ] is close to 0.5, but a closer look at the
data suggests that this is caused by computing P [X < Y ] as an integral over the
whole support of the distributions F1 and F2 . For example, the two empirical
distribution functions (EDF) in Figure 2 show samples from two normal distributions, both with mean 2, but the data from the first group has standard deviation
1, whereas the standard deviation in the second population is 0.2. The traditional
9
1.0
EDF
0.2
0.4
Fn(x)
0.6
0.8
group 1
group 2
0.0
●
0
●
●
●
●
●
●
●
1
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
2
3
4
x
Figure 2: EDF of two samples from normal distributions with means equal to 2,
and standard deviations equal to 1 (group 1) and 2 (group 2)
two-sample Wilcoxon rank sum test gives a p-value of 0.46, so that it is concluded that P [X < Y ] = 0.5. However, Figure 2 indicates that small values of
x are more likely to result from group 1. When x < 2, the latter statement is
true with probability close to one! Also, when x > 2, it is with very high probability that the response of group 2 is smaller than the respose of group 1. This
is exactly what the percentile modified Wilcoxon test of Gastwirth (1965) could
have detected if p was chosen such that a fraction p of the pooled sample data is
smaller than 2 (p approximately 0.5). With such a p, the null hypothesis of the
percentile modified Wilcoxon test is H0 : P [X < Y |X, Y > 2] = 0.5, and for
the example data set this test rejects H0 with a p-value of zero, resulting in the
conlcusion P [X < Y |X, Y > 2] <<< 0.5. With the same trimming proportion
p, a percentile modified Wilcoxon test could have been computed on the lower instead of the upper observations. Also here the null hypothesis is strongly rejected,
resulting in the conclusion P [X < Y |X, Y < 2] >>> 0.5.
10
5
Simulation Study
We limit the simulation study to some specific alternatives which may be of particular interest. For instance, two distributions which show conditional stochastic
ordenings. As a benchmark comparison we also include some simulations under
the location-shift model. Although there are many tests described in the literature
for the two-sample problem, we only include the Wilcoxon rank sum test and the
two-sample t-test. The latter is only considered because it is often (mis)used in
practice. All tests are performed at the 5% level of significance and all powers are
estimated based on 10.000 simulation runs. Sample sizes of n1 = n2 = 20 and
n1 = n2 = 50 are considered.
As a first alternative, two normal populations in a location-shift model are considered. The first is standard normal, and the second has also variance 1 and a varying
mean from 1 to 1.5 in steps of 0.25. Figure 3 shows the EDF’s and profile plots of
a random sample of n1 = n2 = 20 observations. As a second alternative, also two
normal distributions are considered. Again, the first is standard normal, and the
second has mean 1 and a varying standard deviation from 1 to 3 in steps of 0.5.
Figure 4 shows the EDF’s and profile plots of a random sample of n1 = n2 = 20
observations. Finally, two exponential distributions are considered. The first has
rate 1, and the second has a varying rate between 1 and 3 in steps of 0.5. Since
the rate parameter influences both the mean and the shape in general, the samples
from these distributions are corrected for the mean, i.e. the theoretical mean (inverse of the rate) is substracted from each sample. An example is presented in
Figure 5.
The results of the simulation study are presented in Figures 6, 7 and 8.
For the location-shift alternative, Figure 6 shows that the best powers are obtained
with the t-test, very closely followed by the Wilcoxon test. This is as expected, for
the simulations were performed under all optimality conditions for the t-test. The
powers of the three new tests are smaller, but no power break down is observed.
The powers of the new tests are quite similar. Under the second alternative (Figure 7), the powers of the t-test and the Wilcoxon test decrease as the standard
deviation of the second normal distribution is increased. The opposite behaviour
is observed for the Snup and Sntot tests. The power of the Snlow test is generally low.
Finally, Figure 8 shows very large powers for the Snlow and Sntot tests, moderate
powers for the Snup test, low powers Wilcoxon test, and virtually no power for the
11
−1
−2
profile
−3
1
●
●●
●
●
●
●
●
●●
1
●
0
●
●
●
2
●
●
●
3
●
●
●
●
−1
0
●
●
●
●●
●
●
●
● ●●
●
● ●●
●●
●
0
1
● ●
●
●
●
2
3
●
●
●
●
●
●
●
●
●
2
3
●
●
● ●
●
2
3
●
●●
●●
●●
●
● ●
●
●
●
●
●
●
● ●●
●
●
●●
−1
z
●
●
●●
●
●
●
1
total
●
●
●
●
● ●
●
● ●
●●
● ●● ●
●
●
●
●
●
●
●●
●●
lower
●
−2
●
●
●
●
●
●
●
●
●
z
●
●●
−1
●
●
z
−3
profile
0
●●
●
●
−1
0
−1
●●
●
●
●
●
●
●
●
●
●
12
●
●
●
●
●
●
2 4 6 8
●
●
●
●
●
●
●
●
●
●
●
upper
profile
0.4
0.0
Fn(x)
0.8
EDF
0
1
z
Figure 3: Example data set (two normal distributions with variance 1, and means
0 and 1, respectively). Upper left pannel: EDF’s of the two samples; upper right
pannel: profile of T up ; lower left pannel: profile of T low ; lower right pannel:
profile of T tot .
12
−2
0
2
−0.5
0
−1.5
−2.5
total
−2
0
2
4
8
●●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
2
●
●● ●
2
●
0
0.0
−2
lower
profile
2.0
profile
1.0
−4
−4
z
●
●
●
●
z
●
−6
●
●● ●
●
−6
●
● ●
●
●●
●● ● ●
●●●
●
● ●
●
● ●● ●
●
●
● ●●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●●●●
●●● ●
●●
●● ●
●
4
●
●
●
●
●
●
10
−4
●
●
6
−6
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●●
●●
●●
● ●
● ●
● ●
●
●
●
4
●
●
●
upper
profile
0.4
0.0
Fn(x)
0.8
EDF
4
●
●
●
−6
z
−4
●
●
●
●● ●
−2
0
2
●
4
z
Figure 4: Example data set (two normal distributions with means 0 and 1, and
standard deviations 1 and 3, respectively). Upper left pannel: EDF’s of the two
samples; upper right pannel: profile of T up ; lower left pannel: profile of T low ;
lower right pannel: profile of T tot .
13
EDF
3
2
1
0
0.5
●● ●●
●
●
●
●
●
●
●
−0.5
20
total
●
0.0
10
●
●
0.5
●
●●
●
●●
●
● ●●
●
●● ●●
−0.5
z
0.5
●
5
●
●●
●●
●
●●●●
●
●
●
●
●
●
●
●
●
●
●
15
●
●
0.0
lower
0
0
−2
−3
●
● ●●
●
z
●
●
●
●
●
●
●
●
●
●● ●
●● ●
●●
●
●
−4
●
z
●
●
●
−0.5
●
●
●
●
0.0
●● ●●
−1
●
● ●●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
profile
0.4
●
●
●
●
●
●
●
−0.5
profile
●
profile
●
●
● ●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
0.0
Fn(x)
0.8
●
upper
●
0.0
●
0.5
z
Figure 5: Example data set (two exponential distributions with rates 1 and 3, both
corrected for the mean). Upper left pannel: EDF’s of the two samples; upper
right pannel: profile of T up ; lower left pannel: profile of T low ; lower right pannel:
profile of T tot .
14
1.0
Normal, location shift, n=50
1.0
Normal, location shift, n=20
●
●
●
●
●
0.8
0.8
●
0.6
power
power
●
0.4
●
0.2
0.4
S low
S up
S tot
t
W
0.2
0.6
●
●
●
S low
S up
S tot
t
W
●
●
●
0.0
●
0.0
●
0.0
0.5
1.0
1.5
0.0
0.5
1.0
mu2
1.5
mu2
Figure 6: Estimated powers for the location-shift alternatives for sample sizes
n = 20 (left pannel) and n = 50 (right pannel)
1.0
Normal, n=50
1.0
Normal, n=20
●
0.8
0.8
●
●
0.4
0.6
●
power
●
S low
S up
S tot
t
W
●
0.4
power
0.6
●
●
●
●
●
0.2
0.2
●
●
●
S low
S up
S tot
t
W
●
0.0
0.0
●
1.0
1.5
2.0
2.5
3.0
1.0
sd
1.5
2.0
2.5
3.0
sd
Figure 7: Estimated powers for the normal scale alternatives for sample sizes
n = 20 (left pannel) and n = 50 (right pannel)
15
Exponential, n=50
1.0
1.0
Exponential, n=20
●
●
●
●
●
●
●
●
0.6
S low
S up
S tot
t
W
0.4
0.4
power
●
power
●
0.6
●
0.8
0.8
●
●
0.2
●
●
1.0
S low
S up
S tot
t
W
0.0
0.0
0.2
●
1.5
2.0
2.5
3.0
1.0
rate
1.5
2.0
2.5
3.0
rate
Figure 8: Estimated powers for the exponential alternatives for sample sizes n =
20 (left pannel) and n = 50 (right pannel)
t-test.
In general, for the alternatives considered, the Sntot test seems to be a reasonable
compromise between the Snlow and Snup tests.
16
APPENDIX: PROOF OF THEOREM 1
First, the weak convergence of T up (p) is shown. Part (b) of the theorem is a direct
consequence of part (a) and the continuous mapping theorem.
To proof the weak convergence of T up (p), we will first start with showing the
weak convergence of the statistic
√
√
1 up
1 up
up
Sn (p) = 48σn (p)T (p) = n
W (p) − 2 N1. (n(1 − p) + 1) ,
n2
2n
q
1
N up (p)N2.up (p)(n(1 − p) + 1). This proof consists of two
where σn (p) =
12n3 1.
parts: we need to show (1) finite dimensional convergence to a multivariate normal distribution, and (2) tightness of {Sn }. Once weak convergence of Sn is
established, we use the Skorohod construction and a strong approximation result
for the weighted version T up of Sn .
First, the notation is simplified. Let Nn (p) = N1.up (p) and Wn (p) = W up (p).
(1)
To prove multivariate normality, it is sufficient to apply the multivariate central
limit theorem. It has been shown before that for each p, E {Sn (p)} = 0. Furthermore, as n → ∞, it is easy to show that Var {Sn (p)} = 48σn2 (p) → (1 − p)3 (we
have used limn→∞ Nnn(p) = 12 (1 − p), which is the convergence of a non-random
series, for all inference is conditional on Nn (p)). We only need to calculate the covariance Cov {Sn (p), Sn (q)}. Straightforward, though lengthy calculations give
lim Cov {Sn (p), Sn (q)} = (1 − p)3 ∧ (1 − q)3 ,
(12)
n→∞
where ∧ denotes the minimum operator. Equation 12 shows exactly the covariance
function of a Brownian motion in transformed time (1 − p)3 .
(2)
Tightness is proven if we succeed in showing that for all n ≥ 1, and p1 ≤ p ≤ p2 ,
E (Sn (p1 ) − Sn (p))2 (Sn (p) − Sn (p2 ))2 ≤ (µ(p) − µ(p1 ))(µ(p2 ) − µ(p)),
where µ is a continuous measure on [0, 1]. This is essentially Theorem 6 in Chapter 2 of Shorack and Wellner (1986), with a = 1 and b = 2.
17
By the independence of Sn (p1 ) − Sn (p) and Sn (p) − Sn (p2 ), we have
E (Sn (p1 ) − Sn (p))2 (Sn (p) − Sn (p2 ))2
= E (Sn (p1 ) − Sn (p))2 E (Sn (p) − Sn (p2 ))2 .
Further, we have
E (Sn (p1 ) − Sn (p))2 = Var {Sn (p1 )} − 2Cov {Sn (p1 ), Sn (p)} + Var {Sn (p)}
1
1
1
≤
(1 − p1 )3 − 2 (1 − p)3 + (1 − p)3
4
4
4
3
3
≤ (p − 1) − (p1 − 1)
= µ(p) − µ(p1 ),
where µ(p) = (1 − p)3 is a continuous measure on [0, 1]. The first inequality sign
in the above calculations follow from the observation that the variance of Sn (p) is
a decreasing function of n. Similarly, E (Sn (p) − Sn (p1 ))2 ≤ (p2 − 1)3 − (p −
1)3 = µ(p2 ) − µ(p). Hence,
E (Sn (p1 ) − Sn (p))2 (Sn (p) − Sn (p2 ))2 ≤ (µ(p)−µ(p1 ))(µ(p2 )−µ(p)). (13)
This completes the proof of the weak convergence of Sn to a Brownian motion in
transformed time (1 − p)3 .
For the Skorohod construction the following strong approximation holds
Sn (p) − S(Ψ(p)) a.s.
−→ 0 as n → ∞,
(14)
sup Ψ1/2 (p)
0≤p≤θ
Rθ 1
provided that 0 Ψ(t)
dΨ(t) = − ln(Ψ(θ)) < ∞, which holds because θ < 1 is
assumed. This strong approximation result may be reformulated as
3 up
S((1
−
p)
))
a.s.
−→
sup T (p) −
0 as n → ∞,
(15)
3/2
(1 − p)
0≤p≤θ
from which the weak convergence of T up immediately follows.
APPENDIX: PROOF OF THEOREM 2
w
From Theorem 1 we know that T up −→
S◦Ψup
1/2
Ψup
S◦Ψlow
1/2 , as
Ψlow
T tot and Sntot are
n → ∞. Both part (a) and part (b) follow immediately, because
continuous function of T low and T up (continuous mapping theorem).
18
w
, and similarly, T low −→
References
Eubank, R., LaRiccia, V., & Rosenstein, R. (1987). Test statistics derived as components of Pearson’s phi-squared distance measure. Journal of the American Statistical Association, 82, 816-825.
Gastwirth, J. (1965). Percentile modification of two sample rank tests. Journal of
the American Statistical Association, 18, 1127-1141.
Hollander, M., & Wolfe, D. (1999). Nonparametric statistical methods. New
York: Wiley.
Rayner, J., & Best, D. (2001). A contingency table approach to nonparametric
testing. New York: Chapman and Hall.
Shorack, G., & Wellner, J. (1986). Empirical processes with applications to
statistics. New York, USA: Wiley.
Wilcoxon, F. (1945). Individual comparisons by ranking methods. Biometrics
Bulletin, 1, 80-83.
19