hw1_markdown.Rmd

---
title: "Statistical Learning, Homework-01"
author: "Barba Paolo, Candi Matteo, Costantini Silvia, Vestini Maria Vittoria"
date: '2023-05-10'
output: 
  html_document:
    code_folding: hide
    theme: 
      color-contrast-warnings: false
      bg: "#2B3E50"
      fg: "#B8BCC2"
      primary: "#EA80FC"
      secondary: "#00DAC6"
      base_font:
      google: Prompt
      heading_font:
        google: Proza Libre
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
if (requireNamespace("thematic")) 
  thematic::thematic_rmd(font = "auto")
```

```{r package, message=FALSE, warning=FALSE}

rm(list = ls())

# Set reproducibility -----------------------------------------------------
seed <- 1234
set.seed(seed) 

# Libraries and Data ---------------------------------------------------------------

pkt <- c('glmnet', 'NMOF', 'parallel', 'snow', 'ggplot2', 'reshape2', 'gridExtra','tidyverse')

invisible(lapply(pkt, require, character.only = TRUE))

```

## Purpose & Statistical tools used

The homework's goal is to try to estimate an unknown function
$f(X)$ in order to predict the response variable $Y$ from
independent data $(X_1,Y_1) \dots ,(X_n,Y_n)$, to achieve this goal we will assume:

- $Y_i | X_1 , \dots , X_n \sim p(y | f(X_i))$

- $f$ lies in some fixed class of
functions.


For the Homework's sake we are going to use a penalized linear regression model and truncated power basis function, obtaining
the following: $$ Y_i = f(\underline{x_i}) + \epsilon_i \ for \ i = 1\dots n $$
where the $\epsilon_i$ are $i.i.d$ data from some normal distribution
$N(0,\sigma^2)$.

In order to get the optimal hyperparameters and MSE estimation we will
implement CV and nested CV techniques.

## Be linear in transformed feature space

The basic idea of being linear in a transformed feature space is to map
our parameters into a higher dimensional space with d-dimensions. The
model will then be linear in that space.

Let $\phi(x)$ be a feature mapping function that maps the input data $x$
into a transformed d-dimensional space. Let $w$ be a weight vector in
the transformed space. Thus, the prediction function is given by
$f(x) = w^T\phi(x)$. By mapping the parameters into higher dimensional
space, we can capture more complex structures between the features and
the target variable.

The polynomial models in a d-dimensional feature space assume that there
is a fixed functional form for the relationship between features and the
response variable. And since this relationship is defined globally,
these models may not capture the entire complexity of the data.

Spline regression models, instead, provide more flexible and adaptive
modeling in order to catch up the local structures of the target
variable $Y$.

In order to implement splines, we use truncated power functions
$\mathcal{G}_{d,q}$, defined as:

$$\mathcal{G}_{d,q} = \{g_1(x) \dots g_{d+1}(x),g_{(d+1)+1}(x),\dots,g_{(d+1)+q}(x) \}$$
where:

$\{g_1(x)=1,g_2(x)=x,...,g_{d+1}(x)=x^{d}\}$ and
$\{g_{(d+1)+j}(x)= (x−\xi_{j})^d_+\}$ where $(x)_+ =max\{0,x\}$.

It can be easily shown that the spline function $f(x)$ can be
represented by a linear combination of the power functions
$\mathcal{G}_{d,q}$:

$$f(x)= \sum_{j= 1}^{d+q+1} β_j ·g_j(x)$$

where $\beta_j$ are the coefficients to be estimated, and $g_j(x)$
denotes the $j$-th power function of  $\mathcal{G}_{d,q}$.

## Implementation of the truncated power basis

We have developed our own function of the truncated power basis
$\mathcal{G}_{d,q}$. This function takes as input the degree $d$ of the
truncated power functions, the number of knots $q$, and their positions,
and returns a feature matrix $\mathbb{X}$. The feature matrix
$\mathbb{X}$ will be used to obtain a prediction of the response
variable $Y$.

In order to illustrate our implementation, in the following plots,
results are shown considering $d\ \in \{1 , 3\}$ and
$q\ \in \{3 , 10\}$.

```{r power_function}

# Function used to compute the feature matrix
power_functions <- function(d, q, knots, x){
  X <- matrix(NA, length(x), d+q+1)  # Pre-set the feature matrix 
  for( i in 1:length(x)){            # Loop over the data-points
    
    for( j in 1:(d+q+1) ){           # Loop over the basis used
      
      if ( j <= d+1 ){               # Check if the index belong to the firsts (d+1) indices  
        X[i,j] <- x[i]^(j-1)         # Compute the powers
      }
      
      else                           # The index do not belong to the firsts (d+1) indices  
        if((x[i] - knots[j-(d+1)])^d > 0){      # Check if the power of the difference between the data-point and the knot is positive
          X[i,j] <- (x[i] - knots[j-(d+1)])^d   # Compute the power of the difference
        }
      else                           # the power is not positive
         X[i,j] <- 0                 # Set the value to 0
    }   
  }
  return(X)                         # Return the feature matrix
}

