布兰德-奥特曼情节

Question

我需要制作 Bland-Altman 图，但我找不到合适的 R 或 Python 包。

到目前为止，我在 Statsmodels 中找到了一个函数，在名为 Pingouin 的包中找到了另一个函数。但这些似乎都没有涵盖具有比例偏差或异方差差异的情况。

是否有 R 或 Python 软件包（或任何其他软件）可以开箱即用地执行此操作？也就是说，测试比例偏差和异方差性，并根据这些测试绘制适当的图。

Answer 1

您可以回归均值、线性或 LOESS 的差异。如果斜率明显不同于零，则存在比例偏差。此外还进行了异方差性的 Breusch-Pagan 检验。像这样的东西：

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

使用方法

> ## 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

该功能隐形返回测试结果，您可以关闭绘图和打印。

> 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

---

数据：

> 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