---
title: "Regression to predict Ames, IA Home Sales Prices"
output:
  html_notebook:
    toc: yes
    toc_float: true
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
ggplot2::theme_set(ggplot2::theme_minimal())
```

In this case study, our objective is to predict the sales price of a home. This 
is a _regression_ problem since the goal is to predict any real number across
some spectrum (\$119,201, \$168,594, \$301,446, etc). To predict the sales 
price, we will use numeric and categorical features of the home.

As you proceed, you'll work through the steps we discussed in the last module:

1. Prepare data
2. Balance batch size with a default learning rate
3. Tune the adaptive learning rate optimizer
4. Add callbacks to control training
5. Explore model capacity
6. Regularize overfitting
7. Repeat steps 1-6
8. Evaluate final model results

# Package Requirements

```{r load-pkgs}
library(keras3)     # for deep learning
library(testthat)  # unit testing
library(tidyverse) # for dplyr, ggplot2, etc.
library(rsample)   # for data splitting
library(recipes)   # for feature engineering
library(AmesHousing) # The Ames Iowa Housing Data
```

# Step 0: Our Data

## The Ames housing dataset

For this case study we will use the [Ames housing dataset](http://jse.amstat.org/v19n3/decock.pdf) 
provided by the __AmesHousing__ package.

```{r get-data}
ames <- AmesHousing::make_ames()
dim(ames)
```

## Understanding our data

This data has been partially cleaned up and has no missing data:

```{r}
sum(is.na(ames))
```

But this tabular data is a combination of numeric and categorical data that we
need to address.

```{r ames-structure}
str(ames)
```

The numeric variables are on different scales. For example:

```{r numeric-ranges}
ames %>%
  select(Lot_Area, Lot_Frontage, Year_Built, Gr_Liv_Area, Garage_Cars, Mo_Sold) %>%
  gather(feature, value) %>%
  ggplot(aes(feature, value)) +
  geom_boxplot() +
  scale_y_log10(labels = scales::comma)
```

There are categorical features that could be ordered:

```{r numeric-categories}
ames %>%
  dplyr::select(matches("(Qual|Cond|QC|Qu)$")) %>%
  str()
```

And some of the categorical features have many levels:

```{r}
ames %>%
  select_if(~ is.factor(.) & length(levels(.)) > 8) %>%
  glimpse()
```

Consequently, our first challenge is transforming this dataset into numeric
tensors that our model can use.

# Step 1: Prep the Data

## Create train & test splits

One of the first things we want to do is create a train and test set as you
probably noticed that we do not have a train and test set similar to how MNIST 
was already set up for us. We can use the __rsample__ package to create our
train and test datasets.

__Note__: This will randomly select the 70/30 split so we are randomizing our
data with this process.

```{r}
set.seed(123)
ames_split <- initial_split(ames, prop = 0.7)
ames_train <- analysis(ames_split)
ames_test <- assessment(ames_split)

dim(ames_train)
dim(ames_test)
```

## Vectorize and scaling

All inputs and response values in a neural network must be tensors of either 
floating-point or integer data. Moreover, our feature values should not be
relatively large compared to the randomized initial weights _and_ all our 
features should take values in roughly the same range.

Consequently, we need to ___vectorize___ our data into a format conducive to neural 
networks [ℹ️](http://bit.ly/dl-02#3). For this data set, we'll transform our
data by:

1. removing any zero-variance (or near zero-variance) features
2. condensing unique levels of categorical features to "other"
3. ordinal encoding the quality features
4. normalize numeric feature distributions
5. standardizing numeric features to mean = 0, std dev = 1
6. one-hot encoding remaining categorical features

__Note__: we're using the recipes package (https://tidymodels.github.io/recipes)

```{r}
# Check help for each step, e.g., ?step_nzv
blueprint <- recipe(Sale_Price ~ ., data = ames_train) %>%
  step_nzv(all_nominal()) %>%                                       # step #1
  step_other(all_nominal(), threshold = .01, other = "other") %>%   # step #2
  step_integer(matches("(Qual|Cond|QC|Qu)$")) %>%                   # step #3
  step_YeoJohnson(all_numeric(), -all_outcomes()) %>%               # step #4
  step_center(all_numeric(), -all_outcomes()) %>%                   # step #5
  step_scale(all_numeric(), -all_outcomes()) %>%                    # step #5
  step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE)        # step #6

blueprint
```

This next step computes any relavent information (mean and std deviation of
numeric features, names of one-hot encoded features) on the training data so
there is no information leakage from the test data.

```{r}
prepare <- prep(blueprint, training = ames_train)
prepare
```

We can now vectorize our training and test data. If you scroll through the data
you will notice that all features are now numeric and are either 0/1 (one hot
encoded features) or have mean 0 and generally range between -3 and 3.

```{r}
baked_train <- bake(prepare, new_data = ames_train)
baked_test <- bake(prepare, new_data = ames_test)

# unit testing to ensure all columns are numeric
expect_equal(map_lgl(baked_train, ~ !is.numeric(.)) %>% sum(), 0)
expect_equal(map_lgl(baked_test, ~ !is.numeric(.)) %>% sum(), 0)