```

```{r plot truncated power function, fig.showtext = TRUE, fig.height = 7, fig.width = 12, fig.align = "center", message=FALSE, warning=FALSE , echo = FALSE}

plot_func <- function(d, q){
  knots <- seq(0, 1, length.out = q+2)[2 : (q+1)]
  df <- data.frame(power_functions(d, q, knots, xx))
  n_col <- ncol(df)
  legend_labels <- c()
  for (i in 1:n_col) legend_labels <- c(legend_labels, paste0('g', i, '(x)'))
  colnames(df) <- legend_labels
  df <- cbind('xx' = xx, df)
  df <- melt(df,  id.vars = 'xx')
  
  graph <- ggplot(data=df, aes(x = xx, y = value, color = variable)) + 
    geom_line(linewidth = 1) + 
    ggtitle(paste0('Plot with d = ', d, ' and  q = ', q)) + 
    xlab('x') +
    ylab('gi(x)') +
    theme_light() + 
    theme(panel.grid = element_blank(),
          panel.background = element_rect(fill = NA),
          title = element_text(colour = "white"),
          plot.background = element_rect(fill = NA, colour = "#2B3E50"), 
          plot.title = element_text(size=16),
          axis.line = element_line(color = "white"),
          axis.text = element_text(color = "white", size = 10),
          legend.title = element_blank(),
          legend.text = element_text(size = 10, color = 'white'),
          legend.key = element_rect(color = NA, fill = NA),
          legend.key.size = unit(0.6, "cm"),
          legend.background = element_rect(fill = NA))
      
  return(graph)
}

xx <- seq(0, 1, length.out = 1000)

p1 <- plot_func(1, 3)
p2 <- plot_func(3, 3)
p3 <- plot_func(1, 10)
p4 <- plot_func(3, 10)

grid.arrange(p1, p2, p3, p4, ncol=2)

```

## Comments

For $j \le d +1$, the truncated power function is defined as
$g_j(x) = x^{j-1}$. For $j > d +1$, the truncated power function is
defined as $g_j(x) = (x-\xi_j)^d_+$, where $(x)_+ = \max\{0,x\}$, and
$\xi_j$ is the $(j-d-1)$-th knot.

For $j > d +1$, the truncated power function is equal to zero for data
points below the knot $\xi_j$, and $(x-\xi_j)^d$ for data points above
the knot.

The use of truncated power basis allows us to work with a higher
dimensional feature space without using a high degree polynomial. The
additional dimensions in the feature space come from the truncated power
functions that capture the local structure of the target function around
each knot. This can lead to more accurate modeling of complex
relationships between the input features and the target variable.

An issue that can arise in spline regression is the presence of a
non-singular matrix, which can lead to numerical instability in order to
obtain the parameter estimates.

This problem typically arises when the knots are not well-spaced or when
the degree of the spline is too high relative to the number of
observations. In these cases, the matrix of the spline basis functions
can become rank-deficient, which means that it has fewer linearly
independent rows or columns than its dimensions.

## ChatGPT - Power function

In this section we ask CHAT-GPT to implement its own version of
truncated power basis. Below are reported the conversation we had with
it.

**1° Request:**

" Can you implement a function in R to compute the truncated power
basis, with parameters $d$, $q$, knots. Where $g_i=x^{i-1}$ if $i\le(d+1)$
else $g_{d+1+j} = \max((x-knots[j])^d, 0)$ for $j \in [1,q]$? "

**1 ° Response:**

" Certainly! Here's a function in R that computes the truncated power
basis you described:

```{r GPT, class.source = "fold-show"}

truncated_power_basis <- function(x, d, q, knots) {
  
  # Compute the first d+1 basis functions
  basis <- matrix(0, nrow = length(x), ncol = d+1)
  basis[,1] <- 1
  for (i in 2:(d+1)) {
    basis[,i] <- x^(i-1)
  }
  
  # Compute the remaining basis functions using knots
  for (j in 1:q) {
    basis[,d+1+j] <- pmax((x - knots[j])^d, 0)
  }
  
  return(basis)
}
```

This function takes as input the vector $x$, the integer $d$, the integer $q$,
and the vector knots. It returns a matrix with length(x) rows and $d+1+q$
columns, where each column represents a basis function."

Let's give a try to the ChatBot implementation!

```{r GPtest, error=TRUE, class.source = "fold-show"}

