如何在ggplot中填充置信区间曲线区域（解决曲线交叉问题）

我想用ggplot填充两条函数曲线之间的区域来制作自定义置信区间。我用Monod方程构建了nls拟合函数，基于拟合参数绘制了拟合曲线；同时通过bootstrap方法得到置信区间的数值，绘制了两条曲线作为区间边界，但没能成功填充区域，还发现边界曲线存在交叉问题。请问有没有可行的填充方法，或者更合适的置信区间绘制方式？

代码 #

pr20Y = subset(donnees_tot_g_J4, color == "20Y")

start.values <- list(gm=0.7, s=0.05, k=0.2)
nls1 = nls(growth_rate ~ gm*((quantity-s)/(quantity-s+k)), data = pr20Y, start=start.values)  
summary(nls1)
# gm = 0.65 / s = 0.06 / k = 0.11
nls1.boot <- nlsBoot(nls1, niter = 1000)
summary(nls1.boot)
    # Median       2.5%      97.5%
# gm 0.65636222 0.61994560 0.69597337
# s  0.05735049 0.03765189 0.07138767
# k  0.11705112 0.07833971 0.16174788
growthJ4=ggplot(donnees_tot_g_J4, aes(x=quantity, y=growth_rate, col=color))+
  geom_point(size=3,aes(col=color))+
  stat_function(fun = function (quantity) 0.65*((quantity-0.06)/(quantity-0.06+0.11)), color="darkgoldenrod1", size=1)+
  stat_function(fun = function (quantity) 0.61994560*((quantity-0.03765189)/(quantity-0.03765189+0.07833971)), color="darkgoldenrod1", size=1, alpha=0.5)+
  stat_function(fun =function (quantity) 0.69597337*((quantity-0.07138767)/(quantity-0.07138767+0.16605772)), color="darkgoldenrod1", size=1, alpha=0.5)+
  scale_x_continuous(breaks=c(0.1,0.3,0.6,0.9,1.5), limits=c(0.1,1.5))+
  scale_y_continuous(limits=c(0,0.8))+
  scale_colour_manual(values=c("20S" = "aquamarine1","25S" = "aquamarine3","28S" = "aquamarine4","20Y" = "darkgoldenrod1","25Y" = "darkgoldenrod3", "28Y" = "darkgoldenrod4"))+
  scale_fill_manual(values=c("20S" = "aquamarine1","25S" = "aquamarine3","28S" = "aquamarine4","20Y" = "darkgoldenrod1", "25Y" = "darkgoldenrod3","28Y" = "darkgoldenrod4"))+
  theme_minimal()+
  ggtitle("Taux de croissance à J4 en fonction du traitement de nourriture et de la température")+
  xlab("Quantité") + ylab("Taux de croissance µg/j")+
  theme_grey(base_size = 22)
growthJ4

数据 #

structure(list(name = c("J4_S01AC", "J4_S01CC", "J4_S01EC", "J4_S01FC", 
"J4_S03BC", "J4_S03CC", "J4_S03DC", "J4_S03EC", "J4_S03FC", "J4_S06AC", 
"J4_S06DC", "J4_S06EC", "J4_S06FC", "J4_S06KC", "J4_S06MC", "J4_S06NC", 
"J4_S09AC", "J4_S09BC", "J4_S09CC", "J4_S09DC", "J4_S03AM", "J4_S03BM", 
"J4_S03CM", "J4_S15AM", "J4_S15BM", "J4_S15CM", "J4_S15DM", "J4_S01AAC"
), day = c("J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", 
"J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4", 
"J4", "J4", "J4", "J4", "J4", "J4", "J4", "J4"), quality = c("S", 
"S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", 
"S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", 
"S"), quantity = c(0.1, 0.1, 0.1, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 
0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.9, 0.9, 0.9, 0.9, 0.3, 0.3, 
0.3, 1.5, 1.5, 1.5, 1.5, 0.1), qual_quant = c("S01", "S01", "S01", 
"S01", "S03", "S03", "S03", "S03", "S03", "S06", "S06", "S06", 
"S06", "S06", "S06", "S06", "S09", "S09", "S09", "S09", "S03", 
"S03", "S03", "S15", "S15", "S15", "S15", "S01"), temperature = c(28, 
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 
28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28), time = c(101, 101, 
101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 101, 
101, 102, 102, 102, 102, 107.5, 107.5, 107.5, 107.5, 107.5, 107.5, 
107.5, 109), size = c(1344.366, 1291.962, 1304.811, 1128.264, 
1286.151, 1358.421, 1396.671, 1353.076, 1505.565, 1297.17, 0, 
1243.029, 1323.85, 1368.364, 1506.396, 0, 1663.735, 1632.28, 
2115.303, 1921.46, 1506.581, 1501.196, 1370.99, 1870.489, 1941.425, 
1942.186, 1827.395, 1336.588), weight = c(11, 16, 10, 10, 12, 
12, 13, 16, 14, 12, 11, 12, 10, 10, 15, 25, 46, 35, 66, 46, 20, 
16, 15, 49, 73, 63, 60, 11), growth_rate = c(0.0378568479840844, 0.126892202895244, 
0.0152090656484838, 0.0152090656484838, 0.0585326533474545, 0.0585326533474545, 
0.0775525505364441, 0.126892202895244, 0.0951622248242906, 0.0585326533474545, 
0.0378568479840844, 0.0585326533474545, 0.0152090656484838, 0.0152090656484838, 
0.111556439370444, 0.232939774941222, 0.374129156018743, 0.309825356331568, 
0.459072793066665, 0.374129156018743, 0.169037347727828, 0.119219651092275, 
0.104811166292231, 0.369092608582107, 0.458090402952369, 0.425199566970549, 
0.414306966296783, 0.0350783637283718), color = c("28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S", 
"28S", "28S", "28S", "28S", "28S", "28S", "28S", "28S")), row.names = c(NA, 
-28L), class = c("tbl_df", "tbl", "data.frame"))

