```{r echo=FALSE}
# Set up chunk for all slides
knitr::opts_chunk$set(
  fig.width = 4,
  fig.height = 4,
  fig.align = "center",
  out.width = "60%",
  code.line.numbers = FALSE,
  fig.retina = 3,
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  cache = FALSE,
  dev.args = list(pointsize = 11)
)
```

```{r}
#| echo: true
#| code-fold: true
#| code-summary: "Load the libraries and avoid conflicts"
# Load libraries used everywhere
library(tidyverse)
library(tidymodels)
library(patchwork)
library(mulgar)
library(palmerpenguins)
library(GGally)
library(tourr)
library(MASS)
library(discrim)
library(classifly)
library(detourr)
library(crosstalk)
library(plotly)
library(viridis)
library(colorspace)
library(conflicted)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::slice)
conflicts_prefer(palmerpenguins::penguins)
conflicts_prefer(viridis::viridis_pal)

options(digits=2)
p_tidy <- penguins |>
  select(species, bill_length_mm:body_mass_g) |>
  rename(bl=bill_length_mm,
         bd=bill_depth_mm,
         fl=flipper_length_mm,
         bm=body_mass_g) |>
  filter(!is.na(bl)) |>
  arrange(species) |>
  na.omit()
p_tidy_std <- p_tidy |>
    mutate_if(is.numeric, function(x) (x-mean(x))/sd(x))
```

```{r}
#| echo: false
# Set plot theme
theme_set(theme_bw(base_size = 14) +
   theme(
     aspect.ratio = 1,
     plot.background = element_rect(fill = 'transparent', colour = NA),
     plot.title.position = "plot",
     plot.title = element_text(size = 24),
     panel.background = element_rect(fill = 'transparent', colour = NA),
     legend.background = element_rect(fill = 'transparent', colour = NA),
     legend.key = element_rect(fill = 'transparent', colour = NA)
   )
)
```

## `r emo::ji("target")` Objectives

The goal for this week is learn to fit, diagnose, assess assumptions, and predict from logistic regression models, and linear discriminant analysis models. 

## `r emo::ji("wrench")` Preparation 

- Make sure you have all the necessary libraries installed. There are a few new ones this week!

## Exercises: 

Open your project for this unit called `iml.Rproj`. For all the work we will use the penguins data. Start with splitting it into a training and test set, as follows. 

```{r}
set.seed(1148)
p_split <- initial_split(p_tidy_std, 2/3, strata = species)
p_tr <- training(p_split)
p_ts <- testing(p_split)
```

#### 1. LDA

This problem uses linear discriminant analysis on the penguins data. 

a. Is the assumption of equal variance-covariance reasonable to make for this data? 

::: unilur-solution
You need to look at the data in a tour, using:

```{r}
#| eval: false
animate_xy(p_tidy_std[,2:5], col=p_tidy$species)
```

Use the standardised data, because the measurements are in different sizes, and this is not relevant for this data. 
:::

b. Fit the LDA model to the training data, using this code

```{r}
lda_spec <- discrim_linear() |>
  set_mode("classification") |>
  set_engine("MASS", prior = c(1/3, 1/3, 1/3))
lda_fit <- lda_spec |> 
  fit(species ~ ., data = p_tr)
```

::: unilur-solution

```{r}
#| echo: false
lda_fit
```
:::

c. Compute the confusion matrices for training and test sets, and thus the error for the test set. 

::: unilur-solution
```{r}
#| echo: false
p_tr_pred <- p_tr |>
  mutate(pspecies = predict(lda_fit$fit, p_tr)$class)
p_ts_pred <- p_ts |>
  mutate(pspecies = predict(lda_fit$fit, p_ts)$class)

p_tr_pred |> count(species, pspecies) |>
  group_by(species) |>
  mutate(cl_acc = n[pspecies==species]/sum(n)) |>
  pivot_wider(names_from = pspecies, 
              values_from = n, values_fill=0) |>
  select(species, Adelie, Chinstrap, Gentoo, cl_acc)

p_ts_pred |> count(species, pspecies) |>
  group_by(species) |>
  mutate(cl_acc = n[pspecies==species]/sum(n)) |>
  pivot_wider(names_from = pspecies, 
              values_from = n, values_fill=0) |>
  select(species, Adelie, Chinstrap, Gentoo, cl_acc)

accuracy(p_ts_pred, species, pspecies)$.estimate
```
:::