d <- 1
q <- 3
knots <- seq(0, 1, length.out=q)
x <- seq(0, 1, length.out=1000)

M <- truncated_power_basis(x, d, q, knots)

```

Oopsie Doopsie! Looks like the ChatBot response doesn't work. So we will
try and see how much we have to tell it in order to obtain a working
function, the first thing we will do is show it the error we got.

*For this following part the code will not be directly shown.*

**2° message request:**

" It seems like there is an error with your function, if i try to use it
I get this error: "Error in `[<-`(`*tmp*`, , d + 1 + j, value =
pmax((x - knots[j])\^d, 0)): subscript out of bounds". Could you give me
the correct implementation? "

**2° response:**

" I apologize for the error in my previous implementation. Here's an
updated version of the truncated_power_basis function that should work
correctly:

```{r GPtest2, error=TRUE}

truncated_power_basis <- function(x, d, q, knots) {
  # Compute the first d+1 basis functions
  basis <- matrix(0, nrow = length(x), ncol = d+1)
  basis[,1] <- 1
  for (i in 2:(d+1)) {
    basis[,i] <- x^(i-1)
  }
  
  # Compute the remaining basis functions using knots
  if (q > 0) {
    for (j in 1:q) {
      basis[,d+j+1] <- pmax((x - knots[j])^d, 0)
    }
  }
  
  return(basis)
}

```

... "

Again it doesn't work, so this time we give the ChatBot a hint, since we
realize that the problem is that there isn't enough allocated space for
the matrix, in fact the matrix should have $q+d+1$ columns instead of
$d+1$.

**3° message:**

" It doesn't solve it, I think the problem is in the dimension of the
matrix. "

**3° response:**

" I apologize for the continued issue. You are correct that there is an
issue with the dimension of the matrix. Here is the corrected
implementation of the truncated_power_basis function:

```{r GPtest3, error=TRUE}

truncated_power_basis <- function(x, d, q, knots) {
  # Compute the first d+1 basis functions
  basis <- matrix(0, nrow = length(x), ncol = d+1)
  basis[,1] <- 1
  for (i in 2:(d+1)) {
    basis[,i] <- x^(i-1)
  }
  
  # Compute the remaining basis functions using knots
  if (q > 0) {
    for (j in 1:q) {
      basis[,d+j+1] <- pmax((x - knots[j])^(d*(x >= knots[j])), 0)
    }
  }
  
  return(basis)
}
```

... "

At the end since also this implementation doesn't work we decided to
tell the ChatBot how to solve the issue.

**4° message:**

" Still not good, shouldn't the matrix have $q+d+1$ columns since we are
trying to evaluate $q+d+1$ functions? "😅

**4° response:**

" You are correct. I apologize for the error. The matrix should have
$q + d + 1$ columns since we are computing a total of $q + d + 1$ basis
functions. Here is the corrected implementation:

```{r GPtest4, error=TRUE}
truncated_power_basis <- function(x, d, q, knots) {
  # Compute the first d+1 basis functions
  basis <- matrix(0, nrow = length(x), ncol = d+1)
  basis[,1] <- 1
  for (i in 2:(d+1)) {
    basis[,i] <- x^(i-1)
  }
  
  # Compute the remaining basis functions using knots
  if (q > 0) {
    for (j in 1:q) {
      basis[,d+j+1] <- pmax((x - knots[j])^(d*(x >= knots[j])), 0)
    }
  }
  
  return(basis)
}

```

... "

But again even if we straight up told the Bot what's the problem and how
to solve it the function is not working! At the end the only thing left
to do is just tell it the code it should use.

**5° (and final) message:**

" You are not allocating the memory for the matrix correctly, basis
should be defined as "basis \<- matrix(0, nrow = length(x), ncol =
d+q+1)". "

**5° response:**

" You are correct again. I apologize for the confusion. Here is the
corrected implementation:

```{r GPtest5, error=TRUE, class.source = "fold-show"}
truncated_power_basis <- function(x, d, q, knots) {
  # Compute the first d+1 basis functions
  basis <- matrix(0, nrow = length(x), ncol = d+q+1)
  basis[,1] <- 1
  for (i in 2:(d+1)) {
    basis[,i] <- x^(i-1)
  }
  
  # Compute the remaining basis functions using knots
  if (q > 0) {
    for (j in 1:q) {
      basis[,d+j+1] <- pmax((x - knots[j])^(d*(x >= knots[j])), 0)
    }
  }
  ?pmax
  return(basis)
}
```

$\dots$ "

Finally we obtain a function without errors!

```{r GPtest final, error=TRUE, class.source = "fold-show"}
truncated_power_basis(seq(0,1,length.out=6), 3, 2, c(0.1, 0.5))
```

Sadly though we realize that, from the third response on, the ChatBot
changed a line of code from $basis[,d+j+1] <- pmax((x - knots[j])^d, 0)$
to $basis[,d+j+1] <- pmax((x - knots[j])^{d*(x >= knots[j])}, 0)$ and
now the function indeed does gives an output, but it's not the right
one.

## Comments:

After the chat we had with ChatGPT, our thoughts are that the chatbot
can be very useful for providing a starting point. We believe, however,
that can be risky for non-experts in the field to use the chatbot, since
they may not be able to identify errors or misunderstandings in the
information that the chat provides.

We have to remark that chatbots have limitations in fields that require
specialized knowledge or experience, such as statistics or data
analysis.😉

## Dataset

The Dataset we are going to use to implement the model is about
**Wilkinson Microwave Anisotropy Probe**, also known as **WMAP**. An
explorer mission that measured the difference of temperature across the
space.

The Dataset is composed by $675$ observations of $2$ variable of
interest.

```{r dataset}

test_set_vero <- read.csv("test.csv")
train_set <- read.csv("train.csv")

```

```{r plot dataset, fig.showtext = TRUE, fig.height = 5, fig.width = 10, fig.align = "center", message=FALSE, warning=FALSE , echo = FALSE}

# Plot --------------------------------------------------------------------
colors <- c("Real Data" = "snow")
point_size <- 1.4

ggplot() +
  geom_point(aes(x = train_set$x, y = train_set$y, color = 'Real Data'), size = point_size, shape=16) +
  theme_light() +
  labs(x = "x", y = "y", color = "     Legend", title = 'Train set', shape = "", color="") +
  guides(color = guide_legend(override.aes=list(shape = 16, size = 2))) +
  scale_color_manual(values = colors) +
  theme(panel.border = element_rect(color = "white"),
        legend.title = element_text(size=12, color='white'),
        plot.title = element_text(hjust = 0.5,size=16),
        legend.text = element_text(size = 10, color = 'white'),
        legend.key = element_rect(color = NA, fill = NA),
        legend.background = element_rect(fill = NA),
        panel.grid.major = element_line(linetype = 2),
        panel.grid.minor = element_line(linetype = 0),
        panel.background = element_rect(fill = NA),
        title = element_text(colour = "white"),
        plot.background = element_rect(fill = NA, colour = "#2B3E50"), 
        axis.line = element_line(color = "white"),
        axis.text = element_text(color = "white", size = 10),
        axis.title.y = element_text(angle = 0, vjust = 0.5))

```

## Vanilla Cross Validation

In order to estimate the coefficients $\beta_j$ in our penalized linear
model, we are going to perform a *k-fold* cross validation.

Since the variance of the response variable is not constant across the
range of predictor variables it can be possible that, picking knots
where the functions is very variable, may influence the good predictive
performance of our model. In order to try to resolve this problem, we
fix the lower bound of the position of the knots at $0$ and we add the
hyperparameter $p$ that is representing the upper bound. The knots will
be equispaced in interval $(0,p]$.

Once we have performed the power functions, we train a penalized linear
model, where we tune the type of regularization from Ridge to Lasso
throughout some elastic net.

At the end of the day, we end up with $5$ hyperparameter to tune:

-   $\alpha$ for Ridge, Lasso and elastic net
-   $\lambda$ Regularization parameter
-   $d$ degree of the power functions
-   $q$ the number of knots
-   $p$ positions of the upper bound

Since the data shown that for a higher $x$, the $y$ are way more variable,
in the train of our model, we will use into the model weights
proportional to the variance given $x$.

We will down-weight observations with high variance and up-weight
observations with low variance, which can decrease the error of our
predictions.

Given the knots we estimated the conditional variance by computing the
empirical variance of $y$ in the intervals defined by the knots.

The idea is that after cross validation we shrink the space of our
parameters that we consider sub-optimal. Then we will perform the nested
cross validation in order to select in a more accurate way the
hyperparameters.

```{r functions VCV}

# Function to estimate the weights
compute_weights <- function(knots , dataset){
  # Number of knots
  n <- length(knots)
  # Add the bounds 
  knots <- c(0, knots)
  if(knots[n+1] != 1){
    knots <- c(knots, 1)
  }
  
  # Define the matrix
  xx <- rep(NA, length(dataset$x))
  # Intialize v in case there are no data-points in the first interval
  v <- 1
  
  # For each interval bounded by the knots
  for(i in 1:(n+1)){
    # Check if there are at least 2 data points inside the interval
    if(sum(dataset$x >= knots[i] & dataset$x <= knots[i+1])>1){    
      # Compute the conditional variance 
      v  <- var(dataset$y[(dataset$x >= knots[i]) & (dataset$x <= knots[i+1])]) 
      # Set the weight of the data points equal to the inverse of the variance
      xx[dataset$x >= knots[i] & dataset$x <= knots[i+1]] <- 1/v              
    } 
    # If there is only one data point
    else if(sum(dataset$x >= knots[i] & dataset$x <= knots[i+1]) == 1){  
      # Since we can't compute the variance of a single data point then we use the previously computed variance
      xx[dataset$x >= knots[i] & dataset$x <= knots[i+1]] <- 1/v
    }
  }
  
  return(xx)
}


# Function used for the cross validation
cross_val_func <- function(x){
  set.seed(070720)
  # Define the parameters
  d <- x[1]   
  q <- x[2]
  k <- x[3]
  a <- x[4]
  l <- x[5] 
  p <- x[6]
  
  # size of the fold
  l_folds <- nrow(train_set) / k 
  # Pre-set the scores (MSE)
  score <- rep(NA, k)
  # Random  sample of the indices
  idx <- sample((1:nrow(train_set)),nrow(train_set))
  #Loop over the folds
  for ( i in 1:k){
    # Set the validation and train set
    cv_test <- train_set[idx[((i-1)*l_folds+1): (i*l_folds)],]
    cv_train <- train_set[-idx[((i-1)*l_folds+1): (i*l_folds)],]
    
    # Conpute the knots
    knots <- seq(0, p, length.out=q+1)[2:(q+1)]
    
    # Compute the power function on the training set
    M_cv_train <- power_functions(d = d, q = q, knots = knots, x = cv_train$x)
    
    # Compute the power function on the valifation set
    M_cv_test <-  power_functions(d = d, q = q, knots = knots, x = cv_test$x)

    # compute the weights, based on the conditional variance
    hat_weights <- compute_weights(knots = knots , cv_train)
    
    # Train the model
    cv_model <- glmnet(M_cv_train, 
                       cv_train$y,
                       family = "gaussian", 
                       alpha=a, 
                       lambda=l,
                       weights = hat_weights)
    
    # Preditcion on the validation set
    cv_predictions <- predict(cv_model, M_cv_test)
    
    # Compute the score
    score[i] <- sqrt(mean((cv_test$y-cv_predictions)^2))
    
  }
  # Return the mean MSE  
  return(mean(score))
}
```

```{r cv vanilla, eval = FALSE}
# parameters
k <- c(5)
d_grid <- c(1, 3) 
q_grid <- seq(3, 50, 2)
positions <- seq(0.2, 0.8, 0.1)
lambdas <- 10^seq(-0.5, 0, .05)
alphas <- seq(0, 1, 0.1)
# Set the parameter for the CV
parameters <- list(d_grid, q_grid, k, alphas, lambdas, positions)

# CV vanilla --------------------------------------------------------------
# Select the best combination of parameters
cl = makeCluster(detectCores())
clusterExport(cl, c('train_set', 'power_functions', 'glmnet', 'compute_weights'))
res <- gridSearch(cross_val_func, levels=parameters, method = 'snow', cl=cl)
stopCluster(cl)
best_params <- res$minlevels
names(best_params) <- c('d', 'q', 'k', 'alpha', 'lambda', 'position')
```

After the vanilla cross validation this are the hyperparameters
we obtained:

```{r training with vanilla parameters, warning = FALSE }
load("RData//best_params_vanilla_NODC.RData")

# Prediction --------------------------------------------------------------

# Using the best parameters
d_best <- best_params[1]
q_best <- best_params[2]
k_best <- best_params[3]
a_best <- best_params[4]
l_best <- best_params[5]
p_best <- best_params[6]

# Compute the predictions
knots <- seq(0, p_best, length.out=q_best+1)[2:(q_best+1)]
M_train <- power_functions(d = d_best, q = q_best, knots = knots, x = train_set$x)
M_test <- power_functions( d = d_best, q = q_best, knots = knots, x = test_set_vero$x)
knots_test <- power_functions( d = d_best, q = q_best, knots = knots, x = knots)
hat_weights <- compute_weights(knots , train_set)
final_model <- glmnet(M_train, train_set$y, family ="gaussian", 
                      alpha=a_best, lambda=l_best , weights = hat_weights )
predictions <- predict(final_model,M_test)
round(best_params,3)
```

After training the model with those hyperparameters and predicting the
$Y$ in the test set, we can plot our results.

```{r plot after vanilla, fig.showtext = TRUE, fig.height = 5, fig.width = 10, fig.align = "center",echo = FALSE, message=FALSE, warning=FALSE}

colors <- c("Real Data" = "#F1F1E6", "Predicted" = "#E09F20", "Knots" = "#006054")
point_size <- 1.4
knots_size <- 1
knots_shape <- 18


ggplot() +
  geom_point(aes(x = train_set$x, y = train_set$y, color = 'Real Data'), size = point_size, shape=16) +
  geom_line(aes(x = test_set_vero$x, y = predictions, color = 'Predicted'), linewidth = 2) +
  geom_point(aes(x = knots, y = predict(final_model, knots_test), color = 'Knots'), shape = knots_shape, stroke = 1.7, size = knots_size) +
  theme_light() +
  labs(x = "x", y = "y", color = "     Legend", title = 'Prediction on WMAP data after cross validation', shape = "", color="") +
  scale_color_manual(values = colors) +
  guides(color = guide_legend(override.aes=list(shape = c(knots_shape, NA, 18), size = 2, linetype = c(NA, 1, NA)))) +
  theme(panel.border = element_rect(color = "white"),
        legend.title = element_text(size=12, color='white'),
        plot.title = element_text(hjust = 0.5,size=16),
        legend.text = element_text(size = 10, color = 'white'),
        legend.key = element_rect(color = NA, fill = NA),
        legend.background = element_rect(fill = NA),
        panel.grid.major = element_line(linetype = 2, linewidth=.2),
        panel.grid.minor = element_line(linetype = 0),
        panel.background = element_rect(fill = NA),
        title = element_text(colour = "white"),
        plot.background = element_rect(fill = NA, colour = "#2B3E50"), 
        axis.line = element_line(color = "white"),
        axis.text = element_text(color = "white", size = 10),
        axis.title.y = element_text(angle = 0, vjust = 0.5))

```

### Kaggle-LeaderBoard-Result: 3329.84561 in RMSE

## Update the parameters

As anticipated before, after computing the *k-fold* cross validation, we
are going to select some new parameter, similar to the ones obtained
from the validation, to be tuned via Nested Cross Validation.

```{r update parameter, warning=FALSE}
# Update Parameters -------------------------------------------------------------

# Using the best parameters
d_best <- best_params[1]
q_best <- best_params[2]
k_best <- best_params[3]
a_best <- best_params[4]
l_best <- best_params[5]
p_best <- best_params[6]

#Select the new parameters for the nested cross validation
d <- d_best
q <- q_best+seq(-2,2,1)
k <- k_best
a <- a_best + seq(-0.03, 0, 0.01)
l <- l_best + seq(-0.05, 0.05, 0.025)
p <- p_best

# Set the parameter for the CV
parameters <- list(d, q, k, a, l, p)

```

## Nested Cross Validation

To apply the Nested Cross Validation the first thing we have to do is
implement its main functions:

-   **nested_crossval**, that takes as input the hyperparameters and the
    data, inside this function we will iteratively select the *k-folds*
    randomly, then for each of the folds we call the *inner_crossval*
    function (to obtain the vector $e^{(in)}$), train and test the model
    (to compute the vector $e^{(out)}$). After each iteration the
    function will store $e^{(in)}$ in the vector $es$ and compute the
    wanted quantities $(mean(e^{(in)})-mean(e^{(out)}))^2$ and
    $\frac{var(e^{(out)})}{len(fold)}$, to be stored in *a_list* and
    *b_list* respectively. At the end the function will return an estimator
    of the RMSE and its relative MSE.

-   **inner_crossval**, takes as input $k-1$ out of the $k$ folds
    created in the previous function and compute the cross validation on
    the $k-1$ folds, returning as a output a vector of $k-1$ RMSEs.

```{r nested cross validation, eval=FALSE}

# Secondary function for nested CV
inner_crossval <- function(x, train_set){
  # define the parameters
  d <- x[1]
  q <- x[2]
  K <- x[3]
  a <- x[4]
  l <- x[5] 
  p <- x[6]
  
  # Pre-set the error
  e_in <- c()
  # Loop over k-1 fold
  for(k in (1:(K-1))){
    # Set the indices
    idx <- ((k-1)*l_folds+1): (k*l_folds)
    
    # Validation set
    cv_test <- train_set[idx,]
    # Train set
    cv_train <- train_set[-idx,]
    
    # Define the knots
    knots <- seq(0, p, length.out=q+1)[2:(q+1)]
    
    
    # Compute the power functions
    M_cv_train <- power_functions(d = d, q = q, knots = knots, x = cv_train$x)
    M_cv_test <-  power_functions(d = d, q = q, knots = knots, x = cv_test$x)
    
    # Compute the weights based on the conditional variance
    hat_weights <- compute_weights(knots = knots , cv_train)
      
    # Train the model
    cv_model <- glmnet(M_cv_train, 
                       cv_train$y,
                       family = "gaussian", 
                       alpha=a, 
                       lambda=l,
                       weights = hat_weights)
    
    # prediction
    cv_predictions <- predict(cv_model, M_cv_test)
    
    # Compute the error
    e_temp <- sqrt(cv_test$y-cv_predictions)^2)
    e_in <- c(e_in, e_temp)
  }
  
  return(e_in)
}

