Derivation of Passing–Bablok regression from Kendall’s tau

2.1 Synopsis of the relevant regression methods

In the following, a caret will designate the estimated value, e.g., of the slope m ^ , and a bar the mean, e.g., x ‾ = ∑ i x i / n . Furthermore, X : = { x i } , and Y : = { y i } .

It will be shown how the regression schemes mentioned in the introduction, namely, LSR, LPR, Deming regression, and notably PBR, may be derived from the condition that either Pearson’s or Kendall’s correlation diminishes to 0 for transformed input data.

The following assumptions are made, although for some methods further generalizations are possible. For LSR and Deming regression, both weighted and unweighted variants exist. For the ease of presentation, only the unweighted variants will be discussed.

1. There exists a structural relationship between the expected values of random variables x _i and y _i , denoted as x i * and y i * , respectively, of the form y * = b + m x * . Here, x i = x i * + ξ i and y i = y i * + η i . For LSR and TSR, ξ i = 0 .

2. Error terms for points i ≠ j are independent. ξ i and η i are independent, too.

3. E ( ξ i ) = E ( η i ) = 0 , for all i.

4. The ratio of the variances is independent of i, γ 2 : = σ η i 2 / σ ξ i 2 = constant .

5. The distribution functions of the error terms may depend on i, but for fixed i, the distributions of ξ i and η i are equal up to scaling with the standard deviation^[2], i.e., ∀ i ,   x : P ( ξ i / σ ξ i < x ) = P ( η i / σ η i < x ) .

For a detailed description of the estimation procedures to be considered, the reader is referred to the literature [1], [2], [3], [4], [6]. Here, only a very condensed synopsis of the main assumptions and formulas for the estimated slope are given for the methods of interest, as follows:

1. LSR: In least squares regression of Y on X, the estimate of the slope m ^ is obtained by minimizing the sum of the squared deviations of the y _i from the estimated values, ∑ i ( y i − y ‾ − m ( x i − x ‾ ) ) 2 as

1. LPR: In least product regression, the estimate m ^ is obtained minimizing the sum of the products of the deviations of the y _i and x _i from their estimated values, ∑ i ( y i − y ‾ − m ^ ( x i − x ‾ ) ) ( x i − x ‾ − ( y i − y ‾ ) / m ^ ) as

where the sign of the root is chosen equal to the sign of Pearson’s correlation coefficient of x and y. Sometimes, LPR is called geometric mean regression, as the estimate m ^ may also be expressed as the geometric mean of the slope estimate from a LSR of y on x and the reciprocal slope estimate of regression of x on y. Also the terminology “reduced major axis regression” is common.

1. Deming Reg.: In major axis, or Deming regression, the slope is determined by minimization of both the sum of the squared deviances of x and y from the estimated true values x ^ i * and y ^ i * = y ‾ + m ^ ( x ^ i * − x ‾ ) ,

where the sum is minimized with respect to both m ^ and the x ^ i * . The slope estimate is

where s x y = ∑ i ( x i − x ‾ ) ( y i − y ‾ ) , s x x = ∑ i ( x i − x ‾ ) 2 , and s y y = ∑ i ( y i − y ‾ ) 2 . It can be seen that the slope estimators of LSR and LPR arise as special cases of the Deming regression estimate. Namely, γ = ∞ corresponds to regression of y on x and γ = 0 to regression of x on y. Assuming that γ = m , replacement of γ by its estimator m ^ in the last expression, yields the LPR expression for m ^ .

1. Roos Reg.: As first discussed by Roos [11] and Xu [12], LSR may be considered as a special case from a family of regressions where the squared distance of the points from a regression line is minimized along lines with slope κ, e.g., LSR is a special case with κ = ∞ . The slope estimate is

1. TSR: According to Sen [8], in Theil–Sen regression the slope can be estimated from the condition that Kendall’s correlation coefficient τ between the X and the intercepts Y − m X shall be 0,