d. Plot the training and test data in the discriminant space, using symbols to indicate which set. See if you can mark the misclassified cases, too.

::: unilur-solution

```{r}
#| echo: false
p_pred <- predict(lda_fit$fit, p_tidy_std)$x |>
  as.data.frame()
p_pred <- bind_cols(p_tidy_std, p_pred) 
p_pred <- p_pred |>
  mutate(pspecies = predict(lda_fit$fit, p_tidy_std)$class) |>
  mutate(set = ifelse(1:nrow(p_tidy_std) %in% p_split$in_id, 
                      "train", "test"),
         err = ifelse(species != pspecies, "yes", "no"))

ggplot(p_pred, aes(x=LD1, y=LD2, 
                   colour=species, 
                   shape=set)) +
  geom_point() +
  geom_point(data=p_pred[p_pred$err == "yes",], shape=1, cex=5) +
  scale_color_discrete_divergingx("Zissou 1") 
```
:::

e. Re-do the plot of the discriminant space, to examine the boundary between groups. You'll need to generate a set of random points in the domain of the data, predict their class, and projection into the discriminant space. The `explore()` in the `classifly` package can help you generate the box of random points.

::: unilur-solution

```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 4
p_bnd <- explore(lda_fit$fit, p_tidy_std)
p_bnd_ds <- predict(lda_fit$fit, p_bnd, dimen=2)$x 
p_bnd_ds <- as.data.frame(p_bnd_ds)
p_bnd <- bind_cols(p_bnd, p_bnd_ds)
ggplot(p_bnd) +
  geom_point(aes(x=LD1, y=LD2, 
                 colour=species, 
                 shape=.TYPE), alpha=0.8) + 
  scale_color_discrete_divergingx("Zissou 1") +
  scale_shape_manual(values=c(46, 16)) 
```

:::

f. What happens to the boundary, if you change the prior probabilities? And why does this happen? Change the prior probabilities to be 1.999/3, 0.001/3, 1/3 for Adelie, Chinstrap, Gentoo, respectively. Re-do the plot of the boundaries in the discriminant space. 

::: unilur-solution
```{r}
#| echo: false
#| fig-width: 8
#| fig-height: 4
lda_spec2 <- discrim_linear() |>
  set_mode("classification") |>
  set_engine("MASS", prior = c(1.999/3, 0.001/3, 1/3))
lda_fit2 <- lda_spec2 |> 
  fit(species ~ ., data = p_tr)
p_bnd2 <- explore(lda_fit2$fit, p_tidy_std)
p_bnd2_ds <- predict(lda_fit2$fit, p_bnd2, dimen=2)$x 
p_bnd2_ds <- as.data.frame(p_bnd2_ds)
p_bnd2 <- bind_cols(p_bnd2, p_bnd2_ds)
ggplot(p_bnd2) +
  geom_point(aes(x=LD1, y=LD2, 
                 colour=species, 
                 shape=.TYPE), alpha=0.8) + 
  scale_color_discrete_divergingx("Zissou 1") +
  scale_shape_manual(values=c(46, 16)) 
```

If the prior probabilities are unequal, it gives more importance to some classes. Here the importance of the Adelie penguins has been increased to the detriment of the Chinstrap. So the boundary moves away from the Adelie, which means more often a new penguin would be classified as an Adelie. 
:::

#### 2. Logistic

a. Fit a logistic discriminant model to the training set. You can use this code:

```{r}
log_fit <- multinom_reg() |> 
  fit(species ~ ., 
      data = p_tr)
```

::: unilur-solution
```{r}
log_fit
```
:::

b. Compute the confusion matrices for training and test sets, and thus the error for the test set. You can use this code to make the predictions.

