如何给R语言基础绘图添加置信区间？或转用ggplot实现

一、先明确：计算置信区间的核心前提 #

你当前的p1-p4是基于系数点估计计算的预测值，要添加置信区间，必须先获取这些系数（比如-1.2578、-0.005104等）的协方差矩阵/标准误——这通常可以从你的回归模型（看起来是有序logit模型）中通过vcov(model)提取。下面以有序logit模型为例，演示完整的计算和绘图流程。

二、步骤1：计算p1-p4的置信区间（模拟法推荐） #

模拟法适合处理非线性转换（比如累积logistic分布的差值），结果更稳定：

library(mvtnorm)

# 1. 替换为你实际模型的系数和协方差矩阵
# 示例：用你给出的系数构造（实际中替换为coef(model)和vcov(model)）
coefs <- c(theta1 = -1.2578, theta2 = -0.5250, theta3 = 0.6878, beta = -0.005104)
# 示例标准误，实际替换为模型输出的标准误
vcov_mat <- diag(c(0.1^2, 0.1^2, 0.1^2, 0.001^2))

# 2. 模拟1000组系数，模拟系数的不确定性
set.seed(123)
sim_coefs <- rmvnorm(1000, mean = coefs, sigma = vcov_mat)

# 3. 对每个deficit值，计算模拟的p1-p4
deficit <- seq(from=0, to=300, by=10)
n_deficit <- length(deficit)
sim_p <- matrix(NA, nrow = 1000, ncol = 4*n_deficit)

logistic_cdf <- function(x) 1/(1+exp(-x))

for(i in 1:1000){
  theta1_sim <- sim_coefs[i,1]
  theta2_sim <- sim_coefs[i,2]
  theta3_sim <- sim_coefs[i,3]
  beta_sim <- sim_coefs[i,4]
  xbeta_sim <- deficit * beta_sim
  
  p1_sim <- logistic_cdf(theta1_sim - xbeta_sim)
  p2_sim <- logistic_cdf(theta2_sim - xbeta_sim) - logistic_cdf(theta1_sim - xbeta_sim)
  p3_sim <- logistic_cdf(theta3_sim - xbeta_sim) - logistic_cdf(theta2_sim - xbeta_sim)
  p4_sim <- 1 - logistic_cdf(theta3_sim - xbeta_sim)
  
  sim_p[i,] <- c(p1_sim, p2_sim, p3_sim, p4_sim)
}

# 4. 计算95%置信区间（分位数法）
ci_p1 <- t(apply(sim_p[,1:n_deficit], 2, quantile, c(0.025, 0.975)))
ci_p2 <- t(apply(sim_p[,(n_deficit+1):(2*n_deficit)], 2, quantile, c(0.025, 0.975)))
ci_p3 <- t(apply(sim_p[,(2*n_deficit+1):(3*n_deficit)], 2, quantile, c(0.025, 0.975)))
ci_p4 <- t(apply(sim_p[,(3*n_deficit+1):(4*n_deficit)], 2, quantile, c(0.025, 0.975)))

三、基础绘图添加置信区间 #

用polygon()绘制置信区间阴影：

# 初始化绘图
plot(deficit, p1, type='l', 
     ylab='Probability of nuisance level',
     xlab= "Precipitation deficit",
     main= "(b)",
     lwd=2.5, col= "#ffd47f",
     ylim = c(0, 1)) # 调整y轴范围容纳置信区间

# 依次添加各分组的置信区间和折线
# p1
polygon(c(deficit, rev(deficit)), c(ci_p1[,1], rev(ci_p1[,2])), 
        col=adjustcolor("#ffd47f", alpha.f=0.3), border=NA)
lines(deficit, p1, lwd=2.5, col= "#ffd47f")

# p2
polygon(c(deficit, rev(deficit)), c(ci_p2[,1], rev(ci_p2[,2])), 
        col=adjustcolor("#fbfb7c", alpha.f=0.3), border=NA)
lines(deficit, p2, col='#fbfb7c', lwd=2.5)

# p3
polygon(c(deficit, rev(deficit)), c(ci_p3[,1], rev(ci_p3[,2])), 
        col=adjustcolor("#cfe7ef", alpha.f=0.3), border=NA)
lines(deficit, p3, col="#cfe7ef", lwd=2.5)

# p4
polygon(c(deficit, rev(deficit)), c(ci_p4[,1], rev(ci_p4[,2])), 
        col=adjustcolor("#75a5e9", alpha.f=0.3), border=NA)
lines(deficit, p4, col='#75a5e9', lwd=2.5)

# 添加图例
legend("topleft", title = "Nuisance levels", lty=1, lwd=2, 
       col=c("#ffd47f", '#fbfb7c', "#cfe7ef", '#75a5e9'), 
       legend=c("No nuisance", "Low nuisance levels", "High nuisance levels", "Very high Nuisance levels"))

四、替代方案：用ggplot2实现带置信区间的绘图 #

先整理成长格式数据框，再用geom_ribbon绘制置信区间：

library(ggplot2)

# 整理数据框
df <- data.frame(
  deficit = rep(deficit, 4),
  prob = c(p1, p2, p3, p4),
  lower = c(ci_p1[,1], ci_p2[,1], ci_p3[,1], ci_p4[,1]),
  upper = c(ci_p1[,2], ci_p2[,2], ci_p3[,2], ci_p4[,2]),
  level = rep(c("No nuisance", "Low nuisance levels", "High nuisance levels", "Very high Nuisance levels"), each = length(deficit))
)

# 绘图
ggplot(df, aes(x = deficit, y = prob, color = level)) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = level), alpha = 0.3, color = NA) +
  geom_line(linewidth = 2.5) +
  labs(
    x = "Precipitation deficit",
    y = "Probability of nuisance level",
    title = "(b)",
    color = "Nuisance levels",
    fill = "Nuisance levels"
  ) +
  scale_fill_manual(values = c("#ffd47f", '#fbfb7c', "#cfe7ef", '#75a5e9')) +
  scale_color_manual(values = c("#ffd47f", '#fbfb7c', "#cfe7ef", '#75a5e9')) +
  ylim(c(0, 1)) +
  theme_bw()

关键提示 #

如果没有模型的协方差矩阵，无法计算置信区间——因为置信区间的本质是量化系数估计的不确定性对预测概率的影响。模拟法相比delta方法更适合这类非线性转换场景，结果偏差更小。

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