##########################
# Extra exercise 12      #
# STK2100, spring 2026   #
# Berit Omli Øksnes      #
# Based on zip.nn.r      #
# Last modified 16.04.26 #
##########################

# Clean workspace
rm(list=ls())

# Import packages
library(nnet)      # nnet
library(RSNNS)     # mlp
library(MASS)      # stepAIC

# Set seed for reproducibility
set.seed(2100)


# a) Read data ------------------------------------------------------------ ####
folder     <- "https://www.uio.no/studier/emner/matnat/math/STK2100/data/"
data.path  <- paste0(folder, "spam.data")
train.path <- paste0(folder, "spam.traintest")
spam       <- read.table(data.path)
train.ind  <- read.table(train.path)

n <- nrow(spam)
p <- ncol(spam) - 1
colnames(spam) <- c(paste0("x", 1:p), "cl")
spam$cl <- as.factor(spam$cl)

spam.train <- spam[train.ind==1,]
spam.test  <- spam[train.ind==0,]


# b) Single layer NN classifier with m=10 hidden variables ---------------- ####
# First normalize input variables using the training data
scaler <- scale(spam.train[,1:p])
spam.train[,1:p] <- scale(spam.train[,1:p],
                          center=attr(scaler, "scaled:center"),
                          scale=attr(scaler, "scaled:scale"))
spam.test[,1:p]  <- scale(spam.test[,1:p],
                          center=attr(scaler, "scaled:center"),
                          scale=attr(scaler, "scaled:scale"))

# Settings (trying a total of 3*3=9 models)
m.values     <- c(10, 5, 20)
decay.values <- c(0, 0.1, 0.5)
results.nnet <- data.frame(size=character(),
                           decay=numeric(),
                           error=numeric())

for (m in m.values) {
  for (decay in decay.values) {
    # Train neural network
    spam.nnet <- nnet(cl~., data=spam.train, size=m,
                      decay=decay, MaxNWts=10000, maxit=300)

    # Make predictions and compute error rate
    pred <- predict(spam.nnet, spam.test, type="class")
    err.nnet <- mean(spam.test$cl != pred)

    # Store results
    results.nnet[nrow(results.nnet)+1,] <- list(size=paste0("(", m, ")"),
                                                decay=decay,
                                                error=err.nnet)
  }
}

print(results.nnet, digits=3)
print(results.nnet[which.min(results.nnet[,"error"]),], digits=3)


# c) Train single hidden layer NN with mlp function ----------------------- ####
spamTargets.train <- decodeClassLabels(spam.train[,"cl"])

sizes       <- list(c(10), c(5), c(20))
results.mlp <- data.frame(size=character(),
                          error=numeric())

for (size in sizes) {
  # Train NN
  spam.dnet <- mlp(x = as.matrix(spam.train[,1:p]),
                   y = spamTargets.train,
                   size = c(size),
                   learnFunc="Std_Backpropagation",
                   learnFuncParams = c(0.3),
                   maxit = 100)

  # Make predictions and compute error
  pred.dnet <- 1*(predict(spam.dnet, newdata=spam.test[,1:p])[,2] >= 0.5)
  err.dnet <- mean(pred.dnet != spam.test$cl)

  # Store results
  results.mlp[nrow(results.mlp)+1,] <- list(
    size=paste0("(", paste0(size, collapse=","), ")"),
    error=err.dnet)
}

print(results.mlp, digits=3)
print(results.mlp[which.min(results.mlp[,"error"]),], digits=3)


# d) Multiple layer neural networks --------------------------------------- ####
sizes <- list(c(10,10),
              c(10,10,10),
              c(5,5,5,5))

for (size in sizes) {
  # Train NN
  spam.dnet <- mlp(x = as.matrix(spam.train[,1:p]),
                   y = spamTargets.train,
                   size = size,
                   learnFunc = "Std_Backpropagation",
                   learnFuncParams = c(0.3),
                   maxit = 100)

  # Make predictions and compute error
  pred.dnet <- 1*(predict(spam.dnet, newdata=spam.test[,1:p])[,2] >= 0.5)
  err.dnet <- mean(pred.dnet != spam.test$cl)

  # Store results
  results.mlp[nrow(results.mlp)+1,] <- list(
    size=paste0("(", paste0(size, collapse=","), ")"),
    error=err.dnet)
}

print(results.mlp, digits=3)
print(results.mlp[which.min(results.mlp[,"error"]),], digits=3)


# e) Logistic regression to perform model selection ----------------------- ####
# Fit logistic regression model using all input variables
spam.logistic <- glm(cl~., family="binomial", data=spam.train)

# Perform backward selection
bwd.selection <- stepAIC(spam.logistic, direction="backward")
subset <- attr(summary(bwd.selection)$terms, "term.labels")

# Fitting neural networks with the new subset as input variables
sizes <- list(c(10),
              c(10,10),
              c(10,10,10),
              c(5,5,5,5))

results.newsubset <- data.frame(size=character(),
                                error=numeric())

for (size in sizes) {
  # Train NN
  spam.dnet <- mlp(x = as.matrix(spam.train[,subset]),
                   y = spamTargets.train,
                   size = size,
                   learnFunc = "Std_Backpropagation",
                   learnFuncParams = c(0.3),
                   maxit = 100)

  # Make predictions and compute error
  pred.dnet <- 1*(predict(spam.dnet, newdata=spam.test[,subset])[,2] >= 0.5)
  err.dnet <- mean(pred.dnet != spam.test$cl)

  # Store results
  results.newsubset[nrow(results.newsubset)+1,] <- list(
    size=paste0("(", paste0(size, collapse=","), ")"),
    error=err.dnet)
}

print(results.newsubset, digits=3)
print(results.newsubset[which.min(results.newsubset[,"error"]),], digits=3)


# f) Summarize results ---------------------------------------------------- ####
{
  cat("Best single hidden layer NN: \n")
  print(results.nnet[which.min(results.nnet[,"error"]),],
        digits=3, row.names = FALSE)
  cat("Best deep NN: \n")
  print(results.mlp[which.min(results.mlp[,"error"]),],
        digits=3, row.names = FALSE)
  cat("Best deep NN with subset selection: \n")
  print(results.newsubset[which.min(results.newsubset[,"error"]),],
        digits=3, row.names = FALSE)
  }