解决方案 #

1. 解决曲线交叉问题 #

你当前直接用参数的2.5%和97.5%分位点组合成边界曲线，忽略了参数间的相关性，这是导致曲线交叉的核心原因。正确做法是利用bootstrap得到的所有参数样本，对每个x值计算对应的y值，再取每个x的2.5%和97.5%分位点作为置信区间的上下限，保留参数间的关联，避免曲线交叉。

2. 填充置信区间区域 #

ggplot无法直接用stat_function填充两条曲线间的区域，需要先生成包含x、拟合值、置信上下限的数据集，再用geom_ribbon完成填充。

具体实现代码

library(ggplot2)
library(nlstools) # 用于nlsBoot

# 提取bootstrap参数样本
boot_params <- as.data.frame(nls1.boot$coefboot)

# 生成覆盖绘图范围的x序列
x_seq <- seq(0.1, 1.5, length.out = 100)

# 定义Monod方程函数
monod_fun <- function(x, gm, s, k) {
  gm * ((x - s) / (x - s + k))
}

# 计算每个bootstrap样本在x_seq上的预测值
boot_preds <- apply(boot_params, 1, function(params) {
  monod_fun(x_seq, gm = params["gm"], s = params["s"], k = params["k"])
})

# 计算每个x对应的拟合中位数、2.5%和97.5%分位点
pred_df <- data.frame(
  quantity = x_seq,
  fit = apply(boot_preds, 1, median),
  lower = apply(boot_preds, 1, quantile, probs = 0.025),
  upper = apply(boot_preds, 1, quantile, probs = 0.975)
)

# 绘制最终图形
growthJ4 <- ggplot() +
  # 原始散点
  geom_point(data = donnees_tot_g_J4, aes(x = quantity, y = growth_rate, col = color), size = 3) +
  # 置信区间填充区域
  geom_ribbon(data = pred_df, aes(x = quantity, ymin = lower, ymax = upper), fill = "darkgoldenrod1", alpha = 0.3) +
  # 拟合曲线
  geom_line(data = pred_df, aes(x = quantity, y = fit), color = "darkgoldenrod1", size = 1) +
  # 坐标轴与主题设置
  scale_x_continuous(breaks = c(0.1, 0.3, 0.6, 0.9, 1.5), limits = c(0.1, 1.5)) +
  scale_y_continuous(limits = c(0, 0.8)) +
  scale_colour_manual(values = c("20S" = "aquamarine1","25S" = "aquamarine3","28S" = "aquamarine4","20Y" = "darkgoldenrod1","25Y" = "darkgoldenrod3", "28Y" = "darkgoldenrod4")) +
  ggtitle("Taux de croissance à J4 en fonction du traitement de nourriture et de la température") +
  xlab("Quantité") + ylab("Taux de croissance µg/j") +
  theme_grey(base_size = 22)

growthJ4

关键说明 #

• 通过遍历所有bootstrap参数样本计算预测值，再按x维度取分位点，确保置信区间基于参数的联合分布，不会出现交叉。
• geom_ribbon是ggplot中专门用于填充区间的图层，设置alpha参数可以让填充区域不遮挡散点和拟合曲线。
• 注意：当前提供的样本数据中仅包含28S组，若实际数据中无20Y组，需要补充对应数据才能得到有效拟合结果。

内容的提问来源于stack exchange，提问作者Nate Trf