Shading confidence intervals in R - base R if possible

Question

I am comparing two lines that were regressed using LOESS. I want to display the confidence intervals of the two lines clearly, and I am having some difficulties.

I have tried using a variety of line types and colors, but the result is still to busy and messy in my opinion. I think that shading between the confidence intervals might make things clearer, but I am having some difficulties wrapping my head around the problem considering how my coding is structured so far. I have included the produced plot, the data for the two sets Analysis5k and Analysis5kz, and my code up to this point.

I have seen some examples where two polygons were overlapped to show where the confidence intervals overlap which seems like it might be a good way to present the data. If there was a way to draw the polygon in the area that is shared by the two confidence intervals, that might be another good way to present the data.

I understand the basic concept of how the polygon should be done, but the examples I have found have been applied to much more simplistic lines and data. Part of this is my own fault for some poor organization up to this point, but as this step is basically the finishing touch on my data presentation, I really don't want to have to rework everything from ground up.

Any help or insight is greatly appreciated.

UPDATE

I updated the Title. I received some great examples using ggplot, and while I would like to get around to working with ggplot in the future, I have only dealt with base R up to this point. For this particular project would like to try to keep this in base R if possible.

Analysis5k

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

Analysis5kz

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

Code

  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

Here is a solution with base plot based on your code.

The trick with polygon is that you must provide 2 times the x coordinates in one vector, once in normal order and once in reverse order (with function rev) and you must provide the y coordinates as a vector of the upper bounds followed by the lower bounds in reverse order.

We use the adjustcolor function to make standard colors transparent.

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)