# Main function for the nested CV
nested_crossval <- function(x){
  # set the parameters
  d <- x[1]
  q <- x[2]
  K <- x[3]
  a <- x[4]
  l <- x[5] 
  p <- x[6]
  R <- 250
  
  # Define the length of the folds
  l_folds <<- nrow(train_set) / K
  
  # Define the vectors needed
  es <- c()
  a_list <- rep(NA, R*K)
  b_list <- rep(NA, R*K)
  
  # For R times
  for(r in (1:R)){
    # Select the indexes for the folds randomly
    idx <- sample((1:nrow(train_set)),nrow(train_set))
    
    # For each of the folds
    for(k in (1:K)){ 
      # Define train and validation selecting the folds
      cv_test <- train_set[idx[((k-1)*l_folds+1): (k*l_folds)],]
      cv_train <- train_set[-idx[((k-1)*l_folds+1): (k*l_folds)],]
      
      # inner cross
      e_in <- inner_crossval(x, cv_train)
      
      # Outer cross
      knots <- seq(0, p, length.out=q+1)[2:(q+1)]
      
      # Map in the feature matrices
      M_cv_train <- power_functions(d = d, q = q, knots = knots, x = cv_train$x)
      M_cv_test <-  power_functions(d = d, q = q, knots = knots, x = cv_test$x)
      
      # Compute the weights
      hat_weights <- compute_weights(knots = knots , cv_train)
      # Train the model
      cv_model <- glmnet(M_cv_train, 
                         cv_train$y,
                         family = "gaussian", 
                         alpha=a, 
                         lambda=l,
                         weights = hat_weights)
      # Compute the predictions
      cv_predictions <- predict(cv_model, M_cv_test)
      # Compute the error
      e_out <- sqrt((cv_test$y-cv_predictions)^2)
      
      # Store the needed informations
      es <- c(es, e_in)
      a_list[(r-1)*K+k] <- (mean(e_in)-mean(e_out))^2 
      b_list[(r-1)*K+k] <- (sd(e_out)^2)/l_folds
    }
  }
   # Compute the estimated error
  mse <- mean(a_list)-mean(b_list)
  err <- mean(es)
  
  return(paste(mse,err))
}

