如果可能的话，对 R 中的置信区间进行着色 - 以 R 为基础

Question

我正在比较使用 LOESS 回归的两条线。我想清楚地显示两条线的置信区间，但遇到了一些困难。

我尝试过使用各种线条类型和颜色，但在我看来，结果仍然是忙碌和混乱。我认为置信区间之间的阴影可能会让事情变得更清晰，但考虑到到目前为止我的编码结构，我在解决这个问题时遇到了一些困难。我已经包含了生成的图、两组 Analysis5k 和 Analysis5kz 的数据以及到目前为止的代码。

我见过一些示例，其中两个多边形重叠以显示置信区间重叠的位置，这似乎可能是呈现数据的好方法。如果有一种方法可以在两个置信区间共享的区域中绘制多边形，那么这可能是呈现数据的另一种好方法。

我了解多边形应该如何完成的基本概念，但我发现的示例已应用于更简单的线条和数据。到目前为止，部分原因是我自己对一些糟糕的组织的错误，但由于这一步基本上是我的数据演示的最后一步，我真的不想从头开始重新设计所有内容。

非常感谢任何帮助或见解。

更新

我更新了标题。我收到了一些使用 ggplot 的很好的例子，虽然我想在将来开始使用 ggplot，但到目前为止我只处理过基础 R。对于这个特定的项目，如果可能的话，希望尝试将其保留在基础 R 中。

分析5k

Period  15p5    Total_5plus
-4350   0.100529101 12.6
-3900   0.4 20
-3650   0.0625  9.6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3650   0.174757282 20.6
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-3400   0.208333333 24
-4233   0.184027778 19.2
-3650   0.285714286 12.6
-4350   0.166666667 6

分析5kz

Period  15p5    Total_5plus
-4350   0.100529101 12.6
-4350   0   5
-3900   0.4 20
-3650   0.0625  9.6
-3400   0   6
-3900   0.126984127 16.8
-3958   0.133333333 5
-4350   0.150943396 10.6
-3400   0.146341463 8.2
-3650   0.255319149 9.4
-3400   0.222222222 9
-3500   0.245014245 39
-3600   0.125   8
-3650   0   28
-3808   0.1 20
-3900   0.160493827 18
-3958   0.238095238 7
-4058   0.2 5
-3500   0   25
-3500   0.086956522 28.75
-4117   0.141414141 6.6
-4350   0.171038825 31.76666667
-4350   0.166666667 6
-3650   0.143798024 30.36666667
-2715   0.137931034 7.25
-4350   0.235588972 26.6
-3500   0.228840125 79.75
-4350   0.041666667 8
-3500   0   5
-3650   0.174757282 20.6
-3800   0   9
-2715   0.377777778 11.25
-3500   0.2 7.5
-3650   0.078947368 7.6
-4117   0   8
-4350   0   8
-3400   0.208333333 24
-4233   0.184027778 19.2
-3025   0   7
-3650   0.285714286 12.6
-4350   0.166666667 6

