as_tibble(GAGurine)

GAG.lm <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) + I(Age^5) + 
              I(Age^6) + I(Age^7) + I(Age^8), data=MASS::GAGurine)
anova(GAG.lm)

GAG.lm2 <- lm(GAG ~ Age + I(Age^2) + I(Age^3) + I(Age^4) + I(Age^5) + 
              I(Age^6), data=GAGurine)

mypred <- data.frame(Age = seq(0, 17, len = 200))
mypred$GAG <- predict(GAG.lm2, newdata = mypred)
ggplot(GAGurine, aes(x = Age, y = GAG)) +
  geom_point(alpha=0.7) +  # Scatterplot of points
  geom_line(data = mypred, aes(x = Age, y = GAG), color = "blue") +  # Polynomial regression line
  labs(title = "Degree 6 polynomial")

attach(GAGurine)
GAG_ns5  <- lm(GAG ~ splines::ns(Age, df = 5))
GAG_ns10 <- lm(GAG ~ splines::ns(Age, df = 10))
GAG_ns20 <- lm(GAG ~ splines::ns(Age, df = 20))
GAG_sms  <- smooth.spline(Age, GAG)
xx <- seq(min(Age), max(Age), length.out = 200)
detach(GAGurine)
predicted05 <- predict(GAG_ns5,  newdata = data.frame(Age=xx))
predicted10 <- predict(GAG_ns10, newdata = data.frame(Age=xx))
predicted20 <- predict(GAG_ns20, newdata = data.frame(Age=xx))
predicted_sms <- predict(GAG_sms,  newdata = data.frame(Age=xx))
predicted_sms <- rename(as_tibble(predicted_sms), Age = x, GAG = y)
predicted_sms$df <- "Smoothing"
mypred <- data.frame(Age = xx, GAG = predicted05, df = "5")
mypred <- rbind(mypred, data.frame(Age = xx, GAG = predicted10, df = "10"))
mypred <- rbind(mypred, data.frame(Age = xx, GAG = predicted20, df = "20"))
mypred <- rbind(mypred, predicted_sms)
ggplot(GAGurine) + 
  geom_point(aes(x=Age, y=GAG), alpha=0.5, size=0.5) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=df, linetype=df)) +
  labs(title = "Splines") 

attach(GAGurine)
mypred <- as.data.frame(loess.smooth(Age, GAG))
names(mypred) <- c("Age", "GAG")
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  # default: predict(loess(y ~ x))
  geom_smooth(se=FALSE, color="magenta") + 
  geom_line(data=mypred, aes(x=Age, y=GAG), color="blue") +
  labs(title = "loess") 

attach(GAGurine)
mypred1 <- as.data.frame(ksmooth(Age, GAG, "normal", bandwidth = 1))
names(mypred1) <- c("Age", "GAG")
mypred1$bandwidth <- "1"
mypred2 <- as.data.frame(ksmooth(Age, GAG, "normal", bandwidth = 5))
names(mypred2) <- c("Age", "GAG")
mypred2$bandwidth <- "5"
mypred <- rbind(mypred1, mypred2)
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=bandwidth)) +
  labs(title = "ksmooth") 

attach(GAGurine)
h <- KernSmooth::dpill(Age, GAG)  # bandwidth selection
mypred1 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 0, bandwidth = h))
names(mypred1) <- c("Age", "GAG")
mypred1$degree <- "0"
mypred2 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 1, bandwidth = h))
names(mypred2) <- c("Age", "GAG")
mypred2$degree <- "1"
mypred3 <- as.data.frame(KernSmooth::locpoly(Age, GAG, degree = 2, bandwidth = h))
names(mypred3) <- c("Age", "GAG")
mypred3$degree <- "2"
mypred <- rbind(mypred1, mypred2, mypred3)
detach(GAGurine)
ggplot(GAGurine, aes(x=Age, y=GAG)) + 
  geom_point(alpha=0.75) + 
  geom_line(data=mypred, aes(x=Age, y=GAG, color=degree)) +
  labs(title = "locpoly") 

