Block-Passing-Bablok Regression

R implementation of Block-Passing-Bablok (Block-PBR) regression for method comparison studies with repeated measurements, as described in:

F. Baumdicker, U. Hölker (2020). Method comparison with repeated measurements — Passing–Bablok regression for grouped data with errors in both variables. Statistics & Probability Letters, 164, 108801. https://doi.org/10.1016/j.spl.2020.108801

Freely accessible preprint: https://arxiv.org/abs/1905.07649

Background

Passing-Bablok regression (PBR) is a non-parametric method widely used in clinical biochemistry and laboratory medicine to compare two measurement methods. It estimates the slope β and intercept α of the linear relationship y = α + βx by taking the median of all pairwise slopes between measurements.

When the data contain repeated measurements (multiple observations from the same sample or group), the within-group slopes have an expected value of 0 rather than β, and including them biases the slope estimate toward zero. Block-PBR corrects this by excluding within-group slopes from the median calculation and provides a corresponding asymptotic confidence interval that accounts for group sizes and inter-group overlap.

Files

File	Description
mcPaBa_custom.r	Modified mc.paba() function from the mcr R package (originally Roche Diagnostics GmbH, GPL-3), extended with a groupsizes parameter to implement the block modification and corrected confidence intervals.
useBlockPaBa_example.R	Example script that simulates grouped data and demonstrates both classical PBR and Block-PBR, including a PDF plot of the results.

Requirements

Install the mcr package from CRAN:

install.packages("mcr")

Both scripts must be in the same working directory.

Usage

The example script useBlockPaBa_example.R walks through the full workflow:

1. Simulate (or load) data — the script simulates three groups of sizes 40, 180, and 20 with a true slope of 0.8. Replace the simulated data with your own measurements.
2. Run classical PBR via mcreg() from the mcr package.
3. Compute the angle matrix with mc.calcAngleMat(), then set within-group entries to NA to exclude them.
4. Run Block-PBR by passing the masked angle matrix to mc.paba() together with the group sizes.

library(mcr)
source("mcPaBa_custom.r")
calcDiff <- function(x, y) { x - y }

# --- replace with your data ---
# xr: measurements from reference method
# yr: measurements from test method
# groups: integer vector assigning each observation to a group

groupnames <- unique(groups)

# compute angle matrix and mask within-group slopes
myAngleMat <- mc.calcAngleMat(xr, yr)
for (name in groupnames) {
  idx <- which(groups == name)
  myAngleMat[as.matrix(expand.grid(idx, idx))] <- NA
}

groupsizes <- as.vector(table(groups))

# fit Block-Passing-Bablok
classicresult <- mcreg(xr, yr, method.reg = "PaBa", method.ci = "analytical",
                       slope.measure = "radian")
blockPBresult <- classicresult
blockPBresult@para <- mc.paba(myAngleMat, xr, yr,
                              groupsizes = groupsizes,
                              slope.measure = "radian")

print(blockPBresult@para)

The estimated slope, intercept, and their 95% confidence intervals are in blockPBresult@para. Expected numerical output and a reference PDF plot for the simulated data are provided in the expected_output/ folder for comparison.

Advanced: overlapping groups

By default, mc.paba() assumes that the groups do not overlap on the x-axis, which simplifies the confidence interval formula (Remark 3.8 in the paper). If measurements from different groups share a common range of x-values, this assumption is violated and the overlap between groups should be accounted for in the CI calculation.

How severe is the omission? As shown in Remark 3.8 of the paper, omitting the overlap correction overestimates the variance of the estimator, which leads to CIs that are wider than necessary. For the equivalence test β = 1 this means the test is conservative: CI coverage meets or exceeds the nominal level, but power to detect a true deviation from β = 1 is reduced. The simulations in Table 1 and Figure 1 of the paper show that in practice this power loss is negligible when β is close to 1 — the typical null hypothesis in method comparison studies — even under high overlap. However, when β is far from 1 (e.g. β = 0.2) and overlap is high, the power reduction becomes substantial and providing an overlap estimate is advisable.

The example script includes a helper function compute_overlapmatrix() that computes the empirical overlap matrix q_{ij} described in the paper using the empirical CDF. To use it, pass the result as the overlap argument to mc.paba():

overlap <- compute_overlapmatrix(groups, groupnames, xr)

blockPBresult@para <- mc.paba(myAngleMat, xr, yr,
                              groupsizes = groupsizes,
                              overlap = overlap,
                              slope.measure = "radian")

Note that the empirical CDF estimate of q_{ij} may be imprecise when the number of repeated measurements per group is small. In such cases, users may want to consider their own overlap estimates, for example based on a parametric model for the distribution of x-measurements within each group. See Section 3 of the paper for the theoretical definition of q_{ij}.

Citation

If you use this code or the method, please cite the paper:

F. Baumdicker, U. Hölker (2020).
Method comparison with repeated measurements — Passing–Bablok regression
for grouped data with errors in both variables.
Statistics & Probability Letters, 164, 108801.
https://doi.org/10.1016/j.spl.2020.108801

A freely accessible preprint is available on arXiv: https://arxiv.org/abs/1905.07649

License

The modifications in mcPaBa_custom.r are by Franz Baumdicker. The underlying mc.paba() function is copyright (C) 2011 Roche Diagnostics GmbH and distributed under the GNU General Public License v3 or later. This code is distributed under the same license.