Adding Error Bands to a Smooth in a Scatterplot

require(mgcv)
require(lattice)
data(Earthquake,package="MEMSS")
Earthquake$Magnitude <- equal.count(Earthquake$Richter, 3, overlap = 0.1)

coef <- coef(lm(log2(accel) ~ log2(distance), data = Earthquake))
panel.errorbar <- function(x,y,...)
  {
    fit <- gam(y~s(x))
    x0 <- seq(min(x),max(x),len=100)
    y0 <- predict(fit,newdata=data.frame(x=x0),se.fit=T)
    lpolygon(c(x0,rev(x0)),c(y0$fit-2*y0$se.fit,rev(y0$fit+2*y0$se.fit)),col="grey85",border=F,...)
    llines(x0,y0$fit,...)
  }
myPanel <- function(...)
  {
    panel.errorbar(...)
    panel.xyplot(...)
    panel.abline(reg = coef)
  }

xyplot(accel ~ distance | Magnitude, data = Earthquake,scales = list(log = 2),
       col.line = "grey", lwd = 2, xlab = "Distance From Epicenter (km)",
       ylab = "Maximum Horizontal Acceleration (g)",panel=myPanel,layout=c(3,1))

• Allow a fitting method to be specified by a function; any method that produces a fit with a predict method that accepts the se.fit = TRUE argument then includes an se.fit term in the result should work.
• Allow the user to specify how the band is represented (dashed lines or filled area, colors, etc.).