如何使用R语言绘制带预测区间和置信区间的回归线

先修正现有代码的问题 #

• 变量命名颠倒：你调用predict()时interval = "confidence"返回的是置信区间，被你赋值给了pred_interval；interval="prediction"返回的是预测区间，被你赋值给了conf_interval，二者含义和命名完全相反
• 绘制区间前需要先按自变量Latitude排序，否则绘制的区间多边形会出现路径错乱

基础R绘图完整实现 #

# 拟合线性模型
lm_fit <- lm(Number.of.species ~ Latitude, data = dta)
# 计算99%置信区间和预测区间
ci <- predict(lm_fit, level = 0.99, interval = "confidence")
pi <- predict(lm_fit, level = 0.99, interval = "prediction")
# 按纬度排序，避免绘图错乱
dta_sorted <- dta[order(dta$Latitude), ]
ci_sorted <- ci[order(dta$Latitude), ]
pi_sorted <- pi[order(dta$Latitude), ]

# 布置画布，2行1列分别放两张图
par(mfrow = c(2,1), mar = c(4,4,2,1))

# 第一张：带99%预测区间的回归线
plot(
  dta$Latitude, 
  dta$Number.of.species, 
  pch = 1, 
  ylim = c(0, 180), 
  xlim = c(37, 40),
  main = "回归线 + 99%预测区间",
  xlab = "纬度",
  ylab = "物种数量"
)
# 绘制预测区间半透明填充
polygon(
  x = c(dta_sorted$Latitude, rev(dta_sorted$Latitude)),
  y = c(pi_sorted[,2], rev(pi_sorted[,3])),
  col = rgb(0.2,0.6,0.8,0.3),
  border = NA
)
# 加回归线
abline(lm_fit, lwd = 2, col = "darkblue")

# 第二张：带99%置信区间的回归线
plot(
  dta$Latitude, 
  dta$Number.of.species, 
  pch = 1, 
  ylim = c(0, 180), 
  xlim = c(37, 40),
  main = "回归线 + 99%置信区间",
  xlab = "纬度",
  ylab = "物种数量"
)
# 绘制置信区间半透明填充
polygon(
  x = c(dta_sorted$Latitude, rev(dta_sorted$Latitude)),
  y = c(ci_sorted[,2], rev(ci_sorted[,3])),
  col = rgb(0.8,0.2,0.2,0.3),
  border = NA
)
# 加回归线
abline(lm_fit, lwd = 2, col = "darkred")

ggplot2 简化实现（可选） #

如果你可以使用ggplot2包，代码会更简洁，不需要手动排序和绘制多边形：

library(ggplot2)
# 带预测区间的图
ggplot(dta, aes(x = Latitude, y = Number.of.species)) +
  geom_point() +
  geom_smooth(method = "lm", level = 0.99, interval = "prediction", fill = "lightblue") +
  coord_cartesian(xlim = c(37,40), ylim = c(0,180)) +
  labs(title = "回归线 + 99%预测区间", x = "纬度", y = "物种数量")

# 带置信区间的图
ggplot(dta, aes(x = Latitude, y = Number.of.species)) +
  geom_point() +
  geom_smooth(method = "lm", level = 0.99, interval = "confidence", fill = "pink") +
  coord_cartesian(xlim = c(37,40), ylim = c(0,180)) +
  labs(title = "回归线 + 99%置信区间", x = "纬度", y = "物种数量")

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