rm(list = ls()) library(rpanel) library(splines) library(mgcv) library(caret) #library(lattice) #library(latticeExtra) #----------------------------------------------------------------------- # (Global) Polynomial regression. rp.poly <- function(y, x, er = 0) { annotations <- function(m0) { mtext(side = 3, adj = 0, line = 1.5, text = sprintf("X rank: %i", m0$rank)) mtext(side = 3, adj = 0, line = 0.5, text = sprintf("Residual DF: %i", m0$df.residual)) press <- sum((residuals(m0)/(1 - hatvalues(m0)))^2) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("PRESS: %0.2f", press)) sm0 <- summary(m0) lastcoef <- sprintf("High term p-value: %0.5f", sm0$coeff[length(coef(m0)), 4]) mtext(side = 3, adj = 1, line = 2.5, text = lastcoef) mtext(side = 3, adj = 1, line = 1.5, text = sprintf("R² (adj. R²): %0.2f (%0.2f)", 100 * sm0$r.squared, 100 * sm0$adj.r.squared)) } draw.poly <- function(panel) { with(panel, { m0 <- lm(y ~ poly(x, degree = degree)) xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) xx <- seq(xr[1], xr[2], length.out = 201) cb <- predict(m0, newdata = data.frame(x = xx), interval = "confidence") plot(y ~ x, xlim = xr, ylim = yr) matlines(xx, cb, lty = c(1, 2, 2), col = 1) annotations(m0) }) return(panel) } maxd <- length(unique(x)) - 1 panel <- rp.control(x = x, y = y, er = er) rp.doublebutton(panel, variable = degree, action = draw.poly, showvalue = TRUE, step = 1, initval = 1, range = c(1, maxd), title = "Degree of polynomial") rp.do(panel, action = draw.poly) } with(cars, rp.poly(x = speed, y = dist, er = 0.2)) with(faithful, rp.poly(x = eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.poly(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- # Regression splines. rp.spline <- function(x, y, er = 0, ...) { plotfit <- function(x, y, m0, er) { xr <- extendrange(x, f = er) xx <- seq(xr[1], xr[2], length.out = 201) cb <- predict(m0, newdata = data.frame(x = xx), interval = "confidence") matlines(xx, cb, lty = c(1, 2, 2), col = 1) } annotations <- function(m0) { mtext(side = 3, adj = 0, line = 1.5, text = sprintf("X rank: %i", m0$rank)) mtext(side = 3, adj = 0, line = 0.5, text = sprintf("Residual DF: %i", m0$df.residual)) press <- sum((residuals(m0)/(1 - hatvalues(m0)))^2) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("PRESS: %0.2f", press)) sm0 <- summary(m0) mtext(side = 3, adj = 1, line = 1.5, text = sprintf("R² (adj. R²): %0.2f (%0.2f)", 100 * sm0$r.squared, 100 * sm0$adj.r.squared)) } draw.df.degree <- function(panel) { with(panel, { xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) plot(y ~ x, xlim = xr, ylim = yr, ...) nk <- as.integer(df1) - as.integer(degree1) nk <- max(c(0, nk)) if (quantile) { s <- seq(0, 1, length.out = 2 + nk) q <- quantile(x, probs = s[-c(1, nk + 2)]) abline(v = q, col = "gray50", lty = 2) } switch(base1, bs = { m0 <- lm(y ~ bs(x, degree = as.integer(degree1), df = as.integer(df1))) plotfit(x, y, m0, er) annotations(m0) mtext(side = 3, line = 0.5, text = sprintf( "Number of internal nodes: %i", nk)) }, ns = { m0 <- lm(y ~ ns(x, df = as.integer(df1))) plotfit(x, y, m0, er) annotations(m0) mtext(side = 3, line = 0.5, text = sprintf( "Number of internal nodes: %i", nk)) }) }) return(panel) } draw.degree.knots <- function(panel) { with(panel, { xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) plot(y ~ x, xlim = xr, ylim = yr, ...) abline(v = k, col = 2, lty = 2) switch(base2, bs = { if (is.na(k)) { m0 <- lm(y ~ bs(x, degree = as.integer(degree2))) } else { m0 <- lm(y ~ bs(x, degree = as.integer(degree2), knots = k)) } plotfit(x, y, m0, er) annotations(m0) }, ns = { if (is.na(k)) { m0 <- lm(y ~ ns(x, df = as.integer(degree2))) } else { m0 <- lm(y ~ ns(x, knots = k)) } plotfit(x, y, m0, er) annotations(m0) }) }) return(panel) } reset.draw <- function(panel) { k <- locator(type = "p", pch = 19, col = 2)$x if (length(k) == 0) { panel$k <- NULL #NA } else { panel$k <- k abline(v = k, col = 2, lty = 2) } panel$degree2 <- 1 rp.do(panel, draw.degree.knots) return(panel) } panel <- rp.control(x = x, y = y, k = NA, ...) rp.notebook(panel, tabnames = c("DF", "Knots"), tabs = c("Control DF", "Choose knots"), width = 280, height = 200) rp.radiogroup(panel, variable = base1, vals = c("bs", "ns"), action = draw.df.degree, title = "Base spline", parentname = "DF") rp.checkbox(panel, variable = quantile, action = draw.df.degree, title = "Show quantiles of x", parentname = "DF") rp.doublebutton(panel, variable = df1, action = draw.df.degree, showvalue = TRUE, step = 1, initval = 3, range = c(1, 10), title = "Degrees of freedom", parentname = "DF") rp.doublebutton(panel, variable = degree1, action = draw.df.degree, showvalue = TRUE, step = 1, initval = 3, range = c(1, 10), title = "Polynomial degree", parentname = "DF") rp.radiogroup(panel, variable = base2, vals = c("bs", "ns"), action = draw.degree.knots, title = "Base spline", parentname = "Knots") rp.doublebutton(panel, variable = degree2, action = draw.degree.knots, showvalue = TRUE, step = 1, initval = 3, range = c(1, 10), title = "Polynomial degree", parentname = "Knots") rp.button(panel, action = reset.draw, title = "Click to select knots", parentname = "Knots") rp.do(panel, action = draw.df.degree) } with(cars, rp.spline(x = speed, y = dist, er = 0.2, xlab = "Velocidade", ylab = "Distância")) with(faithful, rp.spline(x = eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.spline(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- rp.base <- function(x) { draw <- function(panel) { with(panel, { switch(base, bs = { X <- model.matrix(~bs(x, df = df, degree = degree)) matplot(x = x, y = X[, -1], type = "l") lines(x = x, y = rowSums(X[, -1]), col = 1, lty = 1, lwd = 2) }, ns = { X <- model.matrix(~ns(x, df = df)) matplot(x = x, y = X[, -1], type = "l") # lines(x = x, y = rowSums(X[, -1]), # col = 1, lty = 1, lwd = 2) }) }) return(panel) } panel <- rp.control(x = x) rp.radiogroup(panel, variable = base, vals = c("bs", "ns"), action = draw, title = "Base spline") rp.doublebutton(panel, variable = df, action = draw, showvalue = TRUE, step = 1, initval = 3, range = c(1, 10), title = "Degrees of freedom") rp.doublebutton(panel, variable = degree, action = draw, showvalue = TRUE, step = 1, initval = 3, range = c(1, 10), title = "Polynomial degree") } x <- seq(0, 10, 0.1) rp.base(x) #----------------------------------------------------------------------- # Smoothing splines. rp.smooth.spline <- function(x, y, er = 0, ...) { plotfit <- function(x, y, m0, er, ...) { xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) xx <- seq(xr[1], xr[2], length.out = 201) plot(y ~ x, xlim = xr, ylim = yr, ...) cb <- predict(m0, x = xx) lines(cb$x, cb$y, lty = 1, col = 1) } annotations <- function(m0, y) { mtext(side = 3, adj = 0, line = 1.5, text = sprintf("Equiv. DF: %0.1f", m0$df)) mtext(side = 3, adj = 0, line = 0.5, text = sprintf("Minimized crit.: %0.3f", m0$crit)) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("Smoothing parameter: %0.2f", m0$spar)) r2 <- 100 * cor(fitted(m0), y)^2 mtext(side = 3, adj = 1, line = 1.5, text = sprintf("R²: %0.2f", r2)) } draw.sspline <- function(panel) { with(panel, { switch(control, df = { m0 <- smooth.spline(x = x, y = y, df = df) }, spar = { m0 <- smooth.spline(x = x, y = y, spar = spar) }) plotfit(x, y, m0, er, ...) annotations(m0, y) }) return(panel) } draw.default <- function(panel) { m0 <- smooth.spline(x = x, y = y) plotfit(x, y, m0, er, ...) annotations(m0, y) mtext(side = 3, line = -1.5, text = "Default fit", col = 4) return(panel) } panel <- rp.control(x = x, y = y, er = er) rp.do(panel, action = draw.default) rp.radiogroup(panel, variable = control, vals = c("spar", "df"), labels = c("Smoothing parameter", "Equivalent degrees of freedom"), title = "Argument in control", action = draw.sspline) rp.slider(panel, variable = spar, from = 0, to = 1, initval = 0.5, action = draw.sspline, title = "Smoothing parameter", resolution = 0.05, showvalue = TRUE) rp.doublebutton(panel, variable = df, range = c(1.5, length(x) - 1), initval = 1.5, step = 0.5, showvalue = TRUE, action = draw.sspline, title = "Equivalent DF") rp.button(panel, action = draw.default, title = "Default fit") } with(cars, rp.smooth.spline(speed, dist, er = 0.2, xlab = "Velocidade", ylab = "Comprimento")) with(faithful, rp.smooth.spline(eruptions, waiting, er = 0.2)) with(MASS::mcycle, rp.smooth.spline(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- #----------------------------------------------------------------------- # Local polynomial regression. rp.loess <- function(y, x, n = 20, er = 0, ...) { annotations <- function(m0, y) { mtext(side = 3, adj = 0, line = 1.5, text = sprintf("H trace: %0.2f", m0$trace.hat)) mtext(side = 3, adj = 0, line = 0.5, text = sprintf("Equivalent num. param.: %0.2f", m0$enp)) r2 <- 100 * cor(fitted(m0), y)^2 mtext(side = 3, adj = 1, line = 1.5, text = sprintf("R²: %0.2f", r2)) } draw.loess <- function(panel) { with(panel, { m0 <- loess(y ~ x, span = span, degree = degree, family = "gaussian") erx <- extendrange(x, f = er) xx <- seq(erx[1], erx[2], length.out = 201) yy <- predict(m0, newdata = data.frame(x = xx)) a <- abs(x - x0) if (span < 1) { q <- as.integer(span * length(x)) d <- sort(a)[q] } else { q <- length(x) d <- max(abs(x - x0)) * sqrt(span) } w <- rep(0, length(x)) s <- a <= d w[s] <- (1 - (a[s]/d)^3)^3 i <- as.integer(s) xl <- range(x) xl[1] <- ifelse(x0 - d > xl[1], x0 - d, xl[1]) xl[2] <- ifelse(x0 + d < xl[2], x0 + d, xl[2]) f0 <- function(...) { y.pred <- sum(w * y)/sum(w) segments(xl[1], y.pred, xl[2], y.pred, col = 2, lwd = 2, ...) } f1 <- function(...) { m <- lm(y ~ poly(x, degree = 1), weights = w) y.pred <- predict(m, newdata = list(x = xl)) segments(xl[1], y.pred[1], xl[2], y.pred[2], col = 2, lwd = 2, ...) } f2 <- function(...) { m <- lm(y ~ poly(x, degree = as.integer(degree)), weights = w) x.pred <- seq(xl[1], xl[2], length.out = n) y.pred <- predict(m, newdata = list(x = x.pred)) lines(x.pred, y.pred, col = 2, lwd = 2, ...) } plot(x, y, pch = 2 * (!s) + 1, cex = i * 3 * w + 1, xlim = extendrange(x, f = er), ylim = extendrange(y, f = er), ...) lines(yy ~ xx) abline(v = c(x0, xl), lty = c(2, 3, 3)) annotations(m0, y) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("Number of obs. used/total: %i/%i", sum(w > 0), length(y))) switch(findInterval(degree, c(-Inf, 1, 2, Inf)), `0` = f0(), `1` = f1(), `2` = f2()) }) panel } xr <- extendrange(x, f = er) panel <- rp.control(x = x, y = y, n = n, er = er) rp.doublebutton(panel, variable = degree, action = draw.loess, showvalue = TRUE, step = 1, initval = 0, range = c(0, 2), title = "Degree") rp.slider(panel, variable = x0, action = draw.loess, from = xr[1], to = xr[2], initval = median(range(x)), showvalue = TRUE, title = "x-value") rp.slider(panel, variable = span, action = draw.loess, from = 0, to = 1.5, initval = 0.75, showvalue = TRUE, title = "span") rp.do(panel, action = draw.loess) } with(cars, rp.loess(x = speed, y = dist, n = 25, er = 0.2)) with(faithful, rp.loess(x = eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.loess(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- # GAM. rp.gam <- function(y, x, er = 0, ...) { annotations <- function(m0) { sm0 <- summary(m0) mtext(side = 3, adj = 0, line = 1.5, text = sprintf("Equivalent DF: %0.2f", sm0$edf)) mtext(side = 3, adj = 0, line = 0.5, text = sprintf("Residual DF: %0.1f", sm0$residual.df)) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("Smoothing parameter: %0.2f", m0$smooth[[1]]$sp)) mtext(side = 3, adj = 1, line = 1.5, text = sprintf("R²: %0.2f", 100 * sm0$r.sq)) } plotfit <- function(x, y, m0, er, ...) { b0 <- coef(m0)[1] xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) plot(m0, shift = b0, ylim = yr, xlim = xr, ...) points(y = y, x = x) } draw.gam <- function(panel) { with(panel, { m0 <- mgcv::gam(y ~ s(x, k = k, sp = sp)) plotfit(x, y, m0, er, ...) annotations(m0) }) return(panel) } draw.default <- function(panel) { with(panel, { m0 <- mgcv::gam(y ~ s(x)) plotfit(x, y, m0, er, ...) annotations(m0) mtext(side = 3, line = -1.5, text = "Default fit", col = 4) }) return(panel) } panel <- rp.control(x = x, y = y, er = er, ...) # panel <- rp.control(x = x, y = y, er = er)#, ...) maxd <- length(unique(x)) - 1 rp.doublebutton(panel, variable = k, action = draw.gam, showvalue = TRUE, step = 1, initval = -1, range = c(-1, maxd), title = "Degree of polynomial") rp.slider(panel, variable = sp, from = 0, to = 1, initval = 0.1, action = draw.gam, title = "Smoothing parameter", showvalue = TRUE) rp.button(panel, action = draw.default, title = "Default fit") rp.do(panel, action = draw.default) } with(cars, rp.gam(x = speed, y = dist, er = 0.2, xlab = "Velocidade", ylab = "Comprimento")) with(faithful, rp.gam(eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.gam(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- # K-vizinhos mais próximos (KNN regression), via caret::knnreg(). rp.knn <- function(y, x, er = 0, n.points = 200, ...) { plotfit <- function(x, y, xx, yy, er, ...) { xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) plot(y ~ x, xlim = xr, ylim = yr, ...) lines(xx, yy, lwd = 2, col = 1) } annotations <- function(fit, k, x, y) { yhat <- predict(fit, newdata = data.frame(x = x)) r2 <- 100 * cor(y, yhat)^2 mtext(side = 3, adj = 0, line = 0.5, text = sprintf("R²: %0.2f", r2)) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("k: %i", k)) } draw.knn <- function(panel) { with(panel, { fit <- caret::knnreg(x = data.frame(x = x), y = y, k = k) xr <- extendrange(x, f = er) xx <- seq(xr[1], xr[2], length.out = n.points) yy <- predict(fit, newdata = data.frame(x = xx)) plotfit(x, y, xx, yy, er, ...) annotations(fit, k, x, y) }) return(panel) } draw.default <- function(panel) { with(panel, { k0 <- max(1, round(length(x) / 10)) fit <- caret::knnreg(x = data.frame(x = x), y = y, k = k0) xr <- extendrange(x, f = er) xx <- seq(xr[1], xr[2], length.out = n.points) yy <- predict(fit, newdata = data.frame(x = xx)) plotfit(x, y, xx, yy, er, ...) annotations(fit, k0, x, y) mtext(side = 3, line = -1.5, text = "Default fit (k = n/10)", col = 4) }) return(panel) } panel <- rp.control(x = x, y = y, er = er, n.points = n.points) rp.do(panel, action = draw.default) rp.doublebutton(panel, variable = k, action = draw.knn, showvalue = TRUE, step = 1, initval = max(1, round(length(x) / 10)), range = c(1, length(x)), title = "Número de vizinhos (k)") rp.button(panel, action = draw.default, title = "Default fit") } with(cars, rp.knn(x = speed, y = dist, er = 0.2, xlab = "Velocidade", ylab = "Comprimento")) with(faithful, rp.knn(x = eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.knn(x = times, y = accel, er = 0.2)) #----------------------------------------------------------------------- # Kernel regression (Nadaraya-Watson), via ksmooth(). rp.ksmooth <- function(y, x, er = 0, n.points = 200, ...) { plotfit <- function(x, y, fit, er, ...) { xr <- extendrange(x, f = er) yr <- extendrange(y, f = er) plot(y ~ x, xlim = xr, ylim = yr, ...) lines(fit$x, fit$y, lwd = 2, col = 1) } annotations <- function(fit, x, y) { yhat <- approx(fit$x, fit$y, xout = x)$y r2 <- 100 * cor(y, yhat, use = "complete.obs")^2 mtext(side = 3, adj = 0, line = 0.5, text = sprintf("R²: %0.2f", r2)) } draw.ksmooth <- function(panel) { with(panel, { fit <- ksmooth(x, y, kernel = kernel, bandwidth = bandwidth, n.points = n.points) plotfit(x, y, fit, er, ...) annotations(fit, x, y) mtext(side = 3, adj = 1, line = 0.5, text = sprintf("Bandwidth: %0.2f", bandwidth)) }) return(panel) } draw.default <- function(panel) { with(panel, { fit <- ksmooth(x, y, kernel = "normal", bandwidth = 0.5, n.points = n.points) plotfit(x, y, fit, er, ...) annotations(fit, x, y) mtext(side = 3, line = -1.5, text = "Default fit", col = 4) }) return(panel) } bw.max <- diff(range(x)) panel <- rp.control(x = x, y = y, er = er, n.points = n.points) rp.do(panel, action = draw.default) rp.radiogroup(panel, variable = kernel, vals = c("normal", "box"), title = "Kernel", action = draw.ksmooth) rp.slider(panel, variable = bandwidth, from = bw.max/100, to = bw.max, initval = bw.max/10, resolution = bw.max/100, action = draw.ksmooth, title = "Bandwidth", showvalue = TRUE) rp.button(panel, action = draw.default, title = "Default fit") } with(cars, rp.ksmooth(x = speed, y = dist, er = 0.2, xlab = "Velocidade", ylab = "Comprimento")) with(faithful, rp.ksmooth(x = eruptions, y = waiting, er = 0.2)) with(MASS::mcycle, rp.ksmooth(x = times, y = accel, er = 0.2))