# code based on files "Code regarding section 4.2", "Code regarding section 4.3-4.4",
# "Code regarding section 4.5" and "Code regarding section 4.6"
# (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa)

# call the library splines for fitting splines
library(splines)

# download the data
dataset <- read.table("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", header = TRUE)
attach(dataset)

### code to reproduce figure 4.5a of Azzalini & Scarpa (2012) ###
# plot the data
xi <- seq(min(x), max(x), length = 4)
plot(x, y.yesterday, pch = 17, col = 2, ylab = "y", cex = 2)
abline(v = xi[2], lty = 3, lwd = 2)
abline(v = xi[3], lty = 3, lwd = 2)

# continuous piecewise linear fit
x0 <- seq(min(x), max(x), length = 200)
xx <- sort(c(x0, xi[2:3]))
m1 <- lm(y.yesterday ~  bs(x, knots = xi[2:3], degree = 1)) # degree 1 = linear
fit1<- predict(m1,  data.frame(x = xx))
lines(xx, fit1, col = 3, lwd = 2)

# continuous piecewise quadratic fit
m2 <- lm(y.yesterday ~  bs(x, knots = xi[2:3], degree = 2)) # degree 2 = quadratic
fit2 <- predict(m2,  data.frame(x = xx))
lines(xx, fit2, col = 4, lwd = 2)

# continuous piecewise cubic fit
m3 <- lm(y.yesterday ~  bs(x, knots = xi[2:3], degree = 3))
fit3 <- predict(m3,  data.frame(x = xx))
lines(xx, fit3, col = 5, lwd = 2)

# use this instead of legend
text(x = c(1.05, 0.85, 0.65), y = rep(0.52, 3), as.character(1:3))

# figure 4.5b
plot(c(0.5, 3), c(0, 2), type = "n", xlab = "x",
     ylab = "max(0,x-1.333)")
abline(v = xi[2], lty = 3)
lines(x0, pmax(0, x0 - xi[2]), col = 1, lwd = 2)

detach(dataset)

# example of thin-plate splines
auto <- read.table("http://azzalini.stat.unipd.it/Book-DM/auto.dat", header = TRUE)
attach(auto)
y <- city.distance
x <- engine.size
xlab <- "Engine size (L)"
ylab <- "Curb weight (kg)"
zlab <- "City distance (km/L)"

# figure 4.6
a075 <- smooth.spline(x, y, spar = 0.75)
a090 <- smooth.spline(x, y, spar = 0.95)
a125 <- smooth.spline(x, y, spar = 1.25)
x0 <- seq(min(x), max(x), length = 250)

par(mfrow = c(1, 1))
plot(x, y, xlab = xlab, ylab = zlab,  cex = 2, pch = 16)
lines(predict(a075, x = x0), col = 2, lty = 1, lwd = 3)
lines(predict(a090, x = x0), col = 3, lty = 1, lwd = 3)
lines(predict(a125, x = x0), col = 4, lty = 1, lwd = 3)
legend('topright', legend = c('l  =  0.0007', 
                              'l  =  0.02',
                              'l  =  2.9'),
       col = c(2, 3, 4), lty = c(1, 1, 1))


# code to plot figure 4.7
library(mgcv)

convex.hull <- function(x, eval.points)
{
  ngrid <- nrow(eval.points)
  hull.points <- x[order(x[, 1], x[, 2]), ]
  dh <- diff(hull.points)
  hull.points <- hull.points[c(TRUE, !((dh[, 1] == 0) & (dh[, 2] == 0))), ]
  hull.points <- hull.points[chull(hull.points), ]
  nh <- nrow(hull.points)
  gstep <- matrix(rep(eval.points[2, ] - eval.points[1, ], nh),
                  ncol = 2, byrow = TRUE)
  hp.start <- matrix(rep(eval.points[1, ], nh), ncol = 2,
                     byrow = TRUE)
  hull.points <- hp.start + gstep * round((hull.points -
                                             hp.start)/gstep)
  hull.points <- hull.points[chull(hull.points), ]
  grid.points <- cbind(rep(eval.points[ , 1], ngrid),
                       rep(eval.points[ , 2], rep(ngrid, ngrid)))
  D <- diff(rbind(hull.points, hull.points[1, ]))
  temp <- D[, 1]
  D[, 1] <- D[, 2]
  D[, 2] <- (-temp)
  C <- as.vector((hull.points * D) %*% rep(1, 2))
  C <- matrix(rep(C, ngrid^2), nrow = ngrid^2, byrow = TRUE)
  D <- t(D)
  wy <- ((grid.points %*% D) >= C)
  wy <- apply(wy, 1, all)
  wy[wy] <- 1
  wy[!wy] <- NA
  wy <- matrix(wy, ncol = ngrid)
  wy
}

m1 <- gam(city.distance ~ s(engine.size, curb.weight, k = 18), data = auto)

x1 <- seq(min(engine.size), max(engine.size), length = 50)
x2 <- seq(min(curb.weight), max(curb.weight), length = 50)
x12 <- data.frame(engine.size = as.vector(outer(x1, rep(1, length(x1)))),
                  curb.weight = as.vector(outer(rep(1, length(x2)), x2)))
m1p <- predict(m1, newdata = x12)

# figure 4.7a
matrice.spline <- matrix(m1p, 50, 50)
w1 <- convex.hull(cbind(engine.size, curb.weight), cbind(x1, x2))
persp(x1, x2, matrice.spline * w1, theta = 120, phi = 20, 
      xlab = "Engine size", ylab = "Curb weight", zlab = "City distance",
      ticktype = "detailed")

detach(package:mgcv, unload = TRUE)

y  <- city.distance
x1 <- engine.size
x2 <- curb.weight

library(gam)

a2 <- gam(city.distance ~ s(engine.size) + s(curb.weight),  data = auto)
pa2 <- preplot(a2)

# figure 4.12a
plot(pa2$"s(engine.size)", se = TRUE, ylim = c(-7, 10), lwd = 2,
     xlab = "Engine size", ylab = "s(Engine size)")

# figure 4.12b
plot(pa2$"s(curb.weight)", se = TRUE, ask = TRUE, ylim = c(-7,10), lwd = 2,
     xlab = "Curb weight", ylab = "s(Curb weight)")


library(sm)

x1 <- seq(min(engine.size), max(engine.size), length = 50)
x2 <- seq(min(curb.weight), max(curb.weight), length = 50)
x12<- data.frame(
  engine.size = as.vector(outer(x1, rep(1, length(x1)))),
  curb.weight = as.vector(outer(rep(1, length(x2)), x2)))
a2p <- predict(a2, x12)
sm.options(ngrid = 50)
h<- c(0.5, 150)

# figure 4.13a
matrice.gam <- matrix(a2p, 50, 50)
w1<-convex.hull(cbind(engine.size, curb.weight), cbind(x1, x2))
persp(x1, x2, matrice.gam * w1, theta = 120, phi = 20, 
      xlab = "Engine size", ylab = "Curb weight", zlab = "City distance",
      ticktype = "detailed")

# figure 4.13b
sm <- sm.regression(cbind(engine.size, curb.weight), city.distance, h = h,
                    options=list(xlab = xlab, ylab = ylab, 
                                 zlab = zlab), display = "none")
matrice.sm <- matrix(sm$estimate, 50, 50)
persp(x1, x2, matrice.sm, theta = 120, phi = 20, 
      xlab = "Engine size", ylab = "Curb weight", zlab = "City distance",
      ticktype = "detailed")