Bland-Altman plots

Question

I need to make Bland-Altman plots, but I cannot find a good R or Python package for this.

So far, I've found a single function in Statsmodels and another one in a package called Pingouin. But neither of these seem to cover the case with proportional bias or heteroskedastic differences.

Is there an R or Python package (or any other software) out there that does this out of the box? That is, test for proportional bias and heteroskedasticity and make the appropriate plot based on these tests.

Answer 1

You could regress differences on means, linear or LOESS. If the slope is markedly different from zero, there is a proportional bias. In addition a Breusch-Pagan test for heteroskedasticity. Here an R approach:

bland_altman <- function(x, y, alpha=.05, plot=TRUE, legend=TRUE, print=TRUE) {
  ## Function to create Bland-Altman plot
  z <- qt(p=alpha/2, df=Inf)
  ## Calculate means and differences
  means <- (x + y)/2
  differences <- x - y
  mdif <- mean(differences)
  mdif_se <- sd(differences)/sqrt(length(differences))
  mdif_ci <- mdif + c(-z, z)*sd(differences)
  dif <- c(mean=mdif, se=mdif_se, lower=mdif_ci[1], upper=mdif_ci[2])
  ## Test for proportional bias
  lm_bias <- lm(differences ~ means)
  bias_test <- cbind(summary(lm_bias)$coe, confint(lm_bias))
  ## Test for heteroskedasticity
  heteroskedasticity_test <- lmtest::bptest(lm_bias)
  if (plot) {
    ## Plot
    lclr <- 'grey62'
    pclr <- 'black'
    plot(means, differences, pch=1, col=pclr,
         ylim=range(c(differences, mdif_ci)),
         xlab="Mean", 
         ylab="Difference",
         main="Bland-Altman Plot", 
         panel.first={
           abline(h=mdif, col=lclr)  ## Add mean line
           abline(h=mdif_ci, col=lclr, lty=2)
         })
    abline(lm_bias, col="green", lty=1)
    ## Use a loess smoothing line to show varying spread
    smooth_diff <- loess(differences ~ means)
    pred_diff <- predict(smooth_diff, se=TRUE)
    ord <- order(means)
    lines(means[ord], pred_diff$fit[ord], col='red')
    lines(means[ord], pred_diff$fit[ord] + z*pred_diff$se.fit[ord],
          col="red", lty=2)
    lines(means[ord], pred_diff$fit[ord] - z*pred_diff$se.fit[ord], 
          col="red", lty=2)
    if (legend) {
      legend("topright", legend=c('Mean Difference', 'Prop. bias linear', 
                                  'Prop. bias LOESS'),
             col=c(lclr, "green", "red"), lty=1, cex=.85)
    }
  }
  ## Add legend
  if (print) {
    ## Print results
    message("Mean difference:\n", appendLF=FALSE)
    print(dif)
    message("\nProportional Bias Test:\n", appendLF=FALSE)
    print(bias_test)
    message("\nHeteroskedasticity Test:\n", appendLF=FALSE)
    print(heteroskedasticity_test) 
  }
  return(invisible(list(dif=dif, bias=lm_bias, 
                        heteroskedasticity=heteroskedasticity_test)))
}

Usage

> ## w/o heteroskedasticity
> bland_altman(x, y)
Mean difference:
      mean         se      lower      upper 
  0.903154   1.316731  26.710617 -24.904309 

Proportional Bias Test:
               Estimate Std. Error    t value  Pr(>|t|)       2.5 %     97.5 %
(Intercept)  6.75340556  18.470893  0.3656242 0.7154335 -29.9014806 43.4082917
means       -0.05814348   0.183104 -0.3175434 0.7515069  -0.4215075  0.3052205

Heteroskedasticity Test:

    studentized Breusch-Pagan test

data:  lm_bias
BP = 0.0017019, df = 1, p-value = 0.9671

> ## w/ heteroskedasticity
> bland_altman(x, y_het)
Mean difference:
      mean         se      lower      upper 
 -2.922508   2.313701  42.425189 -48.270206 

Proportional Bias Test:
              Estimate Std. Error   t value     Pr(>|t|)     2.5 %      97.5 %
(Intercept) 115.693668 8.88446507  13.02202 4.238331e-23 98.062737 133.3246003
means        -1.156889 0.08560917 -13.51361 4.061192e-24 -1.326777  -0.9870001

Heteroskedasticity Test:

    studentized Breusch-Pagan test

data:  lm_bias
BP = 6.8571, df = 1, p-value = 0.008829

The function invisibly returns the test results, you may deactivate plotting and printing.

> res <- bland_altman(x, y, plot=FALSE, print=FALSE)
> res
$dif
      mean         se      lower      upper 
  0.903154   1.316731  26.710617 -24.904309 

$bias

Call:
lm(formula = differences ~ means)

Coefficients:
(Intercept)        means  
    6.75341     -0.05814  


$heteroskedasticity

    studentized Breusch-Pagan test

data:  lm_bias
BP = 0.0017019, df = 1, p-value = 0.9671

---

Data:

> set.seed(51676818)
> n <- 100
> x <- rnorm(n, 100, 10)
> y <- rnorm(n, 101, 10)
> y_het <- x + rnorm(n, 0, sd=.2*x)  ## Increasing standard deviation with x