loess.demo <- function(x, y, span = 2/3, degree = 1, nearest = F, xlim=numeric(0), ylim=numeric(0), verbose=F) { # function to demonstrate the locally weighted regression function loess # written by Dr. Greg Snow # Brigham Young University, Department of Statistics # gls@byu.edu # Modified by Henrik Aa. Nielsen, IMM, DTU (han@imm.dtu.dk) miss.xy <- is.na(x) | is.na(y) x <- x[!miss.xy] y <- y[!miss.xy] y <- y[order(x)] x <- x[order(x)] fit.d <- loess(y ~ x, degree = degree, span = span, family = "gaussian", control = loess.control(surface = "direct")) fit.i <- loess(y ~ x, degree = degree, span = span, family = "gaussian") est <- list(x = seq(min(x), max(x), len=500)) est\$y <- predict(fit.i, newdata=data.frame(x = est\$x)) xl <- range(x, est\$x, xlim) xl <- xl + c(-1, 1) * 0.03 * diff(xl) yl <- range(y, est\$y, fitted(fit.d), ylim) yl <- yl + c(-1, 1) * 0.05 * diff(yl) fitPlot <- function(x, y, est, fit.d, xl, yl) { plot(x, y, pch = 3, xlim = xl, ylim = yl) lines(x, fitted(fit.d), col=2) mtext("Exact estimate with linear interpolation between x-values", col = 2, adj = 0.5, line = 0.5) lines(est, col=3) mtext("Estimate obtained using the default interpolation scheme", col = 3, adj = 0.5, line = 2) mtext NULL } fitPlot(x, y, est, fit.d, xl, yl) repeat { x0 <- locator(1)\$x if(length(x0) < 1) break if(nearest) x0 <- unique(x[abs(x - x0) == min(abs(x - x0))]) if(verbose) cat("x0 =", x0, "\n") if( span < 1 ){ q <- as.integer(span * length(x)) d <- sort(abs(x - x0))[q] } else { d <- max( abs(x - x0) ) * sqrt(span) } w <- rep(0, length(x)) s <- abs(x - x0) <= d w[s] <- (1 - (abs(x[s] - x0)/d)^3)^3 fitPlot(x, y, est, fit.d, xl, yl) symbols(x, y, circle = sqrt(w), inches = 0.3, add = T, col= 6) if(degree > 0) lines(x, fitted(lm(y ~ poly(x, degree), weights = w)), col = 8, err = -1) else { ##lines(x, fitted(lm(y ~ 1, weights = w)), col = 8, err = -1) abline(a = sum(w*y)/sum(w), b = 0, col = 8) } abline(v = x0, col = 8) if(x0 - d > xl[1]) abline(v = x0 - d, col = 8, lty = 2) if(x0 + d < xl[2]) abline(v = x0 + d, col = 8, lty = 2) } }