---
title: "Neural Networks"
execute:
  echo: true
  eval: false
  output: false
---

This document collects the code from the slides so you can run it yourself. We predict median house values in Boston with a small neural network and compare it against a linear regression. The `Boston` data ships with the `MASS` package, so there is nothing to download.

```{r}
library(neuralnet)
library(MASS)
library(tidyverse)

set.seed(3456)
```

## The data

`Boston` has 506 rows, one per neighbourhood, and 14 numeric columns. The target is `medv`, the median house value in \$1000s. The other 13 columns are the predictors, e.g. the crime rate (`crim`), the average number of rooms (`rm`) or the share of lower-status population (`lstat`):

```{r}
head(Boston)
summary(Boston)
```

## Train and test split

We keep 75% of the rows for training and hold out the rest to evaluate the models on data they have not seen:

```{r}
train_i <- sample(1:nrow(Boston), round(0.75 * nrow(Boston)))

train_raw <- Boston[train_i, ]
test_raw <- Boston[-train_i, ]
```

## Scaling

The weights of a network start small and random, so a feature measured in the hundreds (`tax`) and one measured between 0 and 1 (`nox`) would not get a fair hearing. We therefore rescale every column to $[0, 1]$ with min-max scaling, $x' = (x - x_{\min}) / (x_{\max} - x_{\min})$. Linear regression does not care about scaling, but we use the same scaled data for both models so their errors are comparable.

The minimum and maximum are taken from the **training** data only and then applied to the test data. If we used the whole data set, information about the test data would leak into training. Because of this, a few test values can end up slightly below 0 or above 1, which is fine.

```{r}
mins <- apply(train_raw, 2, min)
maxs <- apply(train_raw, 2, max)

train <- as.data.frame(scale(train_raw, center = mins, scale = maxs - mins))
test <- as.data.frame(scale(test_raw, center = mins, scale = maxs - mins))

summary(train) # every column now runs from 0 to 1
```

## Fitting the network

`neuralnet` needs the formula written out, `medv ~ crim + zn + ...`; it does not understand `medv ~ .`. We build it from the column names:

```{r}
f <- as.formula(paste("medv ~", paste(setdiff(names(Boston), "medv"),
                                      collapse = " + ")))
f
```

`hidden = c(5, 3)` gives two hidden layers with 5 and 3 neurons. `linear.output = TRUE` leaves the output neuron without an activation function, which is what we want for regression: the prediction can be any number rather than being squished into $(0, 1)$. The hidden neurons use the logistic (sigmoid) activation by default.

```{r}
net <- neuralnet(f, data = train, hidden = c(5, 3), linear.output = TRUE)

plot(net) # weights on the arrows, biases in blue
```

The network has $13 \cdot 5 + 5 \cdot 3 + 3 \cdot 1 = 83$ weights plus $5 + 3 + 1 = 9$ biases, so 92 parameters, all found by gradient descent with backpropagation.

## Comparing against OLS

We predict the test data with the network and with a linear regression fitted on the same training data, and compare the root mean squared error:

```{r}
pred_net <- predict(net, test[, setdiff(names(test), "medv")])
pred_lm <- predict(glm(medv ~ ., data = train), test)

c(net = caret::RMSE(pred_net, test$medv),
  lm = caret::RMSE(pred_lm, test$medv))
```

These errors are on the $[0, 1]$ scale. To read them in \$1000s, we undo the scaling of `medv`:

```{r}
unscale <- function(x) x * (maxs["medv"] - mins["medv"]) + mins["medv"]

c(net = caret::RMSE(unscale(pred_net), test_raw$medv),
  lm = caret::RMSE(unscale(pred_lm), test_raw$medv))
```

Plotting predicted against actual values shows where each model goes wrong. Points on the diagonal are perfect predictions:

```{r}
tibble(
  actual = test_raw$medv,
  `Neural net` = unscale(pred_net[, 1]),
  `Linear model` = unscale(pred_lm)
) %>%
  pivot_longer(-actual, names_to = "model", values_to = "predicted") %>%
  ggplot(aes(x = actual, y = predicted, color = model)) +
  geom_point(alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0) +
  labs(x = "Actual medv", y = "Predicted medv", color = NULL) +
  theme_minimal()
```

## How stable is the result?

The net usually wins here, but only by a little. Both the split and the starting weights are random, so rerun the whole comparison with a few different seeds and see how much the margin moves:

```{r}
compare <- function(seed) {
  set.seed(seed)
  train_i <- sample(1:nrow(Boston), round(0.75 * nrow(Boston)))
  mins <- apply(Boston[train_i, ], 2, min)
  maxs <- apply(Boston[train_i, ], 2, max)
  scaled <- as.data.frame(scale(Boston, center = mins, scale = maxs - mins))
  train <- scaled[train_i, ]
  test <- scaled[-train_i, ]

  net <- neuralnet(f, data = train, hidden = c(5, 3), linear.output = TRUE)
  pred_net <- predict(net, test[, setdiff(names(test), "medv")])
  pred_lm <- predict(glm(medv ~ ., data = train), test)

  tibble(
    seed = seed,
    net = caret::RMSE(pred_net, test$medv),
    lm = caret::RMSE(pred_lm, test$medv)
  )
}

map(1:5, compare) %>% bind_rows()
```

With 13 predictors and 506 rows, this is not where neural networks shine. Their advantage appears when the input is high-dimensional and the structure is genuinely non-linear, as with images, text or audio. Try changing `hidden` (more neurons, more layers) and see whether the test error improves or whether the network starts to overfit.