```

```{r run nested cv , eval = FALSE}
# CV nested --------------------------------------------------------------
cl = makeCluster(detectCores())
clusterExport(cl, c('train_set','compute_weights' ,'inner_crossval', 'power_functions', 'glmnet'))
res <- gridSearch(nested_crossval, levels=parameters, method = 'snow', cl=cl)
stopCluster(cl)

# save the results err and mse in two vectors
return_numeric <- function(x){
  return(as.numeric(strsplit(x, ' ')[[1]]))
} 

vector <- unlist(lapply(res$values, FUN = return_numeric))
n <- length(vector)
mse <- sqrt(vector[seq(n) %% 2 == 1])*1.96
err <- vector[seq(n) %% 2 == 0]
best_params <- res$levels[[which.min(mse)]]

```

After performed the nested cross validation we can compute a confidence
interval for each of the CV RMSE we got, by using the following
computation:

$$
\big( \widehat{RMSE} - z_{1-\alpha/2} \sqrt{\widehat{MSE}}  \ ,  \widehat{RMSE} + z_{1-\alpha/2} \sqrt{\widehat{MSE}} \ \big)
$$

For each combination tested, the following plot show the estimated RMSE
with the $1 - \alpha = 95\%$ confidence interval associated.

```{r plot err with mse CI, fig.showtext = TRUE, fig.height = 5, echo = FALSE,fig.width = 8, fig.align = "center", message=FALSE, warning=FALSE}
mse_err <- read.csv("RData/Nested_NODC.csv")
mse <- mse_err$mse
err <- mse_err$err