N is a normalization factor. When there are no ties, N = 2 / n ( n − 1 ) . The function sgn (x) is the sign of x, if x ≠ 0 , and 0, if x = 0 . If the x _i are ordered so that x j > x i if j > i , then

where S _ij is the slope of the line connecting points i and j.

1. cPBR: In classical Passing–Bablok regression [6] proposed to modify the TSR estimator, to take account of the X not being free of error. If K is the number of the slopes S i j < − 1 and N = n ( n − 1 ) / 2 the total number of pairs, then

where F _S is the empirical cumulative distribution function of the { S i j } . Furthermore, Passing and Bablok made the assumption that m ≈ 1 .

1. ePBR: In a follow-up article Bablok et al. [10] proposed an alternative method to cPBR which is applicable also for slopes ≠ 1 . They gave a recursive definition which can be shown [13] to be equivalent to

In Figure 1, some of these methods are applied to a data set from clinical chemistry containing paired measurements obtained with both PCR assay and an immunoassay. Due to different units, the numerical ranges of the results differ by several orders of magnitude. It can be seen that the LSR, TSR and also the cPBR are sensitive to the choice of the dependent and independent variables, while the LPR and ePBR yield results which are more centered than those of the other methods.

2.2 Coordinate transformations

In this section it will be shown that all these estimators may be derived from a common expression. Since LSR and LPR can be derived from Deming regression, the latter will be the focus. The idea is to rotate the coordinate system so that either Pearson’s correlation coefficient r, or Kendall’s correlation coefficient τ, vanishes [13]. However, this cannot be done directly because, upon rotation, the error terms will in general be no longer independent:

Let Σ be the covariance matrix of the { ( x i ,   y i ) T } , Σ * of the { ( x i * ,   y i * ) T } and Σ i Ξ that of ξ i ,   η i T , then

where Σ i , 11 = E ( ( x i − x i * ) 2 ) , Σ i , 12 = Σ i , 21 = E ( ( x i − x i * ) ( y i − y i * ) ) , and Σ i , 22 = E ( ( y i − y i * ) 2 ) . Corresponding expressions hold for Σ * , e.g., Σ 11 * = ∑ i ( x i * − x ‾ * ) 2 , while Σ i , 11 Ξ = σ ξ i 2 , Σ i , 12 Ξ = Σ i , 21 Ξ = 0 , and Σ i , 22 Ξ = γ 2 σ ξ i 2 . The rotation angle, which diagonalizes Σ is

which is only consistent if γ 2 = 1 .

However, if new x-values are introduced scaling by γ so that x ′ i = γ x i , the covariance matrix of the new error terms ξ ′ = γ ξ and η will become a multiple of the unit matrix. Upon scaling, the original slope m will be transformed into m ′ = m / γ . The coordinate system can now be rotated by the angle ϕ = − arctan   ( m ′ ) , whence

and

Dropping irrelevant constant factors, X ″ : = { y i + x i γ 2 / m } , and Y ″ : = { y i − m x i } are defined.

The condition that the Pearson correlation r ( X ″ ( m ) ,   Y ″ ( m ) ) of the transformed points diminishes to 0 for m = m ^ becomes

which can be simplified to

It can be shown that the solution of this equation yields exactly the Deming regression estimate m ^ .

When Kendall’s τ is substituted for Pearson’s r in the preceding derivation, differences Δ x i j = x j − x i and Δ y i j has to be considered. The errors of these differences, Δ η i j , and Δ ξ i j will also be independent with the same ratio γ ² of the variances. Hence, the equivalent of Eq. (1) becomes

This can also be simplified to

The solution of this equation yields an estimator^[3] which is the equivalent of Deming regression. For this slope estimate to be consistent, it is necessary that the expectation value of Eq. (4) with m ^ replaced by m must be 0,

If the x i * are not random variables, generically, the expectation values of each term in the sum has to be 0 on its own. This condition is fulfilled if either:

1. γ = | m | and the distributions of the m Δ ξ i j and Δ η i j are identical. The term whose expectation is taken in Eq. (5) is anti-symmetric under an interchange of the m Δ ξ i j and Δ η i j whence the expectation must be 0. These are exactly the assumptions made by Passing and Bablok.

