Add a Passing-Bablok regression line

Question

I have to perform many comparisons between different measurement methods and I have to use the Passing-Bablok regression approach.

I would like to take advantage of ggplot2 and faceting, but I don't know how to add a geom_smooth layer based on the Passing-Bablok regression. I was thinking about something like: https://stackoverflow.com/a/59173260/2096356

Furthermore, I would also need to show the regression line equation, with confidence interval for intercept and slope parameters, in each plot.

Edit with partial solution

I've found a partial solution combining the code provided in this post and in this answer.

## Regression algorithm
 passing_bablok.fit <- function(x, y) {
  
    x_name <- deparse(substitute(x))
  
    lx <- length(x)
    l <- lx*(lx - 1)/2
    k <- 0
    S <- rep(NA, lx)
    for (i in 1:(lx - 1)) {
      for (j in (i + 1):lx) {
        k <- k + 1
        S[k] <- (y[i] - y[j])/(x[i] - x[j])
      }
    }
    S.sort <- sort(S)
    N <- length(S.sort)
    neg <- length(subset(S.sort,S.sort < 0))
    K <- floor(neg/2)
    if (N %% 2 == 1) {
      b <- S.sort[(N+1)/2+K]
    } else {
      b <- sqrt(S.sort[N / 2 + K]*S.sort[N / 2 + K + 1])
    }
    a <- median(y - b * x)
    res <- as.vector(c(a,b))

  names(res) <- c("(Intercept)", x_name)
  class(res) <- "Passing_Bablok"
  res
}

## Computing confidence intervals
 passing_bablok <- function(formula, data, R = 100, weights = NULL){
  ret <- boot::boot(
    data = model.frame(formula, data), 
    statistic = function(data, ind) {
      data <- data[ind, ]
      args <- rlang::parse_exprs(colnames(data))
      names(args) <- c("y", "x")
      rlang::eval_tidy(rlang::expr(passing_bablok.fit(!!!args)), data, env = rlang::current_env())
    },
    R=R
  )
  class(ret) <- c("Passing_Bablok", class(ret))
  ret  
}

## Plotting confidence bands
 predictdf.Passing_Bablok <- function(model, xseq, se, level) {
  pred <- as.vector(tcrossprod(model$t0, cbind(1, xseq)))
  if(se) {
    preds <- tcrossprod(model$t, cbind(1, xseq))
    data.frame(
      x = xseq,
      y = pred,
      ymin = apply(preds, 2, function(x) quantile(x, probs = (1-level)/2)),
      ymax = apply(preds, 2, function(x) quantile(x, probs = 1-((1-level)/2)))
    )
  } else {
    return(data.frame(x = xseq, y = pred))
  }
}

An example of usage:

z <- data.frame(x = rnorm(100, mean = 100, sd = 5), 
                y = rnorm(100, mean = 110, sd = 8))

ggplot(z, aes(x, y)) + 
 geom_point() + 
 geom_smooth(method = passing_bablok) + 
 geom_abline(slope = 1, intercept = 0)

So far, I haven't been able to show the regression line equation, with confidence interval for intercept and slope parameters (as +- or in parentheses).

Hello! This question comes across to me like it's asking someone to just do all your work for you. You've been on this site a while now, so I think you already know that you are more likely to receive a useful answer if you provide a reproducible example and show what you have tried so far. — Axeman, Apr 20 '21 at 16:07
Right, with some delay I've added a partial solution to my question. — Ndr, May 21 '21 at 09:05

score 0 · Answer 1 · answered Oct 06 '21 at 22:16

You've arguably done with difficult part with the PaBa regression.

Here's a basic solution using your passing_bablok.fit function:

z <- data.frame(x = 101:200+rnorm(100,sd=10),
                y = 101:200+rnorm(100,sd=8))

mycoefs <- as.numeric(passing_bablok.fit(x = z$x, y=z$y))

paba_eqn <-  function(thecoefs) {
  l <- list(m = format(thecoefs[2], digits = 2),
            b = format(abs(thecoefs[1]), digits = 2))
  if(thecoefs[1] >= 0){
    eq <- substitute(italic(y) == m %.% italic(x) + b,l)
  } else {
    eq <- substitute(italic(y) == m %.% italic(x) - b,l)    
  }
  as.character(as.expression(eq))
}

library(ggplot2)
ggplot(z, aes(x, y)) + 
  geom_point() + 
  geom_smooth(method = passing_bablok) + 
  geom_abline(slope = 1, intercept = 0) + 
  annotate("text",x = 110, y = 220, label = paba_eqn(mycoefs), parse = TRUE)

Note the equation will vary because of rnorm in the data creation..

The solution could definitely be made more slick and robust, but it works for both positive and negative intercepts.

Equation concept sourced from: https://stackoverflow.com/a/13451587/2651663

Add a Passing-Bablok regression line

Edit with partial solution

1 Answers1