```{r}
p_tr_pred <- log_fit |> 
  augment(new_data = p_tr) |>
  rename(pspecies = .pred_class)
p_ts_pred <- log_fit |> 
  augment(new_data = p_ts) |>
  rename(pspecies = .pred_class)
```

::: unilur-solution
```{r}
p_tr_pred |> count(species, pspecies) |>
  group_by(species) |>
  mutate(cl_acc = n[pspecies==species]/sum(n)) |>
  pivot_wider(names_from = pspecies, 
              values_from = n, values_fill=0) |>
  select(species, Adelie, Chinstrap, Gentoo, cl_acc)

p_ts_pred |> count(species, pspecies) |>
  group_by(species) |>
  mutate(cl_acc = n[pspecies==species]/sum(n)) |>
  pivot_wider(names_from = pspecies, 
              values_from = n, values_fill=0) |>
  select(species, Adelie, Chinstrap, Gentoo, cl_acc)

accuracy(p_ts_pred, species, pspecies)$.estimate
```
:::

c. Check the boundaries produced by logistic regression, and how they differ from those of LDA. Using the 2D projection produced by the LDA rule (using equal priors) predict the your set of random points using the logistic model. 

::: unilur-solution
```{r}
p_log_bnd_ds <- log_fit |> 
  augment(new_data = p_bnd) |>
  rename(pspecies = .pred_class)

ggplot() +
  geom_point(
    data=p_log_bnd_ds[p_log_bnd_ds$.TYPE == "simulated",], 
             aes(x=LD1, y=LD2, 
                 colour=pspecies), shape=46, alpha=0.8) + 
  scale_color_discrete_divergingx("Zissou 1") +
  geom_point(data=p_log_bnd_ds[p_log_bnd_ds$.TYPE == "actual",],
             aes(x=LD1, y=LD2, 
                 colour=species), shape=16, alpha=0.8) 
```

One thing that you can notice is that the boundaries are not "crisp", that there is overlap of the coloured points marking the classification regions. This means that the separation from the logistic regression model is not accomplished in the same 2D space as LDA. 
:::
  
#### 3. Interactively explore misclassifications 

Here you are going to use interactive graphics to explore the misclassifications from the linear discriminant analysis. We'll need to use detourr to accomplish this. The code below makes a scatterplot of the confusion matrix, where points corresponding to a class have been spread apart by jittering. This plot is linked to a tour plot. Try:

1. Selecting penguins that have been misclassified, from the display of the confusion matrix. Observe where they are in the data space. Are they in an area where it is hard to distinguish the groups?
2. Selecting neighbouring points in the tour, and examine where they are in the confusion matrix. 

```{r}
#| echo: true
#| eval: false
p_cl <- p_tidy_std |>
  mutate(pspecies = predict(lda_fit$fit, p_tidy_std)$class) |>
  dplyr::select(bl:bm, species, pspecies) |>
  mutate(sp_jit = jitter(as.numeric(species)),
         psp_jit = jitter(as.numeric(pspecies)))
p_cl_shared <- SharedData$new(p_cl)

detour_plot <- detour(p_cl_shared, tour_aes(
  projection = bl:bm,
  colour = species)) |>
    tour_path(grand_tour(2), 
                    max_bases=50, fps = 60) |>
       show_scatter(alpha = 0.9, axes = FALSE,
                    width = "100%", height = "450px")

conf_mat <- plot_ly(p_cl_shared, 
                    x = ~psp_jit,
                    y = ~sp_jit,
                    color = ~species,
                    colors = viridis_pal(option = "D")(3),
                    height = 450) |>
  highlight(on = "plotly_selected", 
              off = "plotly_doubleclick") %>%
    add_trace(type = "scatter", 
              mode = "markers")
  
bscols(
     detour_plot, conf_mat,
     widths = c(5, 6)
 )            
```

#### 4. Exploring the math

Slide 23 of the lecture notes has the steps to go from Bayes rule to the discriminant functions. Explain what was done at each step to get to the next one. 

## `r emo::ji("wave")` Finishing up

Make sure you say thanks and good-bye to your tutor. This is a time to also report what you enjoyed and what you found difficult.