rock.lm <-  lm(log(perm) ~ area + peri + shape, data = rock) 
rock.gam <- gam::gam(log(perm) ~ s(area) + s(peri) + s(shape), 
                     control = gam::gam.control(maxit = 50), 
                     data = rock)
summary(rock.gam)

anova(rock.lm,  rock.gam) 

rock.gam1 <- gam::gam(log(perm) ~ area + peri + s(shape), data = rock) 
anova(rock.lm,  rock.gam1,  rock.gam)

par(mfrow =  c(2, 3), pty = "s") 
plot(rock.gam,  se = T)
plot(rock.gam1,  se = T)

mgcv::gam(log(perm) ~ s(area) + s(peri) + s(shape), data=rock)

head(MASS::cpus)

set.seed(123)
cpus2 <- select(cpus, -c(name, estperf)) # excludes names, authors' predictions
cpus.samp <- sample(seq_len(nrow(cpus2)), 100)
cpus.lm <- lm(log10(perf)~., data=cpus2 |> 
                              dplyr::slice(cpus.samp) )
cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
test.cpus <- function(fit)
  sqrt(sum((log10(cpus2 |> slice(-cpus.samp) |> select(perf)) -
            predict(fit, cpus2 |> slice(-cpus.samp)))^2)/109)
test.cpus(cpus.lm)

cpus.lm2 <- stepAIC(cpus.lm, trace=FALSE)
cpus.lm2$anova

test.cpus(cpus.lm2)

Xin <- as.matrix(cpus2 |> dplyr::slice(cpus.samp) |> select(-perf) )
cpus.mars <- mda::mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) )) 
test2 <- function(fit) {  # test error
    Xp <- as.matrix(cpus2 |> slice(-cpus.samp) |> select(-perf) ) #cpus2[-cpus.samp,  1:6])
    sqrt(sum((log10(cpus2 |> slice(-cpus.samp) |> select(perf) ) - predict(fit, Xp))^2)/109)
}
test2(cpus.mars) 

showcuts <- function(obj) {
    tmp  <-  obj$cuts[obj$sel, ]
    dimnames(tmp)  <- list(NULL,  colnames(Xin)) 
    tmp
}
showcuts(cpus.mars) 

# examine the fitted functions
Xp <- matrix(cpus2 |> slice(cpus.samp) |> select(-perf) |> map_dbl(mean), 
             100, 6, byrow = T)
xr <- sapply(cpus2, range) # variable ranges
df_plt <- list()
for (i in seq_len(ncol(Xp))) {
    Xp1 <- Xp  # hold other variables
    xvals <- seq(xr[1, i], xr[2, i], len = 100)
    Xp1[, i] <- xvals
    df_plt[[i]] <- data.frame(x = xvals, y = predict(cpus.mars, Xp1), 
                              variable = names(cpus2)[i])
}
df_long <- dplyr::bind_rows(df_plt)
ggplot(df_long, aes(x = x, y = y, color = variable)) + geom_line() +
  facet_wrap(~ variable, ncol = 3, scales = "free")

## allow pairwise interaction terms
cpus.mars2 <- mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) ), degree = 2) 
showcuts(cpus.mars2) 

test2(cpus.mars2) 

## allow arbitrary interactions
cpus.mars6 <- mars(Xin, log10(cpus2 |> slice(cpus.samp) |> select(perf) ), degree = 6) 
showcuts(cpus.mars6) 

test2(cpus.mars6) 

summary(rock.ppr) 

par(mfrow = c(3,  2))
plot(rock.ppr)
plot(update(rock.ppr,  bass = 5)) # supsmu control
# use smoothing spline with GCV
plot(rock.ppr2 <- update(rock.ppr,  sm.method = "gcvspline", gcvpen = 2))

summary(rock.ppr2) 

summary(rock1) # to find the ranges of the variables

Xp  <- expand.grid(area = seq(0.1, 1.2, 0.05), peri = seq(0, 0.5, 0.02), 
                   shape = 0.2)
rock.grid <- cbind(Xp, fit = predict(rock.ppr2, Xp)) 
lattice::wireframe(fit ~ area +  peri, rock.grid,
                   screen = list(z = 160,  x = -60),  
                   aspect = c(1,  0.5), drape = T)

