Do this yourself:

<code>
n = 50
b0 = .0823 # completely made up parameters
b1 = 0.153
s = 5.88
x = runif(n, min=12, max=222)
y = b0 + b1*x + rnorm(n,0,s)

fit = glm(y~x)

plot(x,y,xlab='Proxy', ylab="Temperature",ylim=c(-10,40),xlim=c(0,250))
   abline(fit,lty=2,lwd=2)

# new data
nb = 10
x_new = sort(runif(nb, min=12, max=222))
a = summary(fit)

plot(x_new, a$coefficients[1,1] + a$coefficients[2,1]*x_new, xlab="New Proxy", ylab="New Temperature",pch=20,lwd=2,col=4,cex=1.5,ylim=c(-10,40),xlim=c(0,250))

# routine for crude, too narrow Bayesian prediction interval; actual intervals are T
# obs.glm() routine is available in the class notes book at the left; <a href="http://wmbriggs.com/book/">or here</a>.
y_lo = 0
y_hi = 0
for (i in 1:nb){
   newdata = data.frame(x = x_new[i])
   s1 = obs.glm(fit,newdata)
   y_lo[i] = qnorm(0.025,s1$central,s1$scale) # approx; see glm.obs.prob
   y_hi[i] = qnorm(0.975,s1$central,s1$scale) # approx; see glm.obs.prob

}
polygon(c(x_new,x_new[nb:1]),c(y_lo,y_hi[nb:1]) , col="#ebeeca")

# classical parametric bounds
    newdata = data.frame(x = x_new)
    a = predict.glm(fit,newdata,type='response',se.fit=TRUE)
    y_lo = a$fit - 1.96*a$se.fit
    y_hi = a$fit + 1.96*a$se.fit
polygon(c(x_new,x_new[nb:1]),c(y_lo,y_hi[nb:1]) , col="#9b9e8a")

# over plot old data
    a = summary(fit)
    points(x_new, a$coefficients[1,1] + a$coefficients[2,1]*x_new, pch=20,lwd=2,col=4,cex=1.5)
    points(x,y)
</code>