rregressiongroupingextractcoefficients

Extract passing-bablok slopes and intercepts by group using mcreg (mcr package) in R


For about a hundred of different analytes (groups), I would like to compare serum and plasma values using passing-bablok regressions with the mcr::mcreg function, and extract slopes and intercepts (and their confidence intervals) in a dataframe.
I tried to get inspiration from the models proposed models for linear regression, but without success.
These slopes and intercepts are located in the class S4 MCResultBCa object generated by mcreg, which includes the para slot (the numeric matrix with estimates for slope and intercept, and corresponding SD and CI).
Below is an example including three groups of analytes.

Thank you for your help and time.

NB1: please note that I would like to set the mcreg function as follows:

mcreg(x = dat0$serum,
      y = dat0$plasma,
      method.reg = "PaBa",
      method.ci = "bootstrap",
      method.bootstrap.ci = "BCa",
      na.rm = TRUE)  

NB2: considering a S4 MCResultBCa object generated by the mcreg function for e.g. the bm1 analyte, and named e.g. mcreg1, the desired output are given by:

Intercept <- mcreg1@para[1]  
Intercept.LCI <- mcreg1@para[5]  
Intercept.UCI <- mcreg1@para[7]  
Slope <- mcreg1@para[2]  
Slope.LCI <- mcreg1@para[6]  
Slope.UCI <- mcreg1@para[8]  

Desired output:

  analyte nb Intercept Intercept.LCI Intercept.UCI Slope Slope.LCI Slope.UCI
1     bm1 57    -0.451        -2.098         1.262 0.993     0.960     1.023
2     bm2 57    10.651        -9.324        41.424 0.840     0.779     0.898
3     bm3 57    -1.553        -9.174         6.458 1.898     1.389     2.380

Data:

dat0 <- 
structure(list(patient = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 
3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 7L, 7L, 7L, 8L, 8L, 8L, 9L, 
9L, 9L, 10L, 10L, 10L, 12L, 12L, 12L, 13L, 13L, 13L, 14L, 14L, 
14L, 15L, 15L, 15L, 16L, 16L, 16L, 17L, 17L, 17L, 18L, 18L, 18L, 
19L, 19L, 19L, 20L, 20L, 20L, 21L, 21L, 21L, 22L, 22L, 22L, 23L, 
23L, 23L, 28L, 28L, 28L, 29L, 29L, 29L, 30L, 30L, 30L, 31L, 31L, 
31L, 32L, 32L, 32L, 33L, 33L, 33L, 34L, 34L, 34L, 35L, 35L, 35L, 
36L, 36L, 36L, 37L, 37L, 37L, 38L, 38L, 38L, 39L, 39L, 39L, 40L, 
40L, 40L, 41L, 41L, 41L, 42L, 42L, 42L, 43L, 43L, 43L, 44L, 44L, 
44L, 45L, 45L, 45L, 46L, 46L, 46L, 47L, 47L, 47L, 48L, 48L, 48L, 
49L, 49L, 49L, 50L, 50L, 50L, 51L, 51L, 51L, 52L, 52L, 52L, 53L, 
53L, 53L, 55L, 55L, 55L, 56L, 56L, 56L, 57L, 57L, 57L, 58L, 58L, 
58L, 59L, 59L, 59L, 61L, 61L, 61L, 62L, 62L, 62L, 63L, 63L, 63L, 
64L, 64L, 64L, 65L, 65L, 65L), levels = c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
"38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", 
"49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", 
"60", "61", "62", "63", "64", "65", "66"), class = "factor"), 
    analyte = c("bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", "bm2", "bm3", "bm1", 
    "bm2", "bm3"), serum = c(61.6, 904, 10.8, 14.8, 1448, 20.4, 
    30.7, 1396, 12.2, 123, 969, 13.6, 10.9, 1515, 42.9, 62.6, 
    137, 32.2, 74.1, 229, 16.8, 75, 132, 8.56, 69.1, 1974, 16.7, 
    29.2, 972, 16.8, 83.1, 198, 15.7, 30.1, 538, 20, 34, 375, 
    14, 61.7, 768, 29.2, 24.3, 1082, 21, 63.1, 667, 14.4, 33.6, 
    749, 40.1, 36.7, 1118, 17.8, 51.7, 994, 17.1, 60.9, 888, 
    45, 37.5, 1381, 15.1, 45.5, 1229, 9.26, 127, 289, 11.7, 25.2, 
    928, 21.5, 88.6, 1536, 10.3, 97.7, 717, 20.1, 11.2, 637, 
    51.8, 102, 271, 9.15, 70.6, 103, 13, 40.6, 811, 6.53, 120, 
    1374, 24.6, 81, 318, 20.3, 49.5, 366, 26.1, 51.2, 696, 25.9, 
    116, 5063, 31.3, 87.1, 804, 14, 60, 906, 12, 91.3, 597, 11.3, 
    27.9, 831, 27.6, 45.2, 747, 19, 39.3, 179, 19.3, 29.1, 504, 
    18.8, 33.5, 795, 14, 18.9, 1053, 6.46, 13.4, 432, 21.7, 26.3, 
    462, 14.4, 88.2, 937, 30.8, 81.9, 6000.001, 12.3, 45.7, 36.5, 
    16.4, 34.6, 544, 9.27, 43.1, 2814, 19.5, 69.5, 1084, 19.5, 
    72.4, 459, 23.9, 137, 162, 13.2, 29.4, 1910, 24.8, 99.4, 
    517, 10.4, 46.9, 1334, 14.8), plasma = c(57.8, 792, 27.8, 
    17.4, 1436, 26.6, 32.5, 1075, 27.1, 118, 794, 34.4, 10.7, 
    901, 54.4, 70.3, 165, 12.4, 75.2, 193, 34, 74, 116, 15.1, 
    69.5, 1425, 20, 27.5, 804, 23.3, 81, 169, 45.1, 29.3, 413, 
    29.8, 28.9, 326, 22.3, 61, 648, 42.9, 20, 635, 30.2, 65.9, 
    640, 37.2, 32.6, 633, 88.9, 36, 900, 47.3, 46.9, 899, 40.9, 
    55.2, 778, 80.7, 38.5, 1069, 23.1, 42.7, 1167, 22, 119, 261, 
    37, 22.8, 953, 28.2, 88.2, 1264, 31.7, 102, 638, 36.6, 12.7, 
    541, 41.8, 102, 231, 13.7, 68.8, 79.6, 12.5, 41, 770, 10.5, 
    115, 1105, 37.2, 80.3, 292, 41.6, 51.7, 318, 22.7, 49.1, 
    671, 43.9, 117, 3312, 42.5, 85.7, 683, 24.7, 55.3, 792, 24.6, 
    85.6, 532, 38.9, 26.4, 662, 28.7, 43.1, 662, 25.1, 40.9, 
    138, 48.1, 27.3, 471, 21.8, 32.8, 685, 22.9, 18.7, 967, 11.9, 
    11.5, 400, 30.8, 23.7, 427, 27.7, 88.6, 739, 59.6, 76.8, 
    5397, 33.9, 47.8, 24.2, 39.7, 32.6, 516, 22.1, 42.5, 2584, 
    32.5, 70.3, 1087, 36.2, 73, 405, 45.3, 127, 146, 49, 31.8, 
    1251, 45.2, 107, 402, 25.2, 47.3, 1142, 35)), row.names = c(NA, 
-171L), class = "data.frame")  

