test.Rmd

## Metrics

```{r metrics}
metrics <- network %>% evaluate(test_data_vec, test_labels_vec)
```

```{r}
metrics
metrics$acc
# Error rate: incorrect calling
1 - metrics$acc
```

## Predictions

```{r predictions}
network %>% predict_classes(test_data_vec[1:10,])
```

```{r allPredictions}
predictions <- network %>% predict_classes(test_data_vec)
actual <- unlist(test_labels)
totalmisses <- sum(predictions != actual)
```

# Confusion Matrix

```{r confusion, echo = F}
suppressPackageStartupMessages(library(tidyverse))
# library(dplyr)
data.frame(target = actual,
           prediction = predictions) %>% 
  filter(target != prediction) %>% 
  group_by(target, prediction) %>%
  count() %>%
  ungroup() %>%
  mutate(perc = n/nrow(.)*100) %>% 
  filter(n > 1) %>% 
  ggplot(aes(target, prediction, size = n)) +
  geom_point(shape = 15, col = "#9F92C6") +
  scale_x_continuous("Actual Target", breaks = 0:45) +
  scale_y_continuous("Prediction", breaks = 0:45) +
  scale_size_area(breaks = c(2,5,10,15), max_size = 5) +
  coord_fixed() +
  ggtitle(paste(totalmisses, "mismatches")) +
  theme_classic() +
  theme(rect = element_blank(),
        axis.line = element_blank(),
        axis.text = element_text(colour = "black"))

```



# As a Regression Problem

Convert z-scores:

$$z_i=\frac{x_i-\bar{x}}{s}$$

```{r zScores, cache = T}
# parameters for Scaling:
mean <- colMeans(train_data) # mean of each column
std <- apply(train_data, 2, sd) # stdev of each column

# Calculate feature-wise (within-variable) z-scores: (x - mean)/std
train_data <- scale(train_data, center = mean, scale = std)
test_data <- scale(test_data, center = mean, scale = std)
```


Let's see what the data looks like now:

```{r}

train_data %>% 
  as.data.frame() %>% 
  rename_all(funs(abalone_names[-9])) %>% 
  gather() %>% 
  ggplot(aes(key, value)) +
  geom_jitter(shape = 1, alpha = 0.3)

```

# Part 2: Define Network

## Define the network as a function

In contrast to our previous case studies, we're going to call the same model multiple times. So we'll create a function with no arguments that we can call to create our model when ever we want to use it for training. 

Here, I've hardcoded the number of features for this dataset (`13`). To generalize, we could just use `dim(train_data)[2]` to get the number of dimensions from the training set.  

```{r defModel}
build_model <- function() {
  network <- keras_model_sequential() %>% 
    layer_dense(units = 64, activation = "relu", input_shape = 8) %>% 
    layer_dense(units = 64, activation = "relu") %>% 
    layer_dense(units = 1) 
    
  network %>% compile(
    optimizer = "rmsprop", 
    loss = "mse", 
    metrics = c("mae")
  )
}
```

Note two new functions here, the mean squared error:

$$\operatorname{MSE} = \frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{n} = \frac{\sum_{i=1}^n{e_i^2}}{n}$$
and the mean absolute error (MAE):

$$\mathrm{MAE} = \frac{\sum_{i=1}^n\left| y_i-\hat{y_i}\right|}{n} = \frac{\sum_{i=1}^n\left| e_i\right|}{n}$$
where $\hat{y_i}$ is the predicted value, given in our last single-unit layer, and $y_i$ is the actual value, the label.

# Part 3: k-fold cross validation

```{r setkFold, echo = TRUE, results = 'hide'}
k <- 4 # four groups
indices <- sample(1:nrow(train_data)) # randomize the training set before splitting for k-fold cross validation:
folds <- cut(indices, breaks = k, labels = FALSE) # divide the ordered indices into k intervals, labelled 1:k.
```

```{r kfold100, cache = T}
num_epochs <- 100
all_scores <- c() # An empty vector to store the results from evaluation

for (i in 1:k) {
  cat("processing fold #", i, "\n")
  # Prepare the validation data: data from partition # k
  val_indices <- which(folds == i, arr.ind = TRUE) 
  
  # validation set: the ith partition
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Training set: all other partitions
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Call our model function (see above)
  network <- build_model()
  
  # summary(model)
  # Train the model (in silent mode, verbose=0)
  network %>% fit(partial_train_data,
                  partial_train_targets,
                  epochs = num_epochs,
                  batch_size = 1,
                  verbose = 0)
                
  # Evaluate the model on the validation data
  results <- network %>% evaluate(val_data, val_targets, verbose = 0)
  all_scores <- c(all_scores, results$mean_absolute_error)
}  
```

We get 4 mae values

```{r allscores}
all_scores
```

### Training for 500 epochs

