```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


####Analytics Class####

## R programming skillset _ essential pipeline of kernel regression

```{r}
# Step 1 -> Read data into R workstation

#### Read data from a CSV file 
# filename
setwd("C:/Users/shuai/Downloads")
data <- read.csv(file = "KR.csv", header = TRUE)
# str(data)
```

```{r}
# Step 2 -> Data preprocessing
# Create your X matrix (predictors) and Y vector (outcome variable)
X <- data$x
Y <- data$y

# Make sure the outcome variable is legitimate. If it is a continuous variable (regression problem), it should be defined as a "num" variable in R. If it is a binary or a more genernal categorical variable (classification problem), it should be defined as a "factor" variable in R. 

# Create a training data (half the original data size)
train.ix <- sample(nrow(data),floor( nrow(data) * 4/5) )
data.train <- data[train.ix,]
# Create a testing data (half the original data size)
data.test <- data[-train.ix,]

```

```{r}
# Step 3 -> gather a list of candidate models
# Kernel regression: often to compare models with different kernels, different values of the bandwidth of the kernel, different set of variables

# Use different values of bandwidth

# model1: ksmooth(x,y, kernel = "normal", bandwidth=10)
# model2: ksmooth(x,y, kernel = "box", bandwidth=10)

```


```{r}
# Step 4 -> Use 5-fold cross-validation to evaluate all the models

# First, let me use 5-fold cross-validation to evaluate the performance of model1
n_folds = 10 # number of fold (the parameter K as we say, K-fold cross validation)
N <- dim(data.train)[1] # the sample size, N, of the dataset
folds_i <- sample(rep(1:n_folds, length.out = N)) # This randomly creates a labeling vector (1 X N) for the N samples. For example, if N = 16, and I run this function and it returns the value as 5  4  4 10  6  7  6  8  3  2  1  5  3  9  2  1. That means, the first sample is allocated to the 5th fold, the 2nd and 3rd samples are allocated to the 4th fold, etc. 

cv_mse <- NULL # cv_mse aims to make records of the prediction error for each fold
for (k in 1:n_folds) {
  test_i <- which(folds_i == k) # In each iteration of the n_folds iterations, remember, we use one fold of                                      data as the testing data
  data.train.cv <- data.train[-test_i, ] # Then, the remaining n_folds-1 folds' data form our training data
  data.test.cv <- data.train[test_i, ]   # This is the testing data, from the ith fold
  require( 'kernlab' )
  model1 <- ksmooth(data.train.cv$x, data.train.cv$y, kernel = "normal", bandwidth = 20, x.points=data.test.cv[,1]) # (1) Fit the kernel regression model with Gaussian kernel (argument: kernel = "normal") and bandwidth = 0.5; (2) Here, one unique thing about ksmooth is, there is no predict() for it. Rather, it has the argument "x.points=data.test.cv" to specify where you want to predict on 
  y_hat <- model1$y  # Predict on the testing data using the trained model
  true_y <- data.test.cv$y                  # get the true y values for the testing data
  cv_mse[k] <- mean((true_y - y_hat)^2)    # mean((true_y - y_hat)^2): mean squared error (MSE). The small                                              this error, the better your model is
}
mean(cv_mse)


cv_mse <- NULL # cv_mse aims to make records of the prediction error for each fold
for (k in 1:n_folds) {
  test_i <- which(folds_i == k) # In each iteration of the n_folds iterations, remember, we use one fold of                                      data as the testing data
  data.train.cv <- data.train[-test_i, ] # Then, the remaining n_folds-1 folds' data form our training data
  data.test.cv <- data.train[test_i, ]   # This is the testing data, from the ith fold
  require( 'kernlab' )
  model2 <- ksmooth(data.train.cv$x, data.train.cv$y, kernel = "box", bandwidth = 20, x.points=data.test.cv[,1]) # (1) Fit the kernel regression model with Gaussian kernel (argument: kernel = "normal") and bandwidth = 0.5; (2) Here, one unique thing about ksmooth is, there is no predict() for it. Rather, it has the argument "x.points=data.test.cv" to specify where you want to predict on 
  y_hat <- model2$y  # Predict on the testing data using the trained model
  true_y <- data.test.cv$y                  # get the true y values for the testing data
  cv_mse[k] <- mean((true_y - y_hat)^2)    # mean((true_y - y_hat)^2): mean squared error (MSE). The small                                              this error, the better your model is
}
mean(cv_mse)

```

```{r}
#Remark: You could use visual inspection to decide what is the best kernel and bandwidth for your kernel regression as well. Sometimes, visual inspection (sounds subjective) is effective, and can look at the global pattern more clearly than automatic cross-validation that relies on one single index (mse or prediction error)

plot(Y ~ X, col = "gray", lwd = 2)
lines(ksmooth(X,Y, "normal", bandwidth=2),lwd = 3, col = "darkorange")
lines(ksmooth(X,Y, "normal", bandwidth=10),lwd = 3, col = "dodgerblue4")
lines(ksmooth(X,Y, "normal", bandwidth=30),lwd = 3, col = "forestgreen")
lines(ksmooth(X,Y, "normal", bandwidth=100),lwd = 3, col = "black")
legend(x = "topright", legend = c("Kernel Reg (bw = 2)", "Kernel Reg (bw = 10)", "Kernel Reg (bw = 30)","Kernel Reg (bw = 100)"), lwd = rep(3, 4), col = c("darkorange", "dodgerblue4", "forestgreen","black"), 
       text.width = 32, cex = 0.85)
```


```{r}
# Step 5 -> After model selection, use ksmooth() function to build your final model
kr.final <- ksmooth(data.train$x, data.train$y, kernel = "normal", bandwidth = 30, x.points=data.test[,1]) # 

```


```{r}
# Step 6 -> Evaluate the prediction performance of your SVM model
y_hat <- kr.final$y  # Predict on the testing data using the trained model
true_y <- data.test$y                  # get the true y values for the testing data
mse <- mean((true_y - y_hat)^2)    # mean((true_y - y_hat)^2): mean squared error (MSE). The small                                              this error, the better your model is
print(mse)
```
