Do this yourself: ``` 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; or here. 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) ```