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")
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