```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidymodels)
library(janitor)
```

# Background

Reaven and Miller (1979) examined the relationship among blood chemistry measures of glucose tolerance and insulin in 145 non-obese adults. They used the PRIM9 system at the Stanford Linear Accelerator Center to visualize the data in 3D, and discovered a peculiar pattern that looked like a large blob with two wings in different directions.

After further analysis, the subjects were classified as sub-clinical (chemical) diabetics, overt diabetics and normals. This study was influential in defining the stages of development of Type 2 diabetes. Overt diabetes is the most advanced stage, characterized by elevated fasting blood glucose concentration and classical symptoms. Preceding overt diabetes is the latent or chemical diabetic stage, with no symptoms of diabetes but demonstrable abnormality of oral or intravenous glucose tolerance.

[Original paper](https://link.springer.com/content/pdf/10.1007/BF00423145.pdf)

# EDA

First, notice that the classes are not evenly distributed. What does this mean for us when we split?

There are some high correlations between some variables and they are on very different scales. In terms of insulin and glucose, the Overt group has a lot more variability.

```{r}
# The data is included in the heplots library
library(heplots)
data("Diabetes")
#help(Diabetes)

# numeric summaries
library(skimr)
skim(Diabetes)
cor(Diabetes %>% select_if(is.numeric))

# visual summaries
# ggpairs(diabetes, aes(color = group))
```

# Splits and Recipe

```{r}
set.seed(2023)
diab_split <- Diabetes %>%
  initial_split(
    prop = 0.7,
    strata = group
  )
diab_test <- testing(diab_split)
diab_train <- training(diab_split)

recipe <- recipe(group ~ ., data = Diabetes) %>%
  step_normalize(all_predictors())
```

# LDA and QDA

Let's learn to specify a linear and quadratic discriminant analysis model. Make sure you install and load the `parsnip` extension `discrim`, which implements these discriminant methods.

```{r}
library(discrim)
lda_model <- discrim_linear() %>%
  # The default engine is MASS
  # Other engines allow for penalization/regularization
  set_engine("MASS") %>%
  set_mode("classification")

qda_model <- discrim_quad() %>%
  set_engine("MASS") %>%
  set_mode("classification")
```

Let's put these together into a single `workflowset`, along with a KNN model. Notice how each model is identified by a combination of the recipe name and model name that we give in the list.

```{r}
diab_wfs <- workflow_set(
  preproc = list(rec = recipe),
  models = list(
    knn = nearest_neighbor(neighbors=9) %>%
      set_mode("classification") %>%
      set_engine("kknn"),
    lda = lda_model,
    qda = qda_model,
    nb = naive_Bayes() %>%
      set_engine("klaR") %>%
      set_mode("classification")
    # logistic regression had problems with this example,
    # so I'm leaving it out!
    #logreg = logistic_reg(engine = "glm")
  )
)

diab_wfs
```

Let's use cross-validation to see which model performs the best on the training data:

```{r}
folds <- vfold_cv(diab_train, v=5, repeats=5, strata=group)

diab_wfs <- diab_wfs %>% workflow_map(
  # the name of the function to map:
  "fit_resamples", 
  # Options to `fit_resamples()`: 
  resamples = folds, metrics = metric_set(sens, spec, accuracy, roc_auc, f_meas)
  )

# The result variable in the workflowset is now filled
diab_wfs

collect_metrics(diab_wfs)
```
There is also an  `autoplot` method that will plot the results of the resampling. Take a look at this plot and make sure you understand what you are looking at!

```{r}
autoplot(diab_wfs)
```


Notice that the standard errors are higher than our other examples, since each of the folds is somewhat small!

```{r}
diab_results <- diab_wfs %>% 
  extract_workflow("rec_nb") %>%
  fit(data=diab_train) %>%
  augment(new_data=diab_test)
```

<!--
```{r}
conf_mat(diab_results, truth=group, estimate=.pred_class)
sens(diab_results, truth=group, estimate=.pred_class)
kap(diab_results, truth=group, estimate=.pred_class)
f_meas(diab_results, truth=group, estimate=.pred_class)
roc_auc(diab_results, truth=group, c(.pred_Normal, .pred_Chemical_Diabetic,.pred_Overt_Diabetic))
```
-->


# Group project

Your project for next class is to teach us about multiclass classification metrics. Now that you are familiar with this example, you can use it however you would like to demonstrate your group's metric. 

Some questions to consider: Can you give the intuition for what your metric is measuring? What are the values of the metric for these problems? Is there a baseline for what constitutes a "good" value of this metric?  In the context of this problem (classifying diabetes patients), is your metric more or less important than accuracy (or another metric we've learned)?