# 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 necessary libraries
library(splines)
library(gam)
library(sm)


# auxiliary function, just to make plots
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
}


# read the data
auto <- read.table("http://azzalini.stat.unipd.it/Book-DM/auto.dat", header = TRUE)
attach(auto)
xlab <- "Engine size (L)"
ylab <- "Curb weight (kg)"
zlab <- "City distance (km/L)"

# define the variable (y = response, x1 and x2 covariates)
y  <- city.distance
x1 <- engine.size
x2 <- curb.weight

# fit an additive model (gam without specifying "family" fits an additive model)
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)")


# order the values for display reasons
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")