```{r, include = FALSE}
source("../setup.R")
library(ggpubr)
library(kableExtra)
```

## Overview

We will cover:

* Structure of a neural network
* Fitting neural networks
* Diagnosing the fit

## Structure of a neural network {.transition-slide .center}

## Nested logistic regressions

:::: {.columns}
::: {.column}

Remember the logistic function:

\begin{align}
f(x) &= \frac{e^{\beta_0+\sum_{j=1}^p\beta_jx_j}}{1+e^{\beta_0+\sum_{j=1}^p\beta_jx_j}}\\
  &= \frac{1}{1+e^{-(\beta_0+\sum_{j=1}^p\beta_jx_j)}}
\end{align}

Also,

$$\log_e\frac{f(x)}{1 - f(x)} = \beta_0+\sum_{j=1}^p\beta_jx_j$$
:::

::: {.column}
::: {.fragment}

<br><br>
[Above the threshold predict to be 1.]{.smaller} 

```{r}
#| echo: false
#| out-width: 70%
#| fig-width: 4
#| fig-height: 4
library(tidyverse)
x <- seq(-2, 2, 0.1)
y <- exp(1+3*x)/(1+exp(1+3*x))
df2 <- tibble(x, y)
ggplot(df2, aes(x=x, y=y)) + 
  geom_line() +
  geom_hline(yintercept=0.5, colour="orange") +
  annotate("text", x=0.84, y=0.55, label="Activation threshold ??", colour="orange") +
  geom_hline(yintercept=c(0,1), linetype=2)
```

:::

:::
::::

## Linear regression as a network

:::: {.columns}
::: {.column}
$$\widehat{y} =b_0+\sum_{j=1}^pb_jx_j$$

Drawing as a network model: 

$p$ [inputs]{.monash-orange2} (predictors), multiplied by [weights]{.monash-orange2} (coefficients), summed, add a [constant]{.monash-orange2}, predicts [output]{.monash-orange2} (response). 
:::
::: {.column}
![](../images/reg_nn.png){width=80%}

:::
::::

## Single hidden layer NN

:::: {.columns}
::: {.column}

\begin{align}
\widehat{y} =a_{0}+\sum_{k=1}^s(a_{k}(b_{0k}+\sum_{j=1}^pb_{jk}x_j))
\end{align}


:::
::: {.column}

![](../images/nn_annotate.png)
:::
::::

## What does this look like? [(1/2)]{.smallest}

:::: {.columns}
::: {.column style="font-size: 70%"}
The architecture allows for combining multiple linear models to generate non-linear classifications. 

```{r}
#| echo: false
#| out-width: 70%
#| fig-width: 4
#| fig-height: 5
w <- read_csv(here::here("data/wiggly.csv"))
ggplot(w, aes(x=x, y=y, colour=class, shape=class)) + 
  geom_point() +
  scale_color_brewer("", palette="Dark2") +
  scale_shape("") +
  theme(legend.position = "bottom") 
```

:::
::: {.column style="font-size: 70%"}

The best fit uses $s=4$, four nodes in the hidden layer. Can you sketch four lines that would split this data well?

```{r}
#| echo: false
#| out-width: 70%
#| fig-width: 4
#| fig-height: 5
load(here::here("data/nnet_many.rda"))
load(here::here("data/nnet_best.rda"))

ggplot(subset(best$output,  node == 1), aes(x, y)) +
  geom_raster(aes(fill = pred)) +
  geom_point(aes(shape = class), data = w) +
  scale_fill_gradient2("", low="#1B9E77", 
                       high="#D95F02", 
                       mid = "white", 
                       midpoint = 0.5,
                       guide = "colourbar",
                       limits = c(0,1)) +
  scale_shape("") +
  theme(legend.position = "bottom",
        legend.text = element_text(size=6)) 
```
:::
::::