2. Or the m Δ ξ i j − Δ η i j are statistically independent of the corresponding Δ y i j * + γ 2 m Δ x i j * + Δ η i j + γ 2 m Δ ξ i j and their expectations be 0. For general m and γ, this implies that the distribution of the γ ξ i and η i has to be rotationally invariant, which is only true if they are independent samples from a common Gaussian distribution with mean E ( η i ) = E ( ξ i ) = 0 . Although requiring a Gaussian distribution may seem quite restrictive, one should note that the two terms become also approximately statistically independent when the errors are small compared to the mean distance between the points, since it follows that Δ y i j * + γ 2 m Δ x i j * + Δ η i j + γ 2 m Δ ξ i j ≈ Δ y i j * + γ 2 m Δ x i j * in most of the terms, which may therefore be pulled out of the expectation. The expectation E ( sgn ( m Δ ξ i j − Δ η i j ) ) = 0 , if the distribution of the ξ i and η i is symmetric.

The two cases shall be discussed separately in the following subsections. It is remarked that in both cases there are two possible solutions which only differ in the sign of the resulting slope. As in LPR, the sign may be fixed as that of Kendall’s correlation between the original X and Y. In the following descriptions, it will be assumed that this sign is positive, which can always be achieved by eventually changing the sign of the X.

2.3 Passing–Bablok regression

The first case leads to ePBR:

Setting γ 2 = m 2 , finding the m ^ for which the Pearson correlation becomes 0 yields the LPR estimate. Analogously, the condition that Kendall’s correlation shall diminish to 0 leads to an equation for m ^

with

this can also be written as

or,

Unlike the cPBR estimator, a scaling of the S _ij with a factor will also scale the estimator m ^ by the same factor. Therefore, this estimator will be called ePBR estimator in the following.

Alternatively, this estimator arises on introduction of the angles ϕ ^ = arctan   m ^ , and ϕ i j = arctan ( S i j ) , as

Therneau [13] could show that this expression corresponds to an estimator proposed by Bablok et al. [10], who gave a recursive approximation for this estimator. This recursion uses the cPBR as a starting point. It can be seen that if in the first bracket of Eq. (6), m ^ is replaced by 1, this will yield exactly the cPBR estimator. If Pearson correlation is used instead of Kendall’s, the equivalent of the cPBR corresponds to the Roos regression estimate with parameter κ = − 1 .

While Eq. (6) yields a consistent estimator, the bias of the ePBR estimator for small sample sizes remains to be determined. It is assumed that the true m ≠ 0 . Multiplying

with ( m / m ^ ) 2 yields

As m Δ x i j and Δ y i j are distributed alike, it can be concluded that m ^ / m and m / m ^ also follow the same distribution. Therefore, ln ( m ^ ) is symmetrically distributed around ln ( m ) . As the variance of m ^ is O ( 1 / n ) (see below), ln ( m ^ ) is an unbiased estimator of ln ( m ) , and m ^ is an asymptotically unbiased estimate of m. An equivalent result holds for the LPR estimator [14].

Considering the robustness properties of the ePBR estimator, it is easy to show that it has the same minimal breakdown point as the TSR estimator. For the TSR estimator the breakdown point is reached earliest, when the slopes between any two points of which at least one is an outlier make up 50% of all slopes and are either all larger or all smaller than all the slopes between any two points which are not outliers ([15], p. 67). Asymptotically, this is the case when the fraction of the outliers reaches 1 − 1 / 2 ≈ 29 % . The calculation remains true, if, in case of the ePBR, the slopes are replaced by their absolute value, whence the breakdown point of the ePBR is also 29 %.

2.4 Circular median of the doubled angle

Considering the general solution of Eq. (4) for fixed and known γ, setting ϕ : = arctan   ( m / γ ) , ϕ i j = arctan   ( Δ y i j / ( γ Δ x i j ) ) , and making use of addition theorems for the trigonometric functions, one obtains

