PassingBablok {robslopes}

Passing-Bablok slope and intercept estimator.

Description

Computes the equivariant Passing-Bablok regression. The implemented algorithm was proposed by Raymaekers and Dufey (2022) and runs in an expected O(n log n) time while requiring O(n) storage.

Usage

PassingBablok(x, y, alpha = NULL, verbose = TRUE)

Arguments

x

A vector of predictor values.

y

A vector of response values.

alpha

Determines the order statistic of the target slope, which is equal to [alpha*n*(n-1)], where n denotes the sample size. Defaults to NULL, which corresponds with the (upper) median.

verbose

Whether or not to print out the progress of the algorithm. Defaults to TRUE.

Details

Given two input vectors x and y of length n, the equivariant Passing-Bablok estimator is computed as med_{ij} |(y_i - y_j)/(x_i-x_j)|. By default, the median in this experssion is the upper median, defined as \lfloor (n +2) / 2 \rfloor. By changing alpha, other order statistics of the slopes can be computed.

Value

A list with elements:

intecept

The estimate of the intercept.

slope

The Theil-Sen estimate of the slope.

Author(s)

Jakob Raymaekers

References

Passing, H., Bablok, W. (1983). A new biometrical procedure for testing the equality of measurements from two different analytical methods. Application of linear regression procedures for method comparison studies in clinical chemistry, Part I, Journal of clinical chemistry and clinical biochemistry, 21,709-720.

Bablok, W., Passing, H., Bender, R., Schneider, B. (1988). A general regression procedure for method transformation. Application of linear regression procedures for method comparison studies in clinical chemistry, Part III. Journal of clinical chemistry and clinical biochemistry, 26,783-790.

Raymaekers J., Dufey F. (2022). Equivariant Passing-Bablok regression in quasilinear time. (link to open access pdf)

Examples

# We compare the implemented algorithm against a naive brute-force approach.

bruteForcePB <- function(x, y) {
  
  n <- length(x)
  medind1 <- floor(((n * (n - 1)) / 2 + 2) / 2) # upper median
  medind2 <- floor((n + 2) / 2)
  temp <-  t(sapply(1:n, function(z)  apply(cbind(x, y), 1 ,
                                            function(k) (k[2] - y[z]) /
                                              (k[1] - x[z]))))
  PBslope <- sort(abs(as.vector(temp[lower.tri(temp)])))[medind1]
  PBintercept <- sort(y - x * PBslope)[medind2]
  return(list(intercept = PBintercept, slope = PBslope))
}


n = 100
set.seed(2)
x = rnorm(n)
y = x + rnorm(n)

t0 <- proc.time()
PB.fast <- PassingBablok(x, y, NULL, FALSE)
t1 <- proc.time()
t1 - t0

t0 <- proc.time()
PB.naive <- bruteForcePB(x, y)
t1 <- proc.time()
t1 - t0

PB.fast$slope - PB.naive$slope