代码

  ppi <- 300 
  png("5+ KC shaded CI.png", width=6*ppi, height=6*ppi, res=ppi) 
  library(Hmisc) 
  Analysis5k <- read.csv(file.choose(), header = T) 
  Analysis5kz <- read.csv(file.choose(), header = T)
  par(mfrow = c(1,1), pty = "s", oma=c(1,2,1,1), mar=c(4,4,2,2)) 
  plot(X15p5 ~ Period, Analysis5kz, xaxt = "n", yaxt= "n", ylim=c(-0.2,0.7), xlim=c(-5000,-2500), xlab = "Years B.P.", ylab = expression(''[15]*'p'[5]), main = "") 
  vx <- seq(-5000,-2000, by = 500) 
  vy <- seq(-0.2,0.7, by = 0.1) 
  axis(1, at = vx) 
  axis(2, at = vy) 
  a5k <- order(Analysis5k$Period) 
  a5kz <- order(Analysis5kz$Period)
  Analysis5k.lo <- loess(X15p5 ~ Period, Analysis5k, weights = Total_5plus, span = 0.6) 
  Analysis5kz.lo <- loess(X15p5 ~ Period, Analysis5kz, weights = Total_5plus, span = 0.6)      
  pred5k <- predict(Analysis5k.lo, se = TRUE) 
  pred5kz <- predict(Analysis5kz.lo, se = TRUE)      
  lines(Analysis5k$Period[a5k], pred5k$fit[a5k], col="blue", lwd=2) 
  lines(Analysis5kz$Period[a5kz], pred5kz$fit[a5kz], col="skyblue", lwd=2)          
  lines(Analysis5K$Period[a5K], pred5K$fit[a5K] - qt(0.975, pred5K$df)*pred5K$se[a5K],col="blue",lty=2) 
  lines(Analysis5K$Period[a5K], pred5K$fit[a5K] + qt(0.975, pred5K$df)*pred5K$se[a5K],col="blue",lty=2)      
  lines(Analysis5Kz$Period[a5Kz], pred5Kz$fit[a5Kz] - qt(0.975, pred5Kz$df)*pred5Kz$se[a5Kz],col="skyblue",lty=2) 
  lines(Analysis5Kz$Period[a5Kz], pred5Kz$fit[a5Kz] + qt(0.975, pred5Kz$df)*pred5Kz$se[a5Kz],col="skyblue",lty=2)
  abline(h=0.173, lty=3) 
  abline(v=-4700, lty=3)
  abline(v=-4000, lty=3)
  abline(v=-3000, lty=3)
  minor.tick(nx=5, ny=4, tick.ratio=0.5) 
  dev.off()

Answer 1

这是一个基于您的代码的基本图解决方案。

诀窍polygon是，您必须在一个向量中提供 2 倍的 x 坐标，一次按正常顺序，一次按相反顺序（使用函数rev），并且您必须提供 y 坐标作为上界后跟下界的向量相反的顺序。

我们使用该adjustcolor函数使标准颜色透明。

library(Hmisc) 
ppi <- 300 
par(mfrow = c(1,1), pty = "s", oma=c(1,2,1,1), mar=c(4,4,2,2)) 
plot(X15p5 ~ Period, Analysis5kz, xaxt = "n", yaxt= "n", ylim=c(-0.2,0.7), xlim=c(-5000,-2500), xlab = "Years B.P.", ylab = expression(''[15]*'p'[5]), main = "") 
vx <- seq(-5000,-2000, by = 500) 
vy <- seq(-0.2,0.7, by = 0.1) 
axis(1, at = vx) 
axis(2, at = vy) 
a5k <- order(Analysis5k$Period) 
a5kz <- order(Analysis5kz$Period)
Analysis5k.lo <- loess(X15p5 ~ Period, Analysis5k, weights = Total_5plus, span = 0.6) 
Analysis5kz.lo <- loess(X15p5 ~ Period, Analysis5kz, weights = Total_5plus, span = 0.6)      
pred5k <- predict(Analysis5k.lo, se = TRUE) 
pred5kz <- predict(Analysis5kz.lo, se = TRUE)      

polygon(x = c(Analysis5k$Period[a5k], rev(Analysis5k$Period[a5k])),
        y = c(pred5k$fit[a5k] - qt(0.975, pred5k$df)*pred5k$se[a5k], 
              rev(pred5k$fit[a5k] + qt(0.975, pred5k$df)*pred5k$se[a5k])),
        col =  adjustcolor("dodgerblue", alpha.f = 0.10), border = NA)

polygon(x = c(Analysis5kz$Period[a5kz], rev(Analysis5kz$Period[a5kz])),
        y = c(pred5kz$fit[a5kz] - qt(0.975, pred5kz$df)*pred5kz$se[a5kz], 
              rev( pred5kz$fit[a5kz] + qt(0.975, pred5kz$df)*pred5kz$se[a5kz])),
        col =  adjustcolor("orangered", alpha.f = 0.10), border = NA)

lines(Analysis5k$Period[a5k], pred5k$fit[a5k], col="dodgerblue", lwd=2) 
lines(Analysis5kz$Period[a5kz], pred5kz$fit[a5kz], col="orangered", lwd=2)   

abline(h=0.173, lty=3) 
abline(v=-4700, lty=3)
abline(v=-4000, lty=3)
abline(v=-3000, lty=3)
minor.tick(nx=5, ny=4, tick.ratio=0.5)