[[Wickham et al (2015) Removing the Blindfold](http://onlinelibrary.wiley.com/doi/10.1002/sam.11271/abstract)]{.smallest}

## What does this look like? [(2/2)]{.smallest}

:::: {.columns}
::: {.column}

The models at each of the nodes of the hidden layer. 

```{r}
#| echo: false
#| out-width: 90%
#| fig-width: 5
#| fig-height: 4
ggplot(data=best$hidden, aes(x, y)) +
  geom_tile(aes(fill=pred)) + 
  geom_point(data=w, aes(shape=class)) +
  facet_wrap(~node, ncol=2) + 
  scale_fill_gradient2("", low="#AF8DC3",
                                    mid="#F7F7F7",
                                    high="#7FBF7B",
                                    midpoint=0.5,
                                    limits=c(0,1)) +
  scale_shape("") +
  theme(axis.text = element_blank(), 
        axis.title = element_blank(),
        legend.text = element_text(size=6))
```
:::

::: {.column}
::: {.fragment}

```{r}
#| echo: false
#| out-width: 70%
#| fig-width: 4
#| fig-height: 4
ggplot(data=best$hidden, aes(x, y)) + 
  geom_contour(aes(z=pred, group=node), 
               colour="grey50", 
               size=2, 
               breaks = 0.5) +
  geom_point(data=w, aes(colour=class, 
                 shape=class)) + 
  scale_color_brewer("", palette="Dark2") +
  scale_shape("") +
  theme(legend.position = "bottom",
        legend.text = element_text(size=6)) 


```
:::
:::
::::

## But can be painful to find the best!

These are all the models fitted, using $s=2, 3, 4$ with the fit statistics.

```{r}
#| echo: false
#| out-width: 100%
#| fig-width: 12
#| fig-height: 4
qual <- unique(many[, c("value", "accuracy", "nodes", "id")])
ggplot(qual, aes(x=accuracy, y=value)) +
  geom_point(alpha=0.7, size=3) +
  xlab("Accuracy") +
  ylab("Value of fitting criterion") +
  facet_wrap(~nodes)

```

Fitted using the R package `nnet`. It's very unstable, and this is still a problem with current procedures.

## Fitting with keras {.transition-slide .center}

## Steps [(1/2)]{.smaller}

:::: {.columns}
::: {.column}

1. Define [architecture]{.monash-blue2}

    - flatten: if you are classifying images, you need to flatten the image into a single row of data, eg 24x24 pixel image would be converted to a row of 576 values. Each pixel is a variable.
    - How many hidden layers do you need? 
    - How many nodes in the hidden layer?
    - Dropout rate: proportion of nodes removed randomly at each update, for regularisation, to reduce number of parameters to be estimated 

:::
::: {.column}
::: {.fragment}
2. Specify [activation]{.monash-blue2}: linear, relu (rectified linear unit), sigmoid, softmax
:::

::: {.fragment}
3. Choose [loss]{.monash-blue2} function: 

    - MSE: differencing predictive probabilities from binary matrix specified response, eg predict=(0.91, 0.07, 0.02) and true=(1,0,0) then loss is (1-0.91)^2=0.0081.
    - cross-entropy: $\sum -p(x)log_e(q(x))$ where $p$ is true, and $q$ is predicted, eg -1 log_e(0.91)=0.094
:::

:::
::::

## Steps [(2/2)]{.smaller}

:::: {.columns}
::: {.column}
4. [Training]{.monash-blue2} process:

    - epochs: number of times the algorithm sees the entire data set
    - batch_size: subset used at each fit
    - validation_split: proportion for hold-out set for computing error rate
    - batch_normalization: each batch is standardised during the fitting, can be helpful even if full data is standardised



:::
::: {.column}
::: {.fragment}
5. Evaluation: 

    - usual metrics: [accuracy]{.monash-blue2}, ROC, AUC (area under ROC curve), confusion table 
    - [visualise]{.monash-blue2} the predictive probabilities
    - examine misclassifications
    - examine specific nodes, to understand what part it plays
    - examine model boundary relative to the observed data
:::
:::
::::

## Example: penguins [(1/5)]{.smaller}

```{r echo=FALSE}
#| message: false
#| eval: false
animate_xy(p_std[,2:5], col=factor(p_tidy$species))
render_gif(p_std[,2:5], 
           grand_tour(), 
           display_xy(half_range=4.2, axes="off",
                      col=p_std$species),
           gif_file="gifs/penguins_lbl.gif",
           frames=500,
           width = 400,
           height = 400,
           loop=FALSE)
```

:::: {.columns}
::: {.column}

<center>
![](../gifs/penguins_lbl.gif)
</center>

- 4D data
- Simple cluster structure to classes
- How many nodes in the hidden layer?

:::
::: {.column}

::: {.fragment}
Choose 2 nodes, because reducing to 2D, like LDA discriminant space, makes for easy classification.

![](https://dicook.github.io/mulgar_book/images/nn-diagram.png)
:::

:::
::::

## Example: penguins [(2/5)]{.smaller}

:::: {.columns}
::: {.column}

```{r eval=FALSE}
#| echo: true
library(keras)
tensorflow::set_random_seed(211)

# Define model
p_nn_model <- keras_model_sequential()
p_nn_model %>% 
  layer_dense(units = 2, activation = 'relu', 
              input_shape = 4) %>% 
  layer_dense(units = 3, activation = 'softmax')
p_nn_model %>% summary

loss_fn <- loss_sparse_categorical_crossentropy(
  from_logits = TRUE)

p_nn_model %>% compile(
  optimizer = "adam",
  loss      = loss_fn,
  metrics   = c('accuracy')
)
```

Note that the `tidymodels` code style does not allow easy extraction of model coefficients. 

:::
::: {.column}

Split the data into training and test, and check it.

```{r echo=FALSE}
set.seed(821)
p_split <- p_std %>% 
  initial_split(prop = 2/3, 
                strata=species)
p_train <- training(p_split)
p_test <- testing(p_split)

# Check training and test split
p_split_check <- bind_rows(
  bind_cols(p_train, type = "train"), 
  bind_cols(p_test, type = "test")) %>%
  mutate(type = factor(type))
```

```{r echo=FALSE, eval=FALSE}
render_gif(p_split_check[,2:5],
           guided_tour(lda_pp(p_split_check$species)),
           display_xy( 
             col=p_split_check$species, 
             pch=p_split_check$type, 
             shapeset=c(16,1),
             cex=1.5,
             axes="bottomleft"), 
           gif_file="gifs/p_split_guided.gif",
           frames=500,
           loop=FALSE
)
```

<center>
![](../gifs/p_split_guided.gif)
</center>

:::
::::

## Example: penguins [(3/5)]{.smaller}

:::: {.columns}
::: {.column}
Fit the model

```{r echo=TRUE}
# Data needs to be matrix, and response needs to be numeric
p_train_x <- p_train %>%
  select(bl:bm) %>%
  as.matrix()
p_train_y <- p_train %>% pull(species) %>% as.numeric() 
p_train_y <- p_train_y-1 # Needs to be 0, 1, 2
p_test_x <- p_test %>%
  select(bl:bm) %>%
  as.matrix()
p_test_y <- p_test %>% pull(species) %>% as.numeric() 
p_test_y <- p_test_y-1 # Needs to be 0, 1, 2
```

```{r echo=TRUE, eval=FALSE}
#| message: false
# Fit model
p_nn_fit <- p_nn_model %>% 
  keras::fit(
    x = p_train_x, 
    y = p_train_y,
    epochs = 200,
    verbose = 0
  )
```


:::

::: {.column}

How many parameters need to be estimated?

::: {.fragment style="font-size: 60%;"}
Four input variables, two nodes in the hidden layer and a three column binary matrix for output. This corresponds to 5+5+3+3+3=19 parameters.
:::

::: {.fragment}
<br><br>
```{r echo=FALSE}
#| message: false
p_nn_model <- load_model_tf("../data/penguins_cnn")
p_nn_model
```
:::

:::

::::

## Example: penguins [(4/5)]{.smaller}

:::: {.columns}
::: {.column}

Evaluate the fit

```{r echo=TRUE}
p_nn_model %>% 
  evaluate(p_test_x, p_test_y, verbose = 0)
```


Confusion matrices for training and test

```{r echo=FALSE}
# Predict training and test set
p_train_pred <- p_nn_model %>% 
  predict(p_train_x, verbose = 0)
p_train_pred_cat <- levels(p_train$species)[
  apply(p_train_pred, 1,
        which.max)]
p_train_pred_cat <- factor(
  p_train_pred_cat,
  levels=levels(p_train$species))
table(p_train$species, p_train_pred_cat)

p_test_pred <- p_nn_model %>% 
  predict(p_test_x, verbose = 0)
p_test_pred_cat <- levels(p_test$species)[
  apply(p_test_pred, 1, 
        which.max)]
p_test_pred_cat <- factor(
  p_test_pred_cat,
  levels=levels(p_test$species))
table(p_test$species, p_test_pred_cat)
```
[Note: Specifically have chosen settings so fit is not perfect]{.smallest}

:::
::: {.column}

::: {.fragment}
Estimated parameters

```{r}
# Extract hidden layer model weights
p_nn_wgts <- keras::get_weights(p_nn_model, trainable=TRUE)
p_nn_wgts 
```

[Which variables are contributing most to each hidden layer node?]{.smaller}

[Can you write out the model?]{.smaller}

:::

:::
::::

## Example: penguins [(5/5)]{.smaller}

:::: {.columns}
::: {.column}

Check the fit at the hidden layer nodes

```{r echo=FALSE}
# Orthonormalise the weights to make 2D projection
p_nn_wgts_on <- tourr::orthonormalise(p_nn_wgts[[1]])

# Hidden layer
p_train_m <- p_train %>%
  mutate(nn1 = as.matrix(p_train[,2:5]) %*%
           as.matrix(p_nn_wgts_on[,1], ncol=1),
         nn2 = as.matrix(p_train[,2:5]) %*%
           matrix(p_nn_wgts_on[,2], ncol=1))

# Now add the test points on.
p_test_m <- p_test %>%
  mutate(nn1 = as.matrix(p_test[,2:5]) %*%
           as.matrix(p_nn_wgts_on[,1], ncol=1),
         nn2 = as.matrix(p_test[,2:5]) %*%
           matrix(p_nn_wgts_on[,2], ncol=1))
p_train_m <- p_train_m %>%
  mutate(set = "train")
p_test_m <- p_test_m %>%
  mutate(set = "test")
p_all_m <- bind_rows(p_train_m, p_test_m)
ggplot(p_all_m, aes(x=nn1, y=nn2, 
                     colour=species, shape=set)) + 
  geom_point() +
  scale_colour_discrete_divergingx(palette="Zissou 1") +
  scale_shape_manual(values=c(16, 1)) +
  theme_minimal() +
  theme(aspect.ratio=1)
```

::: {.fragment  style="font-size: 70%"}
This is the dimension reduction induced by the model.
:::

::: {.fragment  style="font-size: 50%"}
[Realistically, with a complex neural network, it is too much work to check these nodes.]{.monash-blue2}
:::

:::
::: {.column}

::: {.fragment}
Examine the predictive probabilities

```{r echo=FALSE}
# Set up the data to make the ternary diagram
# Join data sets
colnames(p_train_pred) <- c("Adelie", "Chinstrap", "Gentoo")
colnames(p_test_pred) <- c("Adelie", "Chinstrap", "Gentoo")
p_train_pred <- as_tibble(p_train_pred)
p_train_m <- p_train_m %>%
  mutate(pspecies = p_train_pred_cat) %>%
  bind_cols(p_train_pred) %>%
  mutate(set = "train")
p_test_pred <- as_tibble(p_test_pred)
p_test_m <- p_test_m %>%
  mutate(pspecies = p_test_pred_cat) %>%
  bind_cols(p_test_pred) %>%
  mutate(set = "test")
p_all_m <- bind_rows(p_train_m, p_test_m)

# Add simplex to make ternary
library(geozoo)
proj <- t(geozoo::f_helmert(3)[-1,])
p_nn_v_p <- as.matrix(p_all_m[,c("Adelie", "Chinstrap", "Gentoo")]) %*% proj
colnames(p_nn_v_p) <- c("x1", "x2")
p_nn_v_p <- p_nn_v_p %>%
  as.data.frame() %>%
  mutate(species = p_all_m$species,
         set = p_all_m$set)

simp <- geozoo::simplex(p=2)
sp <- data.frame(cbind(simp$points), simp$points[c(2,3,1),])
colnames(sp) <- c("x1", "x2", "x3", "x4")
sp$species = sort(unique(p_std$species))
ggplot() +
  geom_segment(data=sp, aes(x=x1, y=x2, xend=x3, yend=x4)) +
  geom_text(data=sp, aes(x=x1, y=x2, label=species),
            nudge_x=c(-0.1, 0.15, 0),
            nudge_y=c(0.05, 0.05, -0.05)) +
  geom_point(data=p_nn_v_p, aes(x=x1, y=x2, 
                                colour=species,
                                shape=set), 
             size=2, alpha=0.5) +
  scale_color_discrete_divergingx(palette="Zissou 1") +
  scale_shape_manual(values=c(19, 1)) +
  theme_map() +
  theme(aspect.ratio=1, legend.position = "right")
```
:::

::: {.fragment  style="font-size: 50%"}
[The problem with this model is that Gentoo are too confused with Adelie. This is a structural problem because from the visualisation of the 4D data we know that there is a big gap between the Gentoo and both other species.]{.monash-blue2}
:::

:::
::::

## Want to learn more?

- Work your way through the example of fitting the [fashion MNIST data](https://tensorflow.rstudio.com/tutorials/keras/classification) using tensorflow.

- [Hands-on machine learning](https://bradleyboehmke.github.io/HOML/deep-learning.html) has a lovely step-by-step guide to constructing and fitting. 

- This is a very nice slide set: [A gentle introduction to deep learning in R using Keras](https://lnalborczyk.github.io/slides/vendredi_quanti_2021/vendredi_quantis#1)

- And the tutorials at [TensorFlow for R](https://tensorflow.rstudio.com/install/) have lots of examples.

## Next: Explainable artificial intelligence (XAI) {.transition-slide .center}