colors <- c("RMSE" = "#F1F1E6", "Confidence Interval" = "#FFB238")

ggplot()+
  geom_point(aes(x=1:length(mse),y=err, color='RMSE'), shape = 18, size= 1.8)+
  geom_errorbar(aes(x =1:length(mse), ymin=err-mse, ymax=mse+err, width=0.2, color='Confidence Interval'))+
  theme_bw() +
  scale_y_continuous(limits = c(1100, 2100)) +
  labs(x = "Combination of hyperparameters", y = "CI of RMSE", color = "    Legend", title = 'Confidence Interval of RMSE', shape = "", color="") +
  guides(color = guide_legend(override.aes=list(shape = c(NA, 18), size = 2, linetype = c(1, NA), linewidth=1.4))) +
  scale_color_manual(values = colors) +
  theme(panel.border = element_rect(color = "white"),
      legend.title = element_text(size=12, color='white'),
      plot.title = element_text(hjust = 0.5,size=16),
      legend.text = element_text(size = 10, color = 'white'),
      legend.key = element_rect(color = NA, fill = NA),
      legend.background = element_rect(fill = NA),
      panel.grid.major = element_line(linetype = 0, linewidth=.2),
      panel.grid.minor = element_line(linetype = 0),
      panel.background = element_rect(fill = NA),
      title = element_text(colour = "white"),
      plot.background = element_rect(fill = NA, colour = "#2B3E50"), 
      axis.line = element_line(color = "white"),
      axis.text = element_text(color = "white", size = 10))