baked_train
```

Lastly, we need to create the final feature and response objects for train and 
test data. Since __keras__ and __tensorflow__ require our features & labels to be 
seperate objects we need to separate them. In doing so, our features need to be 
a 2D tensor which is why we apply `as.matrix` and our response needs to be a 
vector which is why we apply `pull`.

```{r}
x_train <- select(baked_train, -Sale_Price) %>% as.matrix()
y_train <- baked_train %>% pull(Sale_Price)

x_test <- select(baked_test, -Sale_Price) %>% as.matrix()
y_test <- baked_test %>% pull(Sale_Price)

# unit testing to x & y tensors have same number of observations
expect_equal(nrow(x_train), length(y_train))
expect_equal(nrow(x_test), length(y_test))
```

Our final feature set now has 188 input variables:

```{r}
dim(x_train)
dim(x_test)
```

# Step 2: Balance batch size with a default learning rate

To get started, let's build a simple model with...

- a single layer model with 128 units in the hidden layer. We have 188 features
  and 1 response node so a good starting point is mean(c(188, 1)) and then round
  up to the nearest value in the $2^s$ range (i.e. 32, 64, 128, 256, 512).
- a basic SGD optimizer
- use a mean square logarithmic error ("msle")
- also track the mean absolute error metric ("mae")
- 20% validation split

Now, start with the default batch size of 32 and then compare with smaller
values (i.e. 16) and larger values (i.e. 128). You're looking to balance the
progression of the loss learning curve and the training speed.

**Question:** Describe your observations how the batch size affects training.

```{r eval=FALSE}
n_feat <- ncol(x_train)

model <- keras_model_sequential() %>% 
  layer_dense(units = ____, activation = "relu", input_shape = ____) %>%
  layer_dense(units = ____)

model %>% compile(
    optimizer = ____,
    loss = ____,
    metrics = ____
  )

history <- model %>% fit(
  x_train,
  y_train,
  ____ ,
  validation_split = 0.2
)
```

```{r}
history
```

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

# Step 3: Tune the adaptive learning rate optimizer

Now go head and start assessing different adaptive learning rates such as:

- SGD+momentum
- RMSprop
- Adamp

Try a variety of learning rates. Recall that we typically start assessing rates
on a logarithmic scale (i.e. 0.1, 0.01, ..., 0.0001).

**Question:** Describe your observations of the effect of varying the learning rate for various optimizers.


```{r eval=FALSE}
model <- keras_model_sequential() %>% 
  layer_dense(units = ____, activation = ____, input_shape = ____) %>%
  layer_dense(units = 1)

model %>% compile(
    optimizer = ____,
    loss = "msle",
    metrics = "mae"
  )

history <- model %>% fit(
  x_train,
  y_train,
  batch_size = ____,
  validation_split = 0.2
)
```

```{r}
history
```

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

# Step 4: Add callbacks to control training

Add the following callbacks and see if your performance improves:

- early stopping with `patience = 3` and `min_delta = 0.00001`
- learning rate reduction upon a plateau with `patience = 1`

**Question:** List the chosen callback functions.

```{r eval=FALSE}
model <- keras_model_sequential() %>% 
  layer_dense(units = 128, activation = "relu", input_shape = ____) %>%
  layer_dense(units = 1)

model %>% compile(
    optimizer = ____,
    loss = "msle",
    metrics = "mae"
  )

history <- model %>% fit(
  x_train,
  y_train,
  batch_size = ____,
  epochs = ____,
  validation_split = 0.2,
  callbacks = list(
    # <add early stopping callback>,
    # <add learning rate reduction callback>
  )
)
```

```{r}
history
```

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

Plotting the learning rate shows that it reduced multiple times during training:

```{r}
plot(history$metrics$lr)
```


# Step 5: Explore model capacity

Now start to explore different widths and depths to your model.

- Assess a single layer with 128, 256, 512, and 1024 nodes
- Assess 1, 2, and 3 hidden layers

**Question:** Describe your observations how the number of nodes and layers affect training.

```{r eval=FALSE}
model <- keras_model_sequential() %>% 
  layer_dense(units = ____, activation = "relu", input_shape = n_feat) %>%
  # <add more layers>
  layer_dense(units = 1)

model %>% compile(
    optimizer = ____,
    loss = "msle",
    metrics = "mae"
  )

history <- model %>% fit(
  x_train,
  y_train,
  batch_size = ____,
  epochs = ____,
  validation_split = 0.2,
  callbacks = list(
    callback_early_stopping(patience = ____, min_delta = _____),
    callback_reduce_lr_on_plateau(patience = ____)
  )
)
```

```{r}
history
```

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

# Step 6: Regularize overfitting

If your model is overfitting, try to add...

- weight decay (i.e. `kernel_regularizer = regularizer_l2(l = xxx)`). Remember,
  we typically start by assessing values on logarithmic scale [0.1, 0.00001].
- dropout (`layer_dropout()`) between each layer. Remember, dropout rates
  typically range from 20-50%.

**Question:** Describe your observations how different regularization methods help to combat overfitting.

```{r eval=FALSE}
model <- ____() %>% 
  layer_dense() %>%
  # <add more layers to match your last model>
  layer_dense(units = 1)

model %>% ____(
    optimizer = ____,
    loss = ____,
    metrics = ____
  )

history <- model %>% ____(
  x_train,
  y_train,
  batch_size = ____,
  epochs = ____,
  validation_split = 0.2,
  callbacks = list(
    callback_____,
    callback_____
  )
)
```

```{r}
history
```
