cat("confbands: a function to plot confidence bands on lm, glm or nls predictions\n") confbands <- function(object, newdata, conf.level = 0.95, col = "black", alpha = 0.3, add = FALSE, ...) { # newdata: a dataframe with one named column (the x variable in model) n <- nrow(newdata) xname <- colnames(newdata) if (inherits(object, "nls")) { pred <- data.frame(newdata, Predict = rep(NA, n), SE = rep(NA, n)) for(k in 1:n) { pred[k, ] <- mydeltaMethod(object, xname = xname, xvalue = newdata[k, 1]) } DF <- summary(object)$df[2] } else if (inherits(object, "glm")) { predi <- predict.glm(object = object, type = "response", newdata = newdata, se.fit = TRUE) pred <- data.frame(newdata, Predict = predi$fit, SE = predi$se.fit) DF <- object$df.residual } else if (inherits(object, "lm")) { predi <- predict.lm(object = object, newdata = newdata, se.fit = TRUE) pred <- data.frame(newdata, Predict = predi$fit, SE = predi$se.fit) DF <- object$df.residual } ttab <- qt(1 - (1 - conf.level)/2, DF) pred$LCL <- pred$Predict - ttab * pred$SE pred$UCL <- pred$Predict + ttab * pred$SE # graph x <- newdata[, 1] if (add == FALSE) { plot(Predict ~ x, data = pred, type = "l", col = col, ...) } else { points(Predict ~ x, data = pred, type = "l", col = col) } polygon(x = c(x, rev(x)), y = c(pred$LCL, rev(pred$UCL)), col = adjustcolor(col, alpha = alpha), border = NA) invisible(pred) } # ------------------------------- # aux function: delta method mydeltaMethod <- function (object, xname, xvalue) { # object: of class nls g <- as.character(as.expression(formula(object)[[3L]])) g <- parse(text = g) vcov. <- vcov(object) object <- coef(object) para <- object para.names <- names(para) q <- length(para) for (i in 1:q) assign(para.names[i], para[i]) assign(xname, xvalue) est <- eval(g) names(est) <- NULL gd <- rep(0, q) for (i in 1:q) gd[i] <- eval(D(g, names(para)[i])) se.est <- as.vector(sqrt(t(gd) %*% vcov. %*% gd)) out <- data.frame(xvalue, Predict = est, SE = se.est) return(out) }