```

## Comments:

Since the plot shown that all combination of hyperparameters have similar $\widehat{RMSE}$, the decision we made to select the best configuration is taking into account the width of the interval. At the end we decided to select the configutations that minimize the $\widehat{MSE}$. The hyperparameters we have selected are shown below.  

```{r training with nested parameters, warning = FALSE}
load("RData//best_params_NODC.RData")
# Prediction --------------------------------------------------------------

# Using the best parameters
d_best <- best_params[1]
q_best <- best_params[2]
k_best <- best_params[3]
a_best <- best_params[4]
l_best <- best_params[5]
p_best <- best_params[6]

# Compute the predictions
knots <- seq(0, p_best, length.out=q_best+1)[2:(q_best+1)]
M_train <- power_functions(d = d_best, q = q_best, knots = knots, x = train_set$x)
M_test <- power_functions( d = d_best, q = q_best, knots = knots, x = test_set_vero$x)
knots_test <- power_functions( d = d_best, q = q_best, knots = knots, x = knots)
hat_weights <- compute_weights(knots , train_set)
final_model <- glmnet(M_train, train_set$y, family ="gaussian", 
                      alpha=a_best, lambda=l_best , weights = hat_weights )
predictions <- predict(final_model, M_test)
round(best_params,3)
```

```{r plot after nested, fig.showtext = TRUE, echo = FALSE,fig.height = 5, fig.width = 10, fig.align = "center", message=FALSE, warning=FALSE}

