plot()与ggplot2绘制Passing-Bablok回归置信区间的差异原因及解决方案

我在R语言中尝试分别使用plot()和ggplot2绘制Passing-Bablok回归直线的置信区间（CI），但发现两种方法呈现的置信区间存在显著差异。以下是我使用的R代码：

# Loading libraries
library(ggplot2)
library(MethComp)

# Hypothethical data set for stackoverflow
set.seed(123) # Set seed for reproducibility
data_for_plot <- data.frame(rater1 = rnorm(100, mean = 50, sd = 10))
data_for_plot$rater2 <- data_for_plot$rater1 + runif(100, min = -20, max = 20)

# Prepare for x and y axis coordinates
common_min <- min(c(data_for_plot$rater1, data_for_plot$rater2))
common_min <- floor(common_min / 5) * 5
common_max <- max(c(data_for_plot$rater1, data_for_plot$rater2))
common_max <- ceiling(common_max / 5) * 5
common_breaks <- seq(from = common_min, to = common_max, length.out = 5)

# Passing-Bablok regression
pb_result <- PBreg(data_for_plot$rater1, data_for_plot$rater2)
pb_intercept <- pb_result$coefficients["Intercept", "Estimate"]
pb_intercept_lci <- pb_result$coefficients["Intercept", "2.5%CI"]
pb_intercept_uci <- pb_result$coefficients["Intercept", "97.5%CI"]
pb_slope <- pb_result$coefficients["Slope", "Estimate"]
pb_slope_lci <- pb_result$coefficients["Slope", "2.5%CI"]
pb_slope_uci <- pb_result$coefficients["Slope", "97.5%CI"]

# Manually creating CI ribbon (Problematic segment)
x_vals <- seq(from = common_min, to = common_max, length.out = 100)
upper_bound <- pb_intercept_uci + pb_slope_uci * x_vals
upper_bound <- pmin(upper_bound, common_max)
lower_bound <- pb_intercept_lci + pb_slope_lci * x_vals
lower_bound <- pmax(lower_bound, common_min)
ribbon_data <- data.frame(
  x = x_vals,
  ymin = lower_bound,
  ymax = upper_bound
)

# Generating plots
gg_plot_result <- ggplot(
  data = data_for_plot, aes(x = rater1, y = rater2)) +
  geom_ribbon(data = ribbon_data, aes(x = x, ymin = ymin, ymax = ymax), 
              fill = "#253494", alpha = 0.2, inherit.aes = FALSE)+
  geom_abline(intercept = 0, slope = 1, 
              color = "orange", size = 0.25, linetype="dashed") +
  geom_abline(intercept = pb_intercept, slope = pb_slope, 
              color = "#253494", linewidth = 1) +  # Plotting the regression line
  geom_point(shape = 21, colour = "black", fill = "white", size = 2.5) +
  labs(title = "Measurement (cm)", x = "Rater 1", y = "Rater 2") +           
  theme_minimal() +                                                      
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
    plot.title = element_text(hjust = 0.5),
    aspect.ratio = 1
  ) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(limits = c(common_min, common_max), breaks = common_breaks) +
  scale_y_continuous(limits = c(common_min, common_max), breaks = common_breaks)

# Generating PB plot directly
pb_plot_result <- plot(
  pb_result, asp = 1,
  xlim = c(common_min, common_max), 
  ylim = c(common_min, common_max),
  main = "Measurement (cm)",
  xlab = "Rater 1",
  ylab = "Rater 2"
)

print(gg_plot_result)

使用plot()生成的图中置信区间范围相对较小，而通过手动计算截距和斜率的置信区间并添加到ggplot2中的图里，置信区间范围明显更大。希望有人能解释差异产生的原因，同时请教其他可基于截距和斜率的置信区间绘制回归置信区间的方法。

差异原因解释 #

• 置信区间类型不同：plot(PBreg)绘制的是回归直线的联合置信区间，它考虑了截距和斜率之间的相关性；而你手动计算的是截距、斜率各自的边际置信区间的简单组合，没有考虑两者的相关性，导致范围被过度放大。
• MethComp的实现逻辑：plot.PBreg内部是基于Passing-Bablok方法的分布特性，直接计算回归直线的联合置信带，并非使用截距和斜率的边际CI拼接。

正确的ggplot2绘制方法 #

方法1：提取plot.PBreg的内置置信带数据 #

这是最简便且准确的方式，直接复用MethComp包已经计算好的联合置信带数据：

# 运行plot但不显示，提取绘图所需的置信带数据
pb_plot_data <- plot(pb_result, plot.it = FALSE)
# 整理置信带数据框
ci_data <- data.frame(
  x = pb_plot_data$x,
  ymin = pb_plot_data$lcl,
  ymax = pb_plot_data$ucl
)

# 绘制正确的ggplot图
gg_plot_correct <- ggplot(data_for_plot, aes(x = rater1, y = rater2)) +
  geom_ribbon(data = ci_data, aes(x = x, ymin = ymin, ymax = ymax), 
              fill = "#253494", alpha = 0.2, inherit.aes = FALSE) +
  geom_abline(intercept = 0, slope = 1, color = "orange", size = 0.25, linetype = "dashed") +
  geom_abline(intercept = pb_intercept, slope = pb_slope, color = "#253494", linewidth = 1) +
  geom_point(shape = 21, colour = "black", fill = "white", size = 2.5) +
  labs(title = "Measurement (cm)", x = "Rater 1", y = "Rater 2") +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.5),
    plot.title = element_text(hjust = 0.5),
    aspect.ratio = 1
  ) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(limits = c(common_min, common_max), breaks = common_breaks) +
  scale_y_continuous(limits = c(common_min, common_max), breaks = common_breaks)

print(gg_plot_correct)

方法2：手动计算联合置信带（进阶） #

如果需要手动实现，需参考Passing-Bablok回归的联合置信区间计算逻辑，利用截距和斜率的协方差矩阵来构建联合置信带。不过由于MethComp已经封装了这一过程，推荐直接使用方法1。

内容的提问来源于stack exchange，提问作者Arthur Yong