- Email: [email protected]
Estimating the covariance of bivariate order statistics with applications
Statistics & Probability Letters 48 (2000) 195 – 203
Estimating the covariance of bivariate order statistics with applications Alan D. Hutson1 Division of Biostatistics, Department of Statistics, Health Science Center, University of Florida, P.O. Box 100212, Gainesville, FL 32610-0212, USA Received April 1999; received in revised form July 1999
Abstract In this note we develop a new analytic bootstrap-type method for estimating the covariance between the bivariate order statistics Xr:n and Ys:n . Correlation and regression estimators follow straightforward from this result. An easy-to-use method for calculating the nonparametric percentile con dence regions for bivariate quantiles is developed and illustrated using biological data. A speci c case of interest is the application to bivariate medians, i.e., correlated marginal medians. c 2000 Elsevier Science B.V. All rights reserved
Keywords: Bootstrap; Concomitants; Bivariate median; Order statistics; Quantile function
1. Introduction Let (X1 ; Y1 ); (X2 ; Y2 ); : : : ; (Xn ; Yn ) denote i.i.d. pairs of data obtained by sampling from the bivariate distribution function F(x; y) with density function f(x; y) continuous everywhere having support in R2 . Also, let X1:n 6X2:n 6 · · · 6Xn:n , Y1:n 6Y2:n 6 · · · 6Yn:n denote the marginal order statistics. Our goal is to estimate analytically the covariance and=or correlation between Xr:n and Ys:n , as well as develop a nonparametric bootstrap-type percentile con dence regions for the bivariate quantiles QX (u1 ) and QY (u2 ). A speci c case that we are interested in is the joint con dence region for the bivariate medians (QX (1=2); QY (1=2)). Note that this new method is not meant as a competitor for estimating bivariate medians, of which there are many choices, e.g. see Small (1990) for an overview of multidimensional medians, and does not have the same properties as many of the bivariate median estimators such as ane equivariance. This is a completely new nonparametric method focusing on the correlation between marginal order statistics obtained from bivariate E-mail address: [email protected] .edu (A.D. Hutson) Work was supported in part by NIH General Clinical Research Center Grant RR0082.
1
c 2000 Elsevier Science B.V. All rights reserved 0167-7152/00/$ - see front matter PII: S 0 1 6 7 - 7 1 5 2 ( 9 9 ) 0 0 2 0 5 - 9
196
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
samples. Furthermore, the method described below is straightforward for the applied statistician to program using standard software packages. It is already known that the exact bootstrap moment estimators of E(Xr:n ), Var(Xr:n ) and Cov(Xr:n ; Xs:n ) 16r6s6n, are analytically tractable and are given by E ∗ (Xr:n ) =
n X
wj(r) Xj:n ;
(1.1)
j=1 n X
Var ∗ (Xr:n ) =
wj(r) (Xj:n − ˆr:n )2 ;
(1.2)
j=1
Cov∗ (Xr:n ; Xs:n ) =
j−1 n X X
wij(rs) (Xi:n − ˆr:n )(Xj:n − ˆs:n ) +
j=2 i=1
where
n X
vj(rs) (Xj:n − ˆr:n )(Xj:n − ˆs:n );
j−1 j n ; r; n − r + 1 − B ; r; n − r + 1 ; B wj(r) = r r n n s−r−1 X
wij(rs) = n Crs
k=0
s−r−1 k
(−1)s−r−1−k s−k −1
vj(rs) = n Crs
k=0
−
j−1 n
s−r−1 k
(−1)s−r−1−k s−k −1
(1.4)
" s−k−1 # s−k−1 i i−1 − n n
j−1 j ; k + 1; n − s + 1 − B ; k + 1; n − s + 1 ; × B n n s−r−1 X
(1.3)
j=1
(1.5)
j j−1 B ; s; n − s + 1 − B ; s; n − s + 1 n n
s−k−1 j−1 j ; k + 1; n − s + 1 − B ; k + 1; n − s + 1 ; B n n
(1.6)
Rx B(x; a; b) = 0 t a−1 (1 − t)b−1 dt is the incomplete beta function and n Crs = n!=[(r − 1)!(s − r − 1)!(n − s)!]; see Hutson and Ernst (2000) for details. The estimators E ∗ (Xr:n ); Var ∗ (Xr:n ) and Cov∗ (Xr:n ; Xs:n ) in conjunction with the new estimator Cov∗ (Xr:n ; Ys:n ) de ned in the next section provide a complete picture of the covariance structure of the vector of order statistics, (X1:n ; X2:n ; : : : ; Xn:n ; Y1:n ; Y2:n ; : : : ; Yn:n ). Maritz (1991) has proposed an exact bootstrap estimator for E(Xi:n Yj:n ) based upon an estimate of the joint ˆ y) for F(x; y) in H (x; y) = P(Xr:n ¡ x; Ys:n ¡ y), c.d.f. for the bivariate order statistics given by replacing F(x; as given in David (1981), where Â1 = F(x; y); Â2 = F(x; +∞) − F(x; y); Â3 = F(+∞; y) − F(x; y); Â4 = 1 − F(x; +∞) − F(+∞; y) + F(x; y); H (x; y) =
n X n X u X i=r j=s k=t
n!  k  i−k Â3j−k Â4n−i−j+k ; k!(i − k)!(j − k)!(n − i − j + k)! 1 2
(1.7)
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
197
t = max(0; i + j − n) and u = min(i; j). The estimated probability mass at (x; y) is then given by Hˆ (x+; y+) − Hˆ (x−; y+) − Hˆ (x+; y−) + Hˆ (x−; y−);
(1.8)
ˆ y) for F(x; y) obtained as per the prescription of Maritz (1991) by substituting the empirical estimator F(x; in (1.7), where x+ implies x + and x− implies x − ( small). The approach of Maritz poses an obvious problem if the data are highly correlated such that i equals [i], where the pair (Xi:n Y[i]:n ); i = 1; 2; : : : ; n, denote the bivariate sample data ordered by X , where Y[i]:n denotes the corresponding ith concomitant statistic; see David (1981) for a description of concomitants. In the case i = [i] it is easy to see that Hˆ (x; y) is 0 for all possible (x; y) combinations. Later on we will illustrate that even in the cases where Hˆ (x; y) is nonzero for some (x; y) pairs the estimate of E(Xi:n Yj:n ) based upon (1.8) has some shortcomings. As an alternative, we propose an easy to calculate bootstrap estimator of the joint probability distribution for ˆ r:n ; Ys:n ) = 1=n Xr:n and Ys:n based upon a simple rank transformation of the joint density function given by f(X when s = [r], and 0, elsewhere. The transformation maintains the appropriate marginal weights given by (1.4) such that the standard bootstrap estimators of the marginal expectations and variances remain unchanged while leading to a new estimator of covariance and=or correlation of Xr:n and Ys:n di erent from that based upon (1.8). Details of the calculations are provided in Section 2. In Section 3 we outline a method for constructing the minimum area joint percentile con dence regions for the bivariate quantiles QX (u1 ) and QY (u2 ) and compare it to the asymptotic con dence ellipse approach and a Bonferroni correction-type approach using a biostatistical example speci cally examining bivariate medians. Furthermore, some crude simulations are carried out in order to examine the coverage probabilities of the the minimum area joint percentile con dence region for bivariate medians. Section 4 brie y describes an approach to order statistic regression. 2. Density and moment estimation In this section we derive the estimated joint probability density function for Xr:n and Ys:n , the bootstrap-type l m Ys:n ), as well as a correlation estimator for bivariate order statistics. Towards this end, let estimator of E(Xr:n (xk:n ; y[k]:n ) denote the kth observed concomitant sample pair, k = 1; 2; : : : ; n ordered by x with corresponding sample concomitant y[k]:n . Let k −1 i−1 ; ; rx1 (i; k) = max FUr:n n n k i ; ; rx2 (i; k) = min FUr:n n n [k] − 1 j−1 ; ; ry1 (j; k) = max FUs:n n n [k] j ; ; ry2 (j; k) = min FUs:n n n
k −1 k [k] − 1 [k] i i−1 j j−1 − ; FUs:n − ; − FUr:n ; − FUs:n b(i; j; k) = min FUr:n n n n n n n n n
for all i; j = 1; 2; : : : ; n, where FUr:n (·) denotes the cumulative distribution function of a beta random variable with parameters r − 1 and n − r + 1. Note that FUr:n (·) is readily calculated using most standard software packages.
198
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
The joint empirical density function for the bivariate order statististics Xr:n and Ys:n evaluated at the point (x; y); (x; y) ∈ {(xi−1:n ; xi:n ] × (yj−1:n ; yj:n ]}, is then given by fˆXr:n Ys:n (x; y) = n
n X
{rx2 (i; k) − rx1 (i; k)}{ry2 (j; k) − ry1 (j; k)}1(b(i; j; k)¿0)
(2.1)
k=1
for all i; j = 1; 2; : : : ; n, where 1(b(i; j; k)¿0) denotes the indicator function. The estimator follows straightforward by rst noting that by de nition the classical empirical density funcˆ r:n ; Ys:n ) = 1=n when s = [r], and 0, elsewhere, has uniform mass in the region (Xr−1:n ; Xr:n ] × tion estimator f(X ˆ r:n ; Ys:n ) to fˆX Y (x; y) is then a matter of subdi(Y[r]−1:n ; Y[r]:n ]. Transforming the empirical estimator f(X r:n s:n viding the mass in each region (Xr−1:n ; Xr:n ] × (Y[r]−1:n ; Y[r]:n ] under the constraint that the marginal totals fˆXr:n (x) = wj (r) for x ∈ (Xj−1:n ; Xj:n ] and fˆYs:n (y) = wj (s) for y ∈ (Yj−1:n ; Yj:n ], where wj (·) is de ned at (1.4), i.e, the marginal empirical distribution functions for Xr:n and Ys:n given by FUr:n (Fˆ X (x)) and FUs:n (Fˆ Y (y)), where Fˆ is the classical estimator, de nes a discrete “grid” which intersects the “grid” given by Fˆ X (x) and Fˆ Y (y), each with steps at the order statistics and mass determined by the intersections of the “grids” as given by (2.1). The numerical eciency of our approach relative to Maritz’s approach can be seen by examing (1.7) and (1.8). His analytical bootstrap approach requires 4un2 computations to estimate the point probability as compared to our approach which requires only n computations. Note that our de nition yields a slightly di erent estimator of fˆXr:n Ys:n (x; y) as compared to (1.8), is symmetric in terms of the correlation structure (positive=negative), always produces the appropriate marginal weights and reduces to the standard density estimator when r = s = n = 1. The covariance estimator follows straightforward and is given by the following: l m Ys:n ) is given by Theorem 2.1. The exact bootstrap-type moment estimator of E(Xr:n l m Ys:n )= E ∗ (Xr:n
n X n X i=1 j=1
l m ˆ xi:n yj:n fXr:n Ys:n (xi:n ; yj:n );
(2.2)
where fˆXr:n Ys:n (x; y) is deÿned at (2:1). Proof. First note that given that the sample order statistic X = (X1:n ; X2:n ; : : : ; Xn:n ) is held xed, then a nonparametric bootstrap replication is generated by taking a sample of size n with replacement from X . This is equivalent to generating a random sample of size n from a uniform(0; 1) distribution, with corresponding order statistic U = (U1:n ; U2:n ; : : : ; Un:n ), and applying the sample quantile function −1 ˆ Q(u) = Fˆ (u) = Xbnuc+1:n
(2.3)
ˆ to each element of U , where 0 ¡ u ¡ 1 and b·c denotes the oor function, i.e., Q(u) = Xi:n in the region given by (i − 1)=n6u ¡ i=n; i = 1; 2; : : : ; n. Then let ZZ l m Ys:n )= xl ym fXr:n Ys:n (x; y) d x dy; (2.4) E(Xr:n where fXr:n Ys:n (x; y) denotes the joint density of Xr:n and Ys:n . Rewriting (2.4) in terms of the quantile functions for Xr:n and Ys:n , respectively, yields Z 1Z 1 fXr:n Ys:n {QXr:n (u); QYs:n (v)} l m du dv; Ys:n )= [QXr:n (u)]l [QYs:n (v)]m E(Xr:n fXr:n {QXr:n (u)}fYs:n {QYs:n (v)} 0 0
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
=
n X n Z X i=1 j=1
×
FUr:n (i=n)
FUr:n ((i−1)=n)
Z
FUs:n ( j=n)
FUs:n (( j−1)=n)
199
[QXr:n (u)]l [QYs:n (v)]m
fXr:n Ys:n {QXr:n (u); QYs:n (v)} du dv: fXr:n {QXr:n (u)}fYs:n {QYs:n (v)}
(2.5)
Substituting the empirical estimators Qˆ Xr:n (u) = Qˆ X (QUr:n (u)), Qˆ Ys:n (u) = Qˆ Y (QUs:n (u)), fˆXr:n {QXr:n (u)} = FUr:n (i=n) − FUr:n ((i −1)=n); fˆYs:n {QYs:n (v)}=FYs:n (j=n)−FUs:n ((j −1)=n) and fˆXr:n Ys:n (x; y) at (2.1) into (2:5) and noting that each of the empirical estimators is constant over the regions of integration yields the result at (2.2). The proof follows similarly to what has been done previously by Hutson and Ernst (2000) in the marginal case. Corollary 2.2. The exact bootstrap-type moment estimator of the correlation of Xr:n and Ys:n denoted (Xr:n ; Ys:n ) is given by ∗ (Xr:n ; Ys:n ) =
E ∗ {(Xr:n − E ∗ (Xr:n ))(Ys:n − E ∗ (Ys:n ))} p ; Var ∗ (Xr:n )Var ∗ (Ys:n )
(2.6)
where the expectations E ∗ (Xr:n ) and E ∗ (Ys:n ) are calculated using (1:1) and the variances Var ∗ (Xr:n ) and Var ∗ (Ys:n ) are calculated using (1:2). Proof. Follows straightforward from Theorem 2.1. The following set of simple Examples contrasts are method with the method of Maritz. Examples. Consider for illustration purposes Example 1b given in Maritz (1991). Let n = 5 observations, (xi:n ; y[i]:n ); i = 1; 2; : : : ; 5, be (1; 2); (2; 3); (3; 1); (4; 4); (5; 5): Then from (1.1) and (1.2) E ∗ (X3:5 )=3; E ∗ (Y3:5 )=3; Var ∗ (X3:5 )=0:9824, Var ∗ (Y3:5 )=0:9824. From Theorem 2.1 and Corollary 2.2 we then have E ∗ (X3:5 Y3:5 )=9:55 and ∗ (X3:5 ; Y3:5 )=0:563. Using the estimator Hˆ (x; y) in conjunction with the point probabilities given by (1.8) we get E ∗ (X3:5 Y3:5 ) = 8:64 and ∗ (X3:5 ; Y3:5 ) = −0:369. Note that the Pearson correlation for the above dataset is ˆ = 0:70 indicating that a positive estimate for ∗ (X3:5 ; Y3:5 ) might be more accurate. Now let the n = 5 observations be (1; 1); (2; 2); (3; 3); (4; 4); (5; 5); then we have E ∗ (X3:5 Y3:5 ) = 8:20 and ∗ (X3:5 ; Y3:5 ) = 0:818. In the case of the Maritz estimator the joint empirical distribution function Hˆ (x; y) does not yield a proper joint empirical density estimator. Therefore, the covariance estimator in this case does not make sense. Note that the observations do not necessarily have to be perfectly correlated for this to happen, only that the ranks of the observations turn out to be similar to what happens in this example. In addition, there appears to be some inherent asymmetry in the approach developed by Maritz (1991) with regards to positively and negatively correlated data, i.e. using the estimator based upon (1.8) ∗6= − ∗ corresponding to the case when y is changed to −y. We conjecture that some type of continuity correction applied to Hˆ (x; y) might provide a better estimator of (Xr:n ; Ys:n ), which will yield results that are close to the results obtained by using our estimator of ∗ (Xr:n ; Ys:n ) de ned at (2.6).
200
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
3. Percentile conÿdence regions Hutson (1999) showed that an accurate analytical (1 − ) × 100% bootstrap-type con dence interval for FX−1 (u) = QX (u); 0 ¡ u ¡ 1, is given by (Qˆ X (QUn0 u:n ( =2)); Qˆ X (QUn0 u:n (1 − =2)));
(3.1)
where QUn0 u:n (1 − =2) denotes the quantile function for ÿ(n0 u; n0 (1 − u)) random variate and by de nition Qˆ X (u) = (1 − )Xbn0 uc:n + Xbn0 uc+1:n ;
(3.2)
where = n0 u − bn0 uc; n0 = n + 1 and b·c denotes the oor function. An obvious candidate for the joint con dence region for the bivariate quantiles (QX (u1 ); QY (u2 )) is to apply the Bonferroni correction to in (3.1) and then de ning the joint con dence region to be rectangular. In this section we show that a tighter nonparametric percentile con dence region can be obtained by trimming fˆXr:n Ys:n (x; y) at (2.1) in a speci c way and therefore taking advantage of the correlation structure between Xr:n and Ys:n in the con dence region estimation procedure. For this estimation procedure de ne Qˆ X (u1 ) = xbnu1 c+1:n and Qˆ Y (u2 ) = ybnu2 c+1:n . The algorithm for calculating the nonparametric percentile con dence region is as follows: 1. Eliminate all pairs (xi:n ; yj:n ) over the i, j combinations where fˆXr:n Ys:n (xi:n ; yj:n ) = 0. 2. Calculate the distance (di; j ) between the point estimate ranks (bnu1 c + 1; bnu2 c + 1) to each of the ranks (i; j) for the pairs remaining from step (1). 3. If 0 ¡ fˆXr:n Ys:n (xi:n ; yj:n ) ¡ then calculate the distance-to-probability ratio ri; j = di; j = fˆXr:n Ys:n (xi:n ; yj:n ). 4. Sort the remaining (xi:n ; yj:n ) pairs by ri; j in descending fashion, i.e. largest to smallest. 5. Trim the sorted pairs based upon each successive ri; j one by one and subtract the corresponding probability fˆXr:n Ys:n (xi:n ; yj:n ) from the total probability of 1 contained in the region. 6. Repeat step (5) until the reduction of the probability from 1 is near or equal to (1 − ). 7. The approximate (1 − ) × 100% percentile con dence region is then given by the convex hull generated from the remaining points, such that (Qˆ X (u1 ); Qˆ Y (u2 )) is an interior point of the region. Note that alternative de nitions where considered and may provide smaller con dence regions for speci c cases. The method de ned above was a trade-o between ease-of-use in terms of programming relative to generating smallest possible con dence region. Furthermore, a continuity correction may be applied to the process by using the linear interpolation quantile function estimator at (3.2) in conjunction with fˆXr:n Ys:n (x; y) at (2.1) and setting r = (n + 1)u1 and s = (n + 1)u2 . Another possibility includes appealing to the joint asymptotic normality of the bivariate quantiles (Qˆ X (u1 ); ˆ QY (u2 )) and generating a con dence ellipse based upon the moment estimators de ned in the previous section and given by the following equation: (x − E ∗ )(y − EY∗ ) (y − EY∗ )2 (x − EX∗ )2 − 2∗ p X ∗ + = (1 − ∗2 ) 22; ; ∗ Var X Var ∗Y Var X Var ∗Y
(3.3)
where EX∗ = E ∗ (Qˆ X (u1 )); Var ∗X = Var ∗ (Qˆ X (u1 )); EY∗ = E ∗ (Qˆ Y (u2 )); Var ∗Y = Var ∗ (Qˆ Y (u2 )); ∗ = ∗ (Qˆ X (u1 ); Qˆ Y (u2 )) and 22; denotes the upper th quantile of a chi-square distribution with 2 degrees of freedom. Example. The following example is taken from an echocardiographic imaging (EI) study, Geiser et al. (1998). The data consists of the fractional area change (FAC) from n = 101 subjects computed as AED − AES ×100%; (3.4) AED where AED and AES are the endocardial areas at end diastole and end systole, respectively. The dataset plotted in Fig. 1 contains the FACs as measured by two distinct echocardiographers (Expert 1(X ) and Expert FAC =
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
201
Fig. 1. Three types of 95% con dence regions for bivariate medians.
2(Y )), the three di erent types of con dence regions overlaying the data and a separate blown-up plot of the con dence regions. The dotted line represents the con dence ellipse, the solid line outline the rectangular shape denotes the Bonferroni region and the solid line outlining the convex hull denotes the nonparametric percentile con dence region for (QX (1=2); QY (1=2)). The bootstrap moment estimates for this dataset are as follows: E ∗ (Qˆ X (1=2)) E ∗ (Qˆ Y (1=2)) Var ∗ (Qˆ X (1=2)) Var ∗ (Qˆ Y (1=2)) ∗ (Qˆ X (1=2); Qˆ Y (1=2)) 43.76
45.29
3.65
1.80
0.86
202
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203 Table 1 Simulation results n
0.0
0.5
0.9
25 51 101
92=100 92=100 90=100
91=100 95=100 96=100
84=100 92=100 92=100
We see from Fig. 1 that the con dence ellipse and the nonparametric con dence regions do not di er substantially and provide a much more precise region as compared to the Bonferroni-type approach to the problem in terms of having a region with the minimal area that captures the correlation structure of the data. Preliminary simulation results: We carried out a crude simulation study for bivariate medians for data sampled from a standard bivariate normal distribution with correlation . Due to the complexity of programming such a simulation we resorted to examing the coverage probabilities for the 95% regions for (QX (1=2); QY (1=2)) via inspection of 100 simulations for samples of size n = 25; 51; 101 and = 0:0; 0:5; 0:9. The coverage results are provided in Table 1. From Table 1 we see that for this basic simulation study the coverage appears to be reasonable for samples of size n = 50 or larger. Obviously, further studies need to be carried out before more general statements can be made about the coverage error of the nonparametric region given alternative underlying distributions.
4. Order statistics regression Suppose that Xr:n and Ys:n are linked by the linear regression relation Ys:n = ÿ0 + ÿ1 Xr:n + :
(4.1)
Then it follows from the theory of least squares that nonparametric bootstrap-type estimators of ÿ0 and ÿ1 are given by minimizing n n X X i=1 j=1
(yj:n − (ÿ0 + ÿ1 xi:n ))2 fˆXr:n Ys:n (xi:n ; yj:n );
(4.2)
with respect to ÿ0 and ÿ1 , where fˆXr:n Ys:n (x; y) is de ned at (2.1). Note that the regression of Ys:n on Xr:n has unique solution and goes through the point (E ∗ (Xr:n ); E ∗ (Ys:n )). The solution to (4.2) is easily carried out using standard software packages via weighted least-squares regression. A special case of (4.1) of interest is the regression of the rth concomitant, Y[r]:n , on the rth order statistic Xr:n , e.g. see David (1981). The approach is easily generalized to a multivariate order statistics regression setup by considering (1.3) in conjunction with (2.2).
Acknowledgements We are grateful to the referees for carefully reading this manuscript and o ering helpful suggestions. Alan Hutson’s work was supported in part by NIH General Clinical Research Center Grant RR0082.
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
203
References David, H.A., 1981. Order Statistics, 2nd Edition. Wiley, New York. Geiser, E.A., Wilson, D.C., Wang, D.X., Conetta, D.A., Murphy, J.D., Hutson, A.D., 1998. Autonomous epicardial and endocardial boundary detection in echocardiography short-axis images. J. Amer. Soc. Echocardiography 11, 338–348. Hutson, A.D., 1999. Calculating nonparametric con dence intervals for quantiles using fractional order statistics. J. Appl. Statist. 26, 339–349. Hutson, A.D., Ernst, M.D., 2000. The exact bootstrap mean and variance of an L-estimator. J. Roy. Statist. Soc. Ser. B 62, 1–6. Maritz, S., 1991. Estimating the covariance matrix of bivariate medians. Statist. Probab. Lett. 12, 305–309. Small, C.G., 1990. A survey of multidimensional medians. Internat. Statist. Inst. Rev. 58, 262–277.
Estimating the covariance of bivariate order statistics with applications Alan D. Hutson1 Division of Biostatistics, Department of Statistics, Health Science Center, University of Florida, P.O. Box 100212, Gainesville, FL 32610-0212, USA Received April 1999; received in revised form July 1999
Abstract In this note we develop a new analytic bootstrap-type method for estimating the covariance between the bivariate order statistics Xr:n and Ys:n . Correlation and regression estimators follow straightforward from this result. An easy-to-use method for calculating the nonparametric percentile con dence regions for bivariate quantiles is developed and illustrated using biological data. A speci c case of interest is the application to bivariate medians, i.e., correlated marginal medians. c 2000 Elsevier Science B.V. All rights reserved
Keywords: Bootstrap; Concomitants; Bivariate median; Order statistics; Quantile function
1. Introduction Let (X1 ; Y1 ); (X2 ; Y2 ); : : : ; (Xn ; Yn ) denote i.i.d. pairs of data obtained by sampling from the bivariate distribution function F(x; y) with density function f(x; y) continuous everywhere having support in R2 . Also, let X1:n 6X2:n 6 · · · 6Xn:n , Y1:n 6Y2:n 6 · · · 6Yn:n denote the marginal order statistics. Our goal is to estimate analytically the covariance and=or correlation between Xr:n and Ys:n , as well as develop a nonparametric bootstrap-type percentile con dence regions for the bivariate quantiles QX (u1 ) and QY (u2 ). A speci c case that we are interested in is the joint con dence region for the bivariate medians (QX (1=2); QY (1=2)). Note that this new method is not meant as a competitor for estimating bivariate medians, of which there are many choices, e.g. see Small (1990) for an overview of multidimensional medians, and does not have the same properties as many of the bivariate median estimators such as ane equivariance. This is a completely new nonparametric method focusing on the correlation between marginal order statistics obtained from bivariate E-mail address: [email protected] .edu (A.D. Hutson) Work was supported in part by NIH General Clinical Research Center Grant RR0082.
1
c 2000 Elsevier Science B.V. All rights reserved 0167-7152/00/$ - see front matter PII: S 0 1 6 7 - 7 1 5 2 ( 9 9 ) 0 0 2 0 5 - 9
196
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
samples. Furthermore, the method described below is straightforward for the applied statistician to program using standard software packages. It is already known that the exact bootstrap moment estimators of E(Xr:n ), Var(Xr:n ) and Cov(Xr:n ; Xs:n ) 16r6s6n, are analytically tractable and are given by E ∗ (Xr:n ) =
n X
wj(r) Xj:n ;
(1.1)
j=1 n X
Var ∗ (Xr:n ) =
wj(r) (Xj:n − ˆr:n )2 ;
(1.2)
j=1
Cov∗ (Xr:n ; Xs:n ) =
j−1 n X X
wij(rs) (Xi:n − ˆr:n )(Xj:n − ˆs:n ) +
j=2 i=1
where
n X
vj(rs) (Xj:n − ˆr:n )(Xj:n − ˆs:n );
j−1 j n ; r; n − r + 1 − B ; r; n − r + 1 ; B wj(r) = r r n n s−r−1 X
wij(rs) = n Crs
k=0
s−r−1 k
(−1)s−r−1−k s−k −1
vj(rs) = n Crs
k=0
−
j−1 n
s−r−1 k
(−1)s−r−1−k s−k −1
(1.4)
" s−k−1 # s−k−1 i i−1 − n n
j−1 j ; k + 1; n − s + 1 − B ; k + 1; n − s + 1 ; × B n n s−r−1 X
(1.3)
j=1
(1.5)
j j−1 B ; s; n − s + 1 − B ; s; n − s + 1 n n
s−k−1 j−1 j ; k + 1; n − s + 1 − B ; k + 1; n − s + 1 ; B n n
(1.6)
Rx B(x; a; b) = 0 t a−1 (1 − t)b−1 dt is the incomplete beta function and n Crs = n!=[(r − 1)!(s − r − 1)!(n − s)!]; see Hutson and Ernst (2000) for details. The estimators E ∗ (Xr:n ); Var ∗ (Xr:n ) and Cov∗ (Xr:n ; Xs:n ) in conjunction with the new estimator Cov∗ (Xr:n ; Ys:n ) de ned in the next section provide a complete picture of the covariance structure of the vector of order statistics, (X1:n ; X2:n ; : : : ; Xn:n ; Y1:n ; Y2:n ; : : : ; Yn:n ). Maritz (1991) has proposed an exact bootstrap estimator for E(Xi:n Yj:n ) based upon an estimate of the joint ˆ y) for F(x; y) in H (x; y) = P(Xr:n ¡ x; Ys:n ¡ y), c.d.f. for the bivariate order statistics given by replacing F(x; as given in David (1981), where Â1 = F(x; y); Â2 = F(x; +∞) − F(x; y); Â3 = F(+∞; y) − F(x; y); Â4 = 1 − F(x; +∞) − F(+∞; y) + F(x; y); H (x; y) =
n X n X u X i=r j=s k=t
n!  k  i−k Â3j−k Â4n−i−j+k ; k!(i − k)!(j − k)!(n − i − j + k)! 1 2
(1.7)
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
197
t = max(0; i + j − n) and u = min(i; j). The estimated probability mass at (x; y) is then given by Hˆ (x+; y+) − Hˆ (x−; y+) − Hˆ (x+; y−) + Hˆ (x−; y−);
(1.8)
ˆ y) for F(x; y) obtained as per the prescription of Maritz (1991) by substituting the empirical estimator F(x; in (1.7), where x+ implies x + and x− implies x − ( small). The approach of Maritz poses an obvious problem if the data are highly correlated such that i equals [i], where the pair (Xi:n Y[i]:n ); i = 1; 2; : : : ; n, denote the bivariate sample data ordered by X , where Y[i]:n denotes the corresponding ith concomitant statistic; see David (1981) for a description of concomitants. In the case i = [i] it is easy to see that Hˆ (x; y) is 0 for all possible (x; y) combinations. Later on we will illustrate that even in the cases where Hˆ (x; y) is nonzero for some (x; y) pairs the estimate of E(Xi:n Yj:n ) based upon (1.8) has some shortcomings. As an alternative, we propose an easy to calculate bootstrap estimator of the joint probability distribution for ˆ r:n ; Ys:n ) = 1=n Xr:n and Ys:n based upon a simple rank transformation of the joint density function given by f(X when s = [r], and 0, elsewhere. The transformation maintains the appropriate marginal weights given by (1.4) such that the standard bootstrap estimators of the marginal expectations and variances remain unchanged while leading to a new estimator of covariance and=or correlation of Xr:n and Ys:n di erent from that based upon (1.8). Details of the calculations are provided in Section 2. In Section 3 we outline a method for constructing the minimum area joint percentile con dence regions for the bivariate quantiles QX (u1 ) and QY (u2 ) and compare it to the asymptotic con dence ellipse approach and a Bonferroni correction-type approach using a biostatistical example speci cally examining bivariate medians. Furthermore, some crude simulations are carried out in order to examine the coverage probabilities of the the minimum area joint percentile con dence region for bivariate medians. Section 4 brie y describes an approach to order statistic regression. 2. Density and moment estimation In this section we derive the estimated joint probability density function for Xr:n and Ys:n , the bootstrap-type l m Ys:n ), as well as a correlation estimator for bivariate order statistics. Towards this end, let estimator of E(Xr:n (xk:n ; y[k]:n ) denote the kth observed concomitant sample pair, k = 1; 2; : : : ; n ordered by x with corresponding sample concomitant y[k]:n . Let k −1 i−1 ; ; rx1 (i; k) = max FUr:n n n k i ; ; rx2 (i; k) = min FUr:n n n [k] − 1 j−1 ; ; ry1 (j; k) = max FUs:n n n [k] j ; ; ry2 (j; k) = min FUs:n n n
k −1 k [k] − 1 [k] i i−1 j j−1 − ; FUs:n − ; − FUr:n ; − FUs:n b(i; j; k) = min FUr:n n n n n n n n n
for all i; j = 1; 2; : : : ; n, where FUr:n (·) denotes the cumulative distribution function of a beta random variable with parameters r − 1 and n − r + 1. Note that FUr:n (·) is readily calculated using most standard software packages.
198
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
The joint empirical density function for the bivariate order statististics Xr:n and Ys:n evaluated at the point (x; y); (x; y) ∈ {(xi−1:n ; xi:n ] × (yj−1:n ; yj:n ]}, is then given by fˆXr:n Ys:n (x; y) = n
n X
{rx2 (i; k) − rx1 (i; k)}{ry2 (j; k) − ry1 (j; k)}1(b(i; j; k)¿0)
(2.1)
k=1
for all i; j = 1; 2; : : : ; n, where 1(b(i; j; k)¿0) denotes the indicator function. The estimator follows straightforward by rst noting that by de nition the classical empirical density funcˆ r:n ; Ys:n ) = 1=n when s = [r], and 0, elsewhere, has uniform mass in the region (Xr−1:n ; Xr:n ] × tion estimator f(X ˆ r:n ; Ys:n ) to fˆX Y (x; y) is then a matter of subdi(Y[r]−1:n ; Y[r]:n ]. Transforming the empirical estimator f(X r:n s:n viding the mass in each region (Xr−1:n ; Xr:n ] × (Y[r]−1:n ; Y[r]:n ] under the constraint that the marginal totals fˆXr:n (x) = wj (r) for x ∈ (Xj−1:n ; Xj:n ] and fˆYs:n (y) = wj (s) for y ∈ (Yj−1:n ; Yj:n ], where wj (·) is de ned at (1.4), i.e, the marginal empirical distribution functions for Xr:n and Ys:n given by FUr:n (Fˆ X (x)) and FUs:n (Fˆ Y (y)), where Fˆ is the classical estimator, de nes a discrete “grid” which intersects the “grid” given by Fˆ X (x) and Fˆ Y (y), each with steps at the order statistics and mass determined by the intersections of the “grids” as given by (2.1). The numerical eciency of our approach relative to Maritz’s approach can be seen by examing (1.7) and (1.8). His analytical bootstrap approach requires 4un2 computations to estimate the point probability as compared to our approach which requires only n computations. Note that our de nition yields a slightly di erent estimator of fˆXr:n Ys:n (x; y) as compared to (1.8), is symmetric in terms of the correlation structure (positive=negative), always produces the appropriate marginal weights and reduces to the standard density estimator when r = s = n = 1. The covariance estimator follows straightforward and is given by the following: l m Ys:n ) is given by Theorem 2.1. The exact bootstrap-type moment estimator of E(Xr:n l m Ys:n )= E ∗ (Xr:n
n X n X i=1 j=1
l m ˆ xi:n yj:n fXr:n Ys:n (xi:n ; yj:n );
(2.2)
where fˆXr:n Ys:n (x; y) is deÿned at (2:1). Proof. First note that given that the sample order statistic X = (X1:n ; X2:n ; : : : ; Xn:n ) is held xed, then a nonparametric bootstrap replication is generated by taking a sample of size n with replacement from X . This is equivalent to generating a random sample of size n from a uniform(0; 1) distribution, with corresponding order statistic U = (U1:n ; U2:n ; : : : ; Un:n ), and applying the sample quantile function −1 ˆ Q(u) = Fˆ (u) = Xbnuc+1:n
(2.3)
ˆ to each element of U , where 0 ¡ u ¡ 1 and b·c denotes the oor function, i.e., Q(u) = Xi:n in the region given by (i − 1)=n6u ¡ i=n; i = 1; 2; : : : ; n. Then let ZZ l m Ys:n )= xl ym fXr:n Ys:n (x; y) d x dy; (2.4) E(Xr:n where fXr:n Ys:n (x; y) denotes the joint density of Xr:n and Ys:n . Rewriting (2.4) in terms of the quantile functions for Xr:n and Ys:n , respectively, yields Z 1Z 1 fXr:n Ys:n {QXr:n (u); QYs:n (v)} l m du dv; Ys:n )= [QXr:n (u)]l [QYs:n (v)]m E(Xr:n fXr:n {QXr:n (u)}fYs:n {QYs:n (v)} 0 0
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
=
n X n Z X i=1 j=1
×
FUr:n (i=n)
FUr:n ((i−1)=n)
Z
FUs:n ( j=n)
FUs:n (( j−1)=n)
199
[QXr:n (u)]l [QYs:n (v)]m
fXr:n Ys:n {QXr:n (u); QYs:n (v)} du dv: fXr:n {QXr:n (u)}fYs:n {QYs:n (v)}
(2.5)
Substituting the empirical estimators Qˆ Xr:n (u) = Qˆ X (QUr:n (u)), Qˆ Ys:n (u) = Qˆ Y (QUs:n (u)), fˆXr:n {QXr:n (u)} = FUr:n (i=n) − FUr:n ((i −1)=n); fˆYs:n {QYs:n (v)}=FYs:n (j=n)−FUs:n ((j −1)=n) and fˆXr:n Ys:n (x; y) at (2.1) into (2:5) and noting that each of the empirical estimators is constant over the regions of integration yields the result at (2.2). The proof follows similarly to what has been done previously by Hutson and Ernst (2000) in the marginal case. Corollary 2.2. The exact bootstrap-type moment estimator of the correlation of Xr:n and Ys:n denoted (Xr:n ; Ys:n ) is given by ∗ (Xr:n ; Ys:n ) =
E ∗ {(Xr:n − E ∗ (Xr:n ))(Ys:n − E ∗ (Ys:n ))} p ; Var ∗ (Xr:n )Var ∗ (Ys:n )
(2.6)
where the expectations E ∗ (Xr:n ) and E ∗ (Ys:n ) are calculated using (1:1) and the variances Var ∗ (Xr:n ) and Var ∗ (Ys:n ) are calculated using (1:2). Proof. Follows straightforward from Theorem 2.1. The following set of simple Examples contrasts are method with the method of Maritz. Examples. Consider for illustration purposes Example 1b given in Maritz (1991). Let n = 5 observations, (xi:n ; y[i]:n ); i = 1; 2; : : : ; 5, be (1; 2); (2; 3); (3; 1); (4; 4); (5; 5): Then from (1.1) and (1.2) E ∗ (X3:5 )=3; E ∗ (Y3:5 )=3; Var ∗ (X3:5 )=0:9824, Var ∗ (Y3:5 )=0:9824. From Theorem 2.1 and Corollary 2.2 we then have E ∗ (X3:5 Y3:5 )=9:55 and ∗ (X3:5 ; Y3:5 )=0:563. Using the estimator Hˆ (x; y) in conjunction with the point probabilities given by (1.8) we get E ∗ (X3:5 Y3:5 ) = 8:64 and ∗ (X3:5 ; Y3:5 ) = −0:369. Note that the Pearson correlation for the above dataset is ˆ = 0:70 indicating that a positive estimate for ∗ (X3:5 ; Y3:5 ) might be more accurate. Now let the n = 5 observations be (1; 1); (2; 2); (3; 3); (4; 4); (5; 5); then we have E ∗ (X3:5 Y3:5 ) = 8:20 and ∗ (X3:5 ; Y3:5 ) = 0:818. In the case of the Maritz estimator the joint empirical distribution function Hˆ (x; y) does not yield a proper joint empirical density estimator. Therefore, the covariance estimator in this case does not make sense. Note that the observations do not necessarily have to be perfectly correlated for this to happen, only that the ranks of the observations turn out to be similar to what happens in this example. In addition, there appears to be some inherent asymmetry in the approach developed by Maritz (1991) with regards to positively and negatively correlated data, i.e. using the estimator based upon (1.8) ∗6= − ∗ corresponding to the case when y is changed to −y. We conjecture that some type of continuity correction applied to Hˆ (x; y) might provide a better estimator of (Xr:n ; Ys:n ), which will yield results that are close to the results obtained by using our estimator of ∗ (Xr:n ; Ys:n ) de ned at (2.6).
200
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
3. Percentile conÿdence regions Hutson (1999) showed that an accurate analytical (1 − ) × 100% bootstrap-type con dence interval for FX−1 (u) = QX (u); 0 ¡ u ¡ 1, is given by (Qˆ X (QUn0 u:n ( =2)); Qˆ X (QUn0 u:n (1 − =2)));
(3.1)
where QUn0 u:n (1 − =2) denotes the quantile function for ÿ(n0 u; n0 (1 − u)) random variate and by de nition Qˆ X (u) = (1 − )Xbn0 uc:n + Xbn0 uc+1:n ;
(3.2)
where = n0 u − bn0 uc; n0 = n + 1 and b·c denotes the oor function. An obvious candidate for the joint con dence region for the bivariate quantiles (QX (u1 ); QY (u2 )) is to apply the Bonferroni correction to in (3.1) and then de ning the joint con dence region to be rectangular. In this section we show that a tighter nonparametric percentile con dence region can be obtained by trimming fˆXr:n Ys:n (x; y) at (2.1) in a speci c way and therefore taking advantage of the correlation structure between Xr:n and Ys:n in the con dence region estimation procedure. For this estimation procedure de ne Qˆ X (u1 ) = xbnu1 c+1:n and Qˆ Y (u2 ) = ybnu2 c+1:n . The algorithm for calculating the nonparametric percentile con dence region is as follows: 1. Eliminate all pairs (xi:n ; yj:n ) over the i, j combinations where fˆXr:n Ys:n (xi:n ; yj:n ) = 0. 2. Calculate the distance (di; j ) between the point estimate ranks (bnu1 c + 1; bnu2 c + 1) to each of the ranks (i; j) for the pairs remaining from step (1). 3. If 0 ¡ fˆXr:n Ys:n (xi:n ; yj:n ) ¡ then calculate the distance-to-probability ratio ri; j = di; j = fˆXr:n Ys:n (xi:n ; yj:n ). 4. Sort the remaining (xi:n ; yj:n ) pairs by ri; j in descending fashion, i.e. largest to smallest. 5. Trim the sorted pairs based upon each successive ri; j one by one and subtract the corresponding probability fˆXr:n Ys:n (xi:n ; yj:n ) from the total probability of 1 contained in the region. 6. Repeat step (5) until the reduction of the probability from 1 is near or equal to (1 − ). 7. The approximate (1 − ) × 100% percentile con dence region is then given by the convex hull generated from the remaining points, such that (Qˆ X (u1 ); Qˆ Y (u2 )) is an interior point of the region. Note that alternative de nitions where considered and may provide smaller con dence regions for speci c cases. The method de ned above was a trade-o between ease-of-use in terms of programming relative to generating smallest possible con dence region. Furthermore, a continuity correction may be applied to the process by using the linear interpolation quantile function estimator at (3.2) in conjunction with fˆXr:n Ys:n (x; y) at (2.1) and setting r = (n + 1)u1 and s = (n + 1)u2 . Another possibility includes appealing to the joint asymptotic normality of the bivariate quantiles (Qˆ X (u1 ); ˆ QY (u2 )) and generating a con dence ellipse based upon the moment estimators de ned in the previous section and given by the following equation: (x − E ∗ )(y − EY∗ ) (y − EY∗ )2 (x − EX∗ )2 − 2∗ p X ∗ + = (1 − ∗2 ) 22; ; ∗ Var X Var ∗Y Var X Var ∗Y
(3.3)
where EX∗ = E ∗ (Qˆ X (u1 )); Var ∗X = Var ∗ (Qˆ X (u1 )); EY∗ = E ∗ (Qˆ Y (u2 )); Var ∗Y = Var ∗ (Qˆ Y (u2 )); ∗ = ∗ (Qˆ X (u1 ); Qˆ Y (u2 )) and 22; denotes the upper th quantile of a chi-square distribution with 2 degrees of freedom. Example. The following example is taken from an echocardiographic imaging (EI) study, Geiser et al. (1998). The data consists of the fractional area change (FAC) from n = 101 subjects computed as AED − AES ×100%; (3.4) AED where AED and AES are the endocardial areas at end diastole and end systole, respectively. The dataset plotted in Fig. 1 contains the FACs as measured by two distinct echocardiographers (Expert 1(X ) and Expert FAC =
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
201
Fig. 1. Three types of 95% con dence regions for bivariate medians.
2(Y )), the three di erent types of con dence regions overlaying the data and a separate blown-up plot of the con dence regions. The dotted line represents the con dence ellipse, the solid line outline the rectangular shape denotes the Bonferroni region and the solid line outlining the convex hull denotes the nonparametric percentile con dence region for (QX (1=2); QY (1=2)). The bootstrap moment estimates for this dataset are as follows: E ∗ (Qˆ X (1=2)) E ∗ (Qˆ Y (1=2)) Var ∗ (Qˆ X (1=2)) Var ∗ (Qˆ Y (1=2)) ∗ (Qˆ X (1=2); Qˆ Y (1=2)) 43.76
45.29
3.65
1.80
0.86
202
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203 Table 1 Simulation results n
0.0
0.5
0.9
25 51 101
92=100 92=100 90=100
91=100 95=100 96=100
84=100 92=100 92=100
We see from Fig. 1 that the con dence ellipse and the nonparametric con dence regions do not di er substantially and provide a much more precise region as compared to the Bonferroni-type approach to the problem in terms of having a region with the minimal area that captures the correlation structure of the data. Preliminary simulation results: We carried out a crude simulation study for bivariate medians for data sampled from a standard bivariate normal distribution with correlation . Due to the complexity of programming such a simulation we resorted to examing the coverage probabilities for the 95% regions for (QX (1=2); QY (1=2)) via inspection of 100 simulations for samples of size n = 25; 51; 101 and = 0:0; 0:5; 0:9. The coverage results are provided in Table 1. From Table 1 we see that for this basic simulation study the coverage appears to be reasonable for samples of size n = 50 or larger. Obviously, further studies need to be carried out before more general statements can be made about the coverage error of the nonparametric region given alternative underlying distributions.
4. Order statistics regression Suppose that Xr:n and Ys:n are linked by the linear regression relation Ys:n = ÿ0 + ÿ1 Xr:n + :
(4.1)
Then it follows from the theory of least squares that nonparametric bootstrap-type estimators of ÿ0 and ÿ1 are given by minimizing n n X X i=1 j=1
(yj:n − (ÿ0 + ÿ1 xi:n ))2 fˆXr:n Ys:n (xi:n ; yj:n );
(4.2)
with respect to ÿ0 and ÿ1 , where fˆXr:n Ys:n (x; y) is de ned at (2.1). Note that the regression of Ys:n on Xr:n has unique solution and goes through the point (E ∗ (Xr:n ); E ∗ (Ys:n )). The solution to (4.2) is easily carried out using standard software packages via weighted least-squares regression. A special case of (4.1) of interest is the regression of the rth concomitant, Y[r]:n , on the rth order statistic Xr:n , e.g. see David (1981). The approach is easily generalized to a multivariate order statistics regression setup by considering (1.3) in conjunction with (2.2).
Acknowledgements We are grateful to the referees for carefully reading this manuscript and o ering helpful suggestions. Alan Hutson’s work was supported in part by NIH General Clinical Research Center Grant RR0082.
A.D. Hutson / Statistics & Probability Letters 48 (2000) 195 – 203
203
References David, H.A., 1981. Order Statistics, 2nd Edition. Wiley, New York. Geiser, E.A., Wilson, D.C., Wang, D.X., Conetta, D.A., Murphy, J.D., Hutson, A.D., 1998. Autonomous epicardial and endocardial boundary detection in echocardiography short-axis images. J. Amer. Soc. Echocardiography 11, 338–348. Hutson, A.D., 1999. Calculating nonparametric con dence intervals for quantiles using fractional order statistics. J. Appl. Statist. 26, 339–349. Hutson, A.D., Ernst, M.D., 2000. The exact bootstrap mean and variance of an L-estimator. J. Roy. Statist. Soc. Ser. B 62, 1–6. Maritz, S., 1991. Estimating the covariance matrix of bivariate medians. Statist. Probab. Lett. 12, 305–309. Small, C.G., 1990. A survey of multidimensional medians. Internat. Statist. Inst. Rev. 58, 262–277.