colors <- c("Real Data" = "#F1F1E6", "Predicted" = "#E09F20", "Knots" = "#006054")
point_size <- 1.4
knots_size <- 1
knots_shape <- 18
line_size <- .6


ggplot() +
  geom_point(aes(x = train_set$x, y = train_set$y, color = 'Real Data'), size = point_size) +
  geom_line(aes(x = test_set_vero$x, y = predictions, color = 'Predicted'), size = 2) +
  geom_point(aes(x = knots, y = predict(final_model, knots_test), color = 'Knots'), shape = knots_shape, stroke = 1.7, size = knots_size) +
  theme_bw() +
  labs(x = "x", y = "y", color = "     Legend", title = 'Prediction on WMAP after nested cross validation', shape = "", color="") +
  scale_color_manual(values = colors) +
  guides(color = guide_legend(override.aes=list(shape = c(knots_shape, NA, 18), size = 2, linetype = c(NA, 1, NA), linewidth=1.8))) +
  theme(panel.border = element_rect(color = "white"),
        legend.title = element_text(size=12, color='white'),
        plot.title = element_text(hjust = 0.5,size=16),
        legend.text = element_text(size = 10, color = 'white'),
        legend.key = element_rect(color = NA, fill = NA),
        legend.background = element_rect(fill = NA),
        panel.grid.major = element_line(linetype = 2, linewidth=.2),
        panel.grid.minor = element_line(linetype = 0),
        panel.background = element_rect(fill = NA),
        title = element_text(colour = "white"),
        plot.background = element_rect(fill = NA, colour = "#2B3E50"), 
        axis.line = element_line(color = "white"),
        axis.text = element_text(color = "white", size = 10),
        axis.title.y = element_text(angle = 0, vjust = 0.5))
```

### Kaggle-LeaderBoard-Result: 3325,03931 in RMSE

# Conclusion & Final Remark

In conclusion, the use of splines for regression can be a powerful tool
for modeling nonlinear relationships between variables. Thanks to the use
of truncated power basis functions, the splines can effectively capture
the complex structure of the relations which can exist between the
variables $X$ and $Y$ without using a high degree polynomial which it
will lead a more accurate and flexible model.