(cpus.ppr <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                 nterms = 2, max.terms = 10, bass = 5)) 

cpus.ppr <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                nterms = 8, max.terms = 10, bass = 5)  # select 8 terms
test.cpus(cpus.ppr) 

ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
    nterms = 2, max.terms = 10, sm.method = "spline")

cpus.ppr2 <- ppr(log10(perf)~., data = cpus2 |> dplyr::slice(cpus.samp),
                 nterms = 7,  max.terms = 10, sm.method = "spline")
test.cpus(cpus.ppr2) 

res3 <- unlist(log10(cpus2 |> dplyr::slice(-cpus.samp) |> select(perf)) - 
                  predict(cpus.ppr, cpus2 |> dplyr::slice(-cpus.samp)) ) 
res2 <- unlist(log10(cpus2 |> slice(-cpus.samp) |> select(perf)) -
                  predict(cpus.lm2, cpus2 |> slice(-cpus.samp)))  # lm
wilcox.test(res2^2, res3^2, paired = T, alternative = "greater") 

nnet(formula, data, weights, size, Wts, linout = F, entropy = F, 
    softmax = F, skip = F, rang = 0.7, decay = 0, maxit = 100,
    trace = T)

rock1 <- rock |> dplyr::mutate(area = area / 10000, peri = peri / 10000)
rock.nn <- nnet::nnet(log(perm) ~ area + peri + shape, data = rock1, 
                      size = 3, decay = 1e-3, linout = T,  skip = T,
                      maxit = 1000, Hess = T)

summary(rock.nn)

sum((log(rock$perm) - predict(rock.nn))^2)

eigen(rock.nn$Hessian, T)$values

attach(cpus2)
cpus3 <- data.frame(syct = syct - 2, mmin = mmin - 3, mmax = mmax - 4,
                    cach = cach / 256,  chmin = chmin / 100,
                    chmax = chmax / 100, perf = perf)
detach(cpus2)
test.cpus <-  function(fit) {
    sqrt(sum((log10(cpus3 |> slice(-cpus.samp) |> select(perf)) - 
              predict(fit, cpus3 |> slice(-cpus.samp)))^2)/109)
}
cpus.nn1 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip = T,  size = 0) # linear model

test.cpus(cpus.nn1)

cpus.nn2 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip=T, size=4, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn2) 

cpus.nn3 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T,
                 skip=T, size=10, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn3)

cpus.nn4 <- nnet(log10(perf)~., cpus3 |> slice(cpus.samp), linout = T, 
                 skip=T, size=25, decay=0.01, maxit = 1000, trace=F)
test.cpus(cpus.nn4)  

CVnn.cpus <- function(formula, data = cpus3,
    size = c(0, 4, 4, 10, 10),
    lambda = c(0, rep(c(0.003, 0.01), 2)),
    nreps = 5, nifold = 10, ...) {
  CVnn1 <- function(formula, data, nreps=1, ri,  ...) {
    truth <- log10(data$perf)
    res <- numeric(length(truth))
    cat("  fold")
    for (i in sort(unique(ri))) {
      cat(" ", i,  sep="")
      for(rep in 1:nreps) {
        learn <- nnet(formula, data[ri !=i,], trace=FALSE, ...)
        res[ri == i] <- res[ri == i] +
                        predict(learn, data[ri == i,])
      }
    }
    cat("\n")
    sum((truth - res/nreps)^2)
  }
  choice <- numeric(length(lambda))
  ri <- sample(nifold, nrow(data), replace = TRUE)
  for(j in seq(along=lambda)) {
    cat("  size =", size[j], "decay =", lambda[j], "\n")
    choice[j] <- CVnn1(formula, data, nreps=nreps, ri=ri,
                       size=size[j], decay=lambda[j], ...)
    }
  cbind(size=size, decay=lambda, fit=sqrt(choice/100))
}
CVnn.cpus(log10(perf) ~ ., data = cpus3 |> slice(cpus.samp),
          linout = TRUE, skip = TRUE, maxit = 1000)