Passing-Bablok slope and intercept estimator.

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.

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)

Author

Jakob Raymaekers

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
#>    user  system elapsed 
#>   0.002   0.000   0.002 

t0 <- proc.time()
PB.naive <- bruteForcePB(x, y)
t1 <- proc.time()
t1 - t0
#>    user  system elapsed 
#>   0.076   0.001   0.077 

PB.fast$slope - PB.naive$slope
#> [1] 0