Here, ⌊ x ⌋ is the floor function, i.e., the largest integer less than or equal to x. ϕ ^ is a standard estimator of circular or, more generally, of directional statistics [16], [17] known as the “circular median” of the doubled angle. Doubling of the angle is standard in circular statistics for data showing diametrically bimodal distributions. The circular median has been shown to be a robust and unbiased location estimator of the slope angle [18]. To visualize the circular median, the (doubled) angles are represented as points on a unit circle and a diameter PQ is drawn which separates the circle into two halves containing the same number of points, cf. Figure 2. The point P, for which the majority of points have less distance than to Q, is then the circular median. For an odd number of points it coincides with a measured angle, while for an even number of points the mean of the nearest two points is used.

Figure 2:

A: Histogram of all unique pairwise slope angles for the example from Figure 1 after scaling with the PB slope estimate. Light grey: slopes which contribute +1 to Kendall’s correlation function; Dark grey: slopes which contribute −1 to the correlation function. B: Histogram for the doubled slope angles (modulo 2π). An equal amount of points lie left (light grey) and right (dark grey) of the line passing through the circular median at π / 2 and the origin of the circle.

Estimation of the slope via the circular median seems to be an attractive and robust alternative to the usual Deming regression, especially if the ratio of the standard errors is known and different from the expected slope. Similar to the way the Deming slope estimator interpolates between LSR estimator for the slope of Y versus X and the inverse of X versus Y depending on the ratio γ, the circular median can be seen to interpolate between the corresponding TSR estimates. The calculation of the PBR estimate is also possible using the circular median algorithm.

While finding the circular median is as complex as finding the usual median as required in the practical computation of the PBR, the concept offers the advantage, that it can be extended to find a regression line when comparing more than two sets of data [19].

It is of interest to compare the ePBR and the circular median method numerically. To that end, 100 points were sampled 1000 times with x * = y * from a uniform distribution, i.e., with m = 1 , and a) equal coefficient of variance c v y = c v x = 0.3 and b) unequal coefficient of variance c v x = 0.3 and c v y = 0.9 , with Gaussian error. In case a), both methods yield practically the same results with mean m ^ = 1.005   ( sd = 0.067 ) for the ePBR and mean m ^ = 1.007   ( sd = 0.084 ) for the median circular slope. However, in case b), the ePBR is biased with mean m ^ = 1.606   ( sd = 0.160 ) , while the mean median circular slope estimate is m ^ = 1.007   ( sd = 0.084 ) .

2.5 Numerical complexity of the calculation of the PBR and circular median estimators

The question remains how efficiently the TSR, PBR, and circular median estimate can be calculated. In case of the TSR, this problem has long been solved and represents a classical problem of computational geometry [20]. While the brute force approach to find the median of all the pairwise slopes between n points has a numerical complexity of O ( n 2 ) , better algorithms have only a complexity of O ( n ln n ) . The problem can be shown to be equivalent to that of computing Kendall’s tau, for which an efficient solution was given by Knight in 1966 [21]. It is based on the observation that the evaluation of Kendall’s τ ( X ″ ,   Y ″ ) is equivalent to that of counting the inversions needed to bring the set Y ″ into the order of X ″ , an operation which can be accomplished in O ( n ln n ) steps using a sorting algorithm like mergesort [22]. x ″ i and y ″ i can be interpreted as intercepts b i ( m ) of lines passing through the points ( x i , y i ) whose slope m is taken as an independent variable. This lends to an interpretation in terms of a duality transformation which maps points ( x i ,   y i ) to lines (parametrized by ( m ,   b i ( m ) ) and lines to points. In the field of image analysis, this transformation is known as Hough transformation [23].

Since Knight’s algorithm has been implemented in many statistical packages and these also take care of the handling of ties, it is an easy task to use these algorithms to calculate the PBR and circular mean estimators. Even in combination with a simple bisection algorithm for root finding it thus becomes possible to find the estimators for data sets with n = 10 7 on an ordinary laptop in a few minutes.