Solution

  • The following code works well, even if the run is a bit long (especially for an iteration on a hundred groups).
    If you find a way to make a shorter code, with a faster run, and if possible avoiding the use of plyr::ddply, that would be great!

    library(dplyr)
    library(plyr)
    library(mcr)
    
    set.seed(123)
    
    for (analyte in unique(dat0$analyte)){
      dat1 <- subset(dat0, dat0$analyte == analyte)
      
    # paba model
    dat1_paba_mod <-
      ddply(dat1, .(analyte), transform,
            fit =
              calcResponse(
                mcreg(
                  x = serum,
                  y = plasma,
                  method.reg = "PaBa",
                  method.ci = "bootstrap",
                  method.bootstrap.ci = "BCa",
                  na.rm = TRUE),
                x.levels = serum)
            )
    
    # paba regression
    dat1_paba_reg <- function(dat) {
      mcreg <-
        mcreg(
          x = dat$serum,
          y = dat$plasma,
          method.reg = "PaBa",
          method.ci = "bootstrap",
          method.bootstrap.ci = "BCa",
          na.rm = TRUE
          )
      nb <- nrow(dat)
      intercept <- round(mcreg@para[1], 3)
      intercept_lci <- round(mcreg@para[5], 3)
      intercept_uci <- round(mcreg@para[7], 3)
      slope <- round(mcreg@para[2], 3)
      slope_lci <- round(mcreg@para[6], 3)
      slope_uci <- round(mcreg@para[8], 3)
      data.frame(
        nb = nb,
        Intercept = intercept,
        Intercept.LCI = intercept_lci,
        Intercept.UCI = intercept_uci,
        Slope = slope,
        Slope.LCI = slope_lci,
        Slope.UCI = slope_uci,
        stringsAsFactors = FALSE
        )
      } 
    
    # do paba_reg by group
      dat1_paba_gp <-
        dat1 %>%
        group_by(analyte) %>%
        do(dat1_paba_reg(.))
      dat1_paba_gp
    } 
    
    dat1_paba_gp
    
    # > dat1_paba_gp
    # # A tibble: 3 × 8
    # # Groups:   analyte [3]
    # analyte    nb Intercept Intercept.LCI Intercept.UCI Slope Slope.LCI Slope.UCI
    # <chr>   <int>     <dbl>         <dbl>         <dbl> <dbl>     <dbl>     <dbl>
    # 1 bm1        57    -0.451         -2.10          1.25 0.993     0.956     1.02 
    # 2 bm2        57    10.7           -9.41         48.9  0.84      0.771     0.899
    # 3 bm3        57    -1.55         -10.6           5.94 1.90      1.39      2.46