PassingBablok {robslopes} | R Documentation |
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 |
verbose |
Whether or not to print out the progress of the algorithm. Defaults to |
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