Let's try training the network for a bit longer: 500 epochs. To keep a record of how well the model did at each epoch, we will modify our training loop to save the per-epoch validation score log:

```{r clearMem}
# Some memory clean-up
K <- backend()
K$clear_session()
```

Train our the models:

```{r kfold500, echo = T, results = 'hide', cache = T}
num_epochs <- 500
all_mae_histories <- NULL # an empty object to cumulatively store the model metrics

for (i in 1:k) {
  cat("processing fold #", i, "\n")
  
  # Prepare the validation data: data from partition # k
  val_indices <- which(folds == i, arr.ind = TRUE)
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Prepare the training data: data from all other partitions
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Build the Keras model (already compiled)
  model <- build_model()
  
  # Train the model (in silent mode, verbose=0)
  history <- model %>% fit(partial_train_data, 
                           partial_train_targets,
                           validation_data = list(val_data, val_targets),
                           epochs = num_epochs, 
                           batch_size = 1, 
                           verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}
```

Calculate the average per-epoch MAE score for all folds:

```{r plot1}
average_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, mean)
)

p <- ggplot(average_mae_history, aes(x = epoch, y = validation_mae))

p + 
  geom_point()

p + 
  geom_smooth(method = 'loess', se = FALSE)
```

According to this plot, it seems that validation MAE stops improving significantly after circa 80 epochs. Past that point, we start overfitting.

Once we are done tuning other parameters of our model (besides the number of epochs, we could also adjust the size of the hidden layers), we can train a final "production" model on all of the training data, with the best parameters, then look at its performance on the test data:

```{r runZ, echo = F, results = 'hide', cache = T}
# Get a fresh, compiled model.
model <- build_model()

# Train it on the entirety of the data.
model %>% fit(train_data, 
              train_targets,
              epochs = 80, 
              batch_size = 16, 
              verbose = 0)

result <- model %>% evaluate(test_data, test_targets)
```

```{r resultsZ}
result
```

We are still off by about `r round(result$mean_absolute_error * 1000)`.

## Alternatives: No Normalization

Let's imagine that we didn't normalize the input variables

```{r setupNone, echo = F, cache = T}
# Obtain the raw data
c(c(train_data, train_targets), c(test_data, test_targets)) %<-% dataset_boston_housing()

# Some memory clean-up
K <- backend()
K$clear_session()

num_epochs <- 500
all_mae_histories <- NULL # an empty object to cumulatively store the model metrics

for (i in 1:k) {
  cat("processing fold #", i, "\n")
  
  # Prepare the validation data: data from partition # k
  val_indices <- which(folds == i, arr.ind = TRUE)
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Prepare the training data: data from all other partitions
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Build the Keras model (already compiled)
  model <- build_model()
  
  # Train the model (in silent mode, verbose=0)
  history <- model %>% fit(partial_train_data, 
                           partial_train_targets,
                           validation_data = list(val_data, val_targets),
                           epochs = num_epochs, 
                           batch_size = 1, 
                           verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}
```

Calculate the average per-epoch MAE score for all folds:

```{r plot2, echo = F, cache = T}
average_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, mean)
)

p <- ggplot(average_mae_history, aes(x = epoch, y = validation_mae))

p + 
  geom_point()

p + 
  geom_smooth(method = 'loess', se = FALSE)
```

The validation MAE stops improving significantly after circa 140 epochs.

```{r runNone, echo = F, results = 'hide', cache = T}
# Get a fresh, compiled model.
model <- build_model()

# Train it on the entirety of the data.
model %>% fit(train_data, 
              train_targets,
              epochs = 140, 
              batch_size = 16, 
              verbose = 0)

result_none <- model %>% evaluate(test_data, test_targets)
```

```{r resultsNone}
result_none
```

Now, without any normalization, we're off by about `r round(result_none$mean_absolute_error * 1000)`, compared to `r round(result$mean_absolute_error * 1000)` previously.

## Alternatives: 0-1 normalization 

How about if we did 0-1 normalization?

$$z_i=\frac{x_i-\min(x)}{\max(x)-\min(x)}$$

```{r setup01, echo = F, cache = T}
# Obtain the raw data
c(c(train_data, train_targets), c(test_data, test_targets)) %<-% dataset_boston_housing()

train_data <- apply(train_data, 2, function(x) (x-min(x))/(max(x)-min(x)))
test_data <- apply(test_data, 2, function(x) (x-min(x))/(max(x)-min(x)))

# Some memory clean-up
K <- backend()
K$clear_session()

num_epochs <- 500
all_mae_histories <- NULL # an empty object to cumulatively store the model metrics

for (i in 1:k) {
  cat("processing fold #", i, "\n")
  
  # Prepare the validation data: data from partition # k
  val_indices <- which(folds == i, arr.ind = TRUE)
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Prepare the training data: data from all other partitions
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Build the Keras model (already compiled)
  model <- build_model()
  
  # Train the model (in silent mode, verbose=0)
  history <- model %>% fit(partial_train_data, 
                           partial_train_targets,
                           validation_data = list(val_data, val_targets),
                           epochs = num_epochs, 
                           batch_size = 1, 
                           verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}
```

Calculate the average per-epoch MAE score for all folds:

```{r plot3, echo = F, cache = T}
average_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, mean)
)

p <- ggplot(average_mae_history, aes(x = epoch, y = validation_mae))

p + 
  geom_point()

p + 
  geom_smooth(method = 'loess', se = FALSE)
```

The validation MAE stops improving significantly after circa 140 epochs.

```{r run01, echo = F, results = 'hide', cache = T}
# Get a fresh, compiled model.
model <- build_model()

# Train it on the entirety of the data.
model %>% fit(train_data, 
              train_targets,
              epochs = 140, 
              batch_size = 16, 
              verbose = 0)

result_01 <- model %>% evaluate(test_data, test_targets)
```

```{r results01}
result_01
```

Here, we're off by about `r round(result_01$mean_absolute_error * 1000)`, compared to `r round(result$mean_absolute_error * 1000)` with z scores.













## old


An example of the data

```{r}
train_data[[1]]
```

## Prepare the data:

Let's look at the first example from the training set. Recall that these are the index positions of the words 

```{r}
train_example <- sort(unique(train_data[[1]]))
train_example
```

Now we have a large matrix, where each row is 10000 elements long. Wherever we have a value in the above data set, the matrix has a 1

```{r}
# Just the first 100 values in the first entry (row)
train_data_vec[1,1:100]
```

We can confirm this by counting the values:

```{r}
sum(train_data_vec[1,]) == length(train_example)
```

The position of the 1s corresponds to the indices above:

```{r}
which(as.logical(train_data_vec[1,]))
```

## Prepare labels:

The `_labels` objects contain the news wire labels. Each newswire can only have one *label* (i.e. "sigle-label"), from a total of 46 possible *classes* (i.e. "multi-class"). The classes are just given numerical values (0 - 45), it doesn't matter what they are actually called, although that information would be helpful in understanding mis-labeling.

```{r strLabelsPre}
str(train_labels)
sort(unique(train_labels))
```

Some classes are very common, which we'll see play out in our confusion matrix below 

```{r plotLabelsPre}
# Note plyr not dplyr here. I'm just using a shortcut
library(ggplot2)
train_labels %>% 
  plyr::count() %>%
  ggplot(aes(x, freq)) +
  geom_col()
```

The distribution of the test and training set should be roughly equivalent, so let's have a look. 

```{r}
data.frame(x = train_labels) %>% 
  group_by(x) %>% 
  summarise(train_freq = 100 * n()/length(train_data)) -> train_labels_df

data.frame(x  = test_labels) %>% 
  group_by(x) %>% 
  summarise(train_freq = 100 * n()/length(test_data)) %>% 
  inner_join(train_labels_df, by="x") %>% 
  gather(key, value, -x) %>% 
  ggplot(aes(x, value, fill = key)) +
  geom_col(position = "dodge") +
  scale_y_continuous("Percentage", limits = c(0,40), expand = c(0,0)) +
  scale_x_continuous("Label", breaks = 0:45, expand = c(0,0)) +
  scale_fill_manual("", labels = c("test","train"), values = c("#AEA5D0", "#54C8B7")) +
  theme_classic() +
  theme(legend.position = c(0.8, 0.8),
        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black"))
```

We treat these just like how we treated the MNIST labels in the previous unit. We make the format match the output we expect to get from softmax so that we can make a direct comparison.

```{r prepLabels}
train_labels_vec <- to_categorical(train_labels)
test_labels_vec <- to_categorical(test_labels)
```

```{r strLabelsPost}
str(train_labels_vec)
str(test_labels_vec)
```

Notice the similiarity to how we prepared our training data with one-hot encoding, both are sparse matrices. Each row in the training data contains a 1 at the position where that word is present, and in the labels, each row contains a maximum of a single 1, indicating the class. 

# Part 2: Define Network

## Define the network

```{r architecture}
network <- keras_model_sequential() %>% 
  layer_dense(units = 64, activation = "relu", input_shape = c(10000)) %>% 
  layer_dense(units = 64, activation = "relu") %>% 
  layer_dense(units = 46, activation = "softmax")

```

## View a summary of the network

```{r summary}
summary(network)
```

## Compile

```{r compile}
network %>% compile(
  optimizer = "rmsprop",
  loss = "categorical_crossentropy",
  metrics = c("accuracy")
)
```

# Part 3: Validate our approach

Let's set apart 1,000 samples in our training data to use as a validation set:

```{r}
index <- 1:1000

val_data_vec <- train_data_vec[index,]
train_data_vec <- train_data_vec[-index,]

val_labels_vec <- train_labels_vec[index,]
train_labels_vec = train_labels_vec[-index,]
```

Now let's train our network for 20 epochs:

```{r echo=TRUE, results = "hide", warning = FALSE}
history <- network %>% fit(
  train_data_vec,
  train_labels_vec,
  epochs = 20,
  batch_size = 512,
  validation_data = list(val_data_vec, val_labels_vec)
)
```

Let's display its loss and accuracy curves:

```{r}
plot(history)
```

The network begins to overfit after nine epochs. Let's train a new network from scratch for nine epochs and then evaluate it on the test set.

```{r, echo=TRUE, results='hide'}
network <- keras_model_sequential() %>% 
  layer_dense(units = 64, activation = "relu", input_shape = c(10000)) %>% 
  layer_dense(units = 64, activation = "relu") %>% 
  layer_dense(units = 46, activation = "softmax")
  
network %>% compile(
  optimizer = "rmsprop",
  loss = "categorical_crossentropy",
  metrics = c("accuracy")
)

history <- network %>% fit(
  train_data_vec,
  train_labels_vec,
  epochs = 9,
  batch_size = 512,
  validation_data = list(val_data_vec, val_labels_vec)
)
```

# Part 4: Using sparse categorical crossentropy

Above, we vectorized the labels, like what we did with the MNIST data set. Alternatively, we could have just used the original integer values. To showcase this, let's create a new network, `network_int`, so that we don't mix up our results. The network architecture is the same:

```{r}
network_int <- keras_model_sequential() %>% 
  layer_dense(units = 64, activation = "relu", input_shape = c(10000)) %>% 
  layer_dense(units = 64, activation = "relu") %>% 
  layer_dense(units = 46, activation = "softmax")
```

Here, the only thing we need to chance is the loss function. `categorical_crossentropy`, expects the labels to follow a categorical encoding, but `sparse_categorical_crossentropy` expects integer labels. 

```{r}
network_int %>% compile(
  optimizer = "rmsprop",
  loss = "sparse_categorical_crossentropy",
  metrics = c("accuracy")
)
```

Before we train the model, let's make a validation set, like we did above. We'll use the original training set for this.

```{r}
val_train_labels <- train_labels[index]
train_labels <- train_labels[-index]
```

Now let's train our model `network_int` using the integer data, instead of the vectorized data:

```{r}
history_int <- network_int %>% fit(
  train_data_vec,
  train_labels,
  epochs = 9,
  batch_size = 512,
  validation_data = list(val_data_vec, val_train_labels)
)
```

This new loss function is mathematically the same as `categorical_crossentropy`. It just has a different interface. When we look at our metrics below we'll use the original model, that accessed the vectorized data. If you want to use `network_int` make sure you use the original integer labels of the test set, `test_labels`, not `test_labels_vec`. 

# Part 5: Check output

Let's return to our original model using the vectorized data:

## Metrics

```{r metrics}
metrics <- network %>% evaluate(test_data_vec, test_labels_vec)
```

```{r}
metrics
metrics$acc
# Error rate: incorrect calling
1 - metrics$acc
```

## Predictions

```{r predictions}
network %>% predict_classes(test_data_vec[1:10,])
```

```{r allPredictions}
predictions <- network %>% predict_classes(test_data_vec)
actual <- unlist(test_labels)
totalmisses <- sum(predictions != actual)
```

# Confusion Matrix

```{r confusion, echo = F}
suppressPackageStartupMessages(library(tidyverse))
# library(dplyr)
data.frame(target = actual,
           prediction = predictions) %>% 
  filter(target != prediction) %>% 
  group_by(target, prediction) %>%
  count() %>%
  ungroup() %>%
  mutate(perc = n/nrow(.)*100) %>% 
  filter(n > 1) %>% 
  ggplot(aes(target, prediction, size = n)) +
  geom_point(shape = 15, col = "#9F92C6") +
  scale_x_continuous("Actual Target", breaks = 0:45) +
  scale_y_continuous("Prediction", breaks = 0:45) +
  scale_size_area(breaks = c(2,5,10,15), max_size = 5) +
  coord_fixed() +
  ggtitle(paste(totalmisses, "mismatches")) +
  theme_classic() +
  theme(rect = element_blank(),
        axis.line = element_blank(),
        axis.text = element_text(colour = "black"))

```

## from applied predictive modelling

- Plot the data to assess the functional relationships between the predictors and the outcome.
- Use scatter plots and correlation plots to understand how the predictors relate to one another.
- Estimate variable importance scores for each predictor. Develop an approach to determining a reduced set of nonredundant predictors.
- Apply principal component analysis to the continuous predictors to determine how many distinct underlying pieces of information are in the data. Would feature extraction help these data?