Unobvious problems of using the R's implementation of the Hodges-Lehmann estimator

The Hodges-Lehmann location estimator (also known as pseudo-median) is a robust, non-parametric statistic used as a measure of the central tendency. For a sample $\mathbf{x} = \{ x_1, x_2, \ldots, x_n \}$, it is defined as follows:

Essentially, it’s the median of the Walsh (pairwise) averages.

For two samples $\mathbf{x} = \{ x_1, x_2, \ldots, x_n \}$ and $\mathbf{y} = \{ y_1, y_2, \ldots, y_m \}$, we can also consider the Hodges-Lehmann location shift estimator:

In R, both estimators are available via the wilcox.test function. Here is a usage example:

set.seed(1729)
x <- rnorm(2000, 5) # A sample of size 2000 from the normal distribution N(5, 1)
y <- rnorm(2000, 2) # A sample of size 2000 from the normal distribution N(2, 1)
wilcox.test(x, conf.int = TRUE)$estimate
# (pseudo)median
#       5.000984
wilcox.test(y, conf.int = TRUE)$estimate
# (pseudo)median
#       1.969096
wilcox.test(x, y, conf.int = TRUE)$estimate
# difference in location
#               3.031782

In most cases, this function works fine. However, there is an unobvious corner case, in which it returns wrong values. In this post, we discuss the underlying problem and provide a correct implementation for the Hodges-Lehmann estimators.

Problem 1: Zero values

Let us consider the one-sample estimator for a sample that contains exactly one zero element:

x <- c(0, 1, 2)
wilcox.test(x, conf.int = TRUE)$estimate
# (pseudo)median 
#            1.5 
# Warning messages:
# 1: In wilcox.test.default(x, conf.int = TRUE) :
#   requested conf.level not achievable
# 2: In wilcox.test.default(x, conf.int = TRUE) :
#   cannot compute exact p-value with zeroes
# 3: In wilcox.test.default(x, conf.int = TRUE) :
#   cannot compute exact confidence interval with zeroes

Obviously, for $\mathbb{x} = \{ 0, 1, 2 \}$, the correct pseudo-median value is $1$. However, wilcox.test returns $1.5$ (a wrong value) and prints a set of warnings about a zero value in the sample. The problem can be found in the source of wilcox.test:

ZEROES <- any(x == 0)
if(ZEROES)
    x <- x[x != 0]

For the one-sample case, wilcox.test removes all zero elements from the sample. This logic is needed to properly perform the Wilcoxon rank-sum test (also known as the Mann–Whitney U test). The Hodges-Lehmann estimation is an additional feature of this function. Unfortunately, this feature is affected by this cleanup of zero values. Therefore, it actually estimated the pseudo-median not for $\mathbf{x} = \{ 0, 1, 2 \}$, but for $\mathbf{x}' = \{ 1, 2 \}$ (which is $1.5$).

Problem 2: Tied values (one sample)

Now let us consider the following sample of four elements with a tie:

The true value of the Hodges-Lehmann location estimator is given by:

Now let us look at the calculated result by wilcox.test:

x <- c(-2.12984, -2.12984, 1.1479, -0.4895)
wilcox.test(x, conf.int = TRUE)$estimate
# (pseudo)median 
#     -0.6514901 
# Warning messages:
# 1: In wilcox.test.default(x, conf.int = TRUE) :
#   requested conf.level not achievable
# 2: In wilcox.test.default(x, conf.int = TRUE) :
#   cannot compute exact p-value with ties
# 3: In wilcox.test.default(x, conf.int = TRUE) :
#   cannot compute exact confidence interval with ties

As we can see, the returned value of $-0.6514901$ significantly differs from the expected value of $-0.90032$. When the sample contains tied values, wilcox.test switched from the exact implementation of Wilcoxon rank-sum test to the approximated one. As a side effect, it also uses a peculiar approximation of the Hodges-Lehmann estimator that leads to another pseudo-median estimation that differs from the explicit equation.

Problem 3: Tied values (two samples)

When we estimate the location shift between two samples, tied values are also an issue. Let us consider the following samples:

The expected value of $\operatorname{HL}(\mathbf{x}, \mathbf{y})$ is $\approx 0.1285858$. Now let us check the output of wilcox.test:

x <- c(1.5274454801712, 1.5274454801712, 0.3)
y <- c(3.3, -1.72972619537396)
wilcox.test(x, y, conf.int = TRUE)$estimate
# difference in location
#               1.503729

We have $\approx 0.1285858$ vs. $\approx 1.503729$. It is a huge difference! The underlying problem is the same as the previous one: existing of tied values forces wilcox.test to switch to the approximated algorithm that returns strange results.

Problem 4: Degenerate samples

When sample ranges are degenerate ($\min(x) = \max(x)$, $\min(y) = \max(y)$), we get a corner case that is not supported in wilcox.test:

x <- c(2, 2)
y <- c(1, 1)
wilcox.test(x, y, conf.int = TRUE)$estimate
# Error in if (f.lower <= 0) return(mumin) : 
#   missing value where TRUE/FALSE needed
# In addition: Warning messages:
# 1: cannot compute confidence interval when all observations are tied 
# 2: cannot compute confidence interval when all observations are tied 

While the actual value of $\operatorname{HL}(\mathbf{x}, \mathbf{y})$ is obviously $1$, wilcox.test fails to provide any result.

The correct implementation

Given the problems discussed above, I recommend against using wilcox.test for Hodges-Lehmann estimations due to its numerous corner cases Instead, I suggest the following simple and reliable implementation that supports all the discussed scenarios:

hl <- function(x, y = NULL) {
  if (is.null(y)) {
    walsh <- outer(x, x, "+") / 2
    median(walsh[lower.tri(walsh, diag = TRUE)])
  } else {
    median(outer(x, y, "-"))
  }
}

Let us check that it works correctly:

x <- c(0, 1, 2)
hl(x)
# [1] 1

x <- c(-2.12984, -2.12984, 1.1479, -0.4895)
hl(x)
# [1] -0.90032

x <- c(1.5274454801712, 1.5274454801712, 0.3)
y <- c(3.3, -1.72972619537396)
hl(x, y)
# [1] 0.1285858

x <- c(2, 2)
y <- c(1, 1)
hl(x, y)
# [1] 1