---
title: "Unit 05 Lab Agenda"
date: "02/24/2025"
format: 
  html: 
    embed-resources: true
    toc: true 
    toc_depth: 4
editor_options: 
  chunk_output_type: console
---


### 1. Review Key

#### Discrepancy between data dictionary and data and possible ways to address.  

- First, a careful descriptive peek into the data helps avoid extra error later. We'd better make sure if the data is distributed reasonably as the dictionary shows.

- How do we handle the decimals in "nominal" variables and do we need to make additional transformations since the data are read in as numeric and it is already ordered 0-4? 
  - We could tackle with this in different ways. In the key, we adopted the method to treat non-integer values as NAs. Another way is to round the decimals into integers. 
  - When we round the values and plan to treat them as numeric, we're assuming equal intervals between levels. However, if the data is not distributed very linearly, deeming those as numeric interval variables may lead to the overlook of comparable relationship between levels.
  - Turning the variables as factors could help us generate dummy features and we could utilize the paired up between level relationships.

#### If we had used the forcats package we would have known NA values were introduced!    

Lets look at our data again
```{r}
library(tidyverse) |> suppressPackageStartupMessages()

path_data <- "application_assignments/unit_05"

data_all <- read_csv(here::here(path_data, "smoking_ema.csv"),
                     show_col_types = FALSE) |> 
  glimpse()
```

Now as a reminder when we classed our factor levels with `factor()` we didn't get any warnings or errors. Lets now class using `forcats::fct()`    

I'm only going to use `craving` for this example   

Also note a few differences between `factor()` and `fct()`    

- `factor()` codes levels alphabetically (unless otherwise specified)
- `fct()` codes levels by first appearance in data (unless otherwise specified)
- If your variable is numeric (like in this case) you must convert it to character 
- With `fct()` you cannot add labels that differ from levels in same command, you must follow with `fct_recode()`

```{r}
#| eval: false
data_all <- data_all |> 
  mutate(craving = fct(as.character(craving), levels = c("0", "1", "2", "3", "4"))) |> 
  glimpse()
```

Here we got errors showing us the unexpected values in our factor!



### 2. Parrallel processing and caching

#### Parallel processing

- Is there any way(such as CUDA) to speed up R Studio? Running the KNN takes a long time.
  - This is beyond the scope of our course. More advanced compute hardware and architectures are more used for neural network based models.

This is the code we use to set up parallel processing
```{r}
cl <- parallel::makePSOCKcluster(parallel::detectCores(logical = FALSE))
doParallel::registerDoParallel(cl)
```

To see how many cores you have you can look at `cl` or type `parallel::detectCores(logical = FALSE)` in your console.

```{r}
cl
```

#### Cache

When we're doing caching, we store the results of the processing steps as objects. And then, we could read them in directly once we run through the process. It could save us a lot of time when working with large scale data and repetitive procedures.

If you want to set up the rerun option globally, assign it to a variable.

```{r}
library(xfun, include.only = "cache_rds")
rerun_setting <- F
```

```{r}
results <- cache_rds(
  expr = {
    Sys.sleep(10)
    results <- c(1, 2, 3)
  }, 
  dir = "cache/",
  file = "cache_demo_xfun",
  rerun = rerun_setting)
```

```{r}
results
```


After running the first time, the rerun will be quicker because the process had been cached in the folder.

And we could also do the cache more manually with `write_rds` ourselves.

```{r}
if (file.exists(here::here("lab/cache/brute_force_method.rds"))){
  
 # file with computations exists, so just load it
 results <- readr::read_rds(here::here("lab/cache/brute_force_method.rds"))
  
} else {
 
  # pretend these computations take a while by calling Sys.sleep ()
  Sys.sleep(10)
  results <- 1 
 
  # write results so we have it the next time we run this code chunk 
  results |>  readr::write_rds(here::here("lab/cache/brute_force_method.rds"))
}

results
```

This way, the cached file could also be uploaded and downloaded for the usage else where (instead of being limited locally with hashing suffixes by `xfun`). More information see John's tutorial: https://jjcurtin.github.io/book_dwvt/cache.html.

### 3. Hyperparameter tuning

#### KNN

(Review the key)

- Can you tell us what to do when our performance estimates are taking too long to load? It's frustrating to have to run it multiple times and wait for a while only to find out I've made an error.

- Also any advice on setting a hyperparameter grid to ensure that we don't have to revisit it to capture a wide enough range of k values?

- When we use KNN, what range do we typically choose?

  - I will suggest a wider range of hyperparameter searching with small but reasonable number of CV sets in the beginning, and smaller range after the first location of optimal range. Especially, pay attention to whether a sweet spot appears. Also a small number of CV sets will help debug before a final recipe & processing steps are ready.
  
- Also, my hyperparameter grid plots looked pretty whacky. There was not much movement in the estimates for changes in k. It could be an error in my code but if we could talk about why changing k should or should not change estimates/SE that would be helpful.
  - When changing the k from a very small value to larger, the performance is usually expected to increase, on account of the reduction in variance. And it gradually reaches an optimal value, where variance and bias are counterbalanced best for the training data. When the k value becomes larger and larger, the algorithm will be capturing more of the overall trend of the data, which overlooks the local variance and brings more bias, and the performance will be decreasing and reaching a plateau.
  - What we expect is to find out the sweet spot balancing the two kinds of reducible errors.

#### LASSO/Elastic Net

**Let's make a glmnet model and tune it's hyperparameter together!**

Brain imaging data could be very high-dimensional. And thus, regularization techniques are very helpful in building neuroimaging decoding models. (Decoding here refers to use the neuroimaging data to "predict" the information about the subjects and stimuli.) here I created a fake dataset simulating fMRI data. Each column will be a voxel (the spatial unit of the imaging, proxy toward a small region of cortex) and each row is the response time to a stimulus when the simulated subject was lying in the scanner watching visual stimuli and making decisions on it.

```{r}
# dependencies
options(conflicts.policy = "depends.ok")
devtools::source_url("https://github.com/jjcurtin/lab_support/blob/main/fun_ml.R?raw=true")
tidymodels_conflictRules()
conflictRules("Matrix", mask.ok = c("expand", "pack", "unpack"))

theme_set(theme_classic())

devtools::source_url("https://github.com/jjcurtin/lab_support/blob/main/fun_plots.R?raw=true")
devtools::source_url("https://github.com/jjcurtin/lab_support/blob/main/fun_eda.R?raw=true")
```

```{r}
# packages
library(tidymodels)
library(ggplot2)

# Set seed for reproducibility
set.seed(42)

# Define parameters
n_trials <- 500   # Number of stimulus-response trials
n_voxels <- 100  # Number of voxels (features)

# Simulate voxel activity (random normal distribution)
voxel_data <- matrix(rnorm(n_trials * n_voxels, mean = 0, sd = 1), 
                     nrow = n_trials, ncol = n_voxels)

# Simulate sparse true coefficients (only some voxels react relevant with response time)
true_coefficients <- rep(0, n_voxels)
true_coefficients[sample(1:n_voxels, 20)] <- rnorm(20, mean = 5, sd = 2)  # 20 predictive voxels

# Generate response time (y) with noise
response_times <- voxel_data %*% true_coefficients + rnorm(n_trials, mean = 500, sd = 50) # RT in ms

# Convert to data frame
brain_data <- as.data.frame(voxel_data)
colnames(brain_data) <- paste0("Voxel_", 1:n_voxels)  # Name columns
brain_data <- 
  brain_data |> 
  mutate(rt = as.numeric(response_times))  # Add response time column

# Print first few rows
head(brain_data[, c(1:5, 100, 101)]) |> view()
```

```{r}
skimr::skim(brain_data)
```

```{r}
set.seed(20210511)
splits <- brain_data |> 
  initial_split(prop = 3/4, strata = "rt", breaks = 4)

data_trn <- splits |> analysis()

data_test <- splits |> assessment()
```

We're going to predict the response time from the voxel activation data.

```{r}
rec <- recipe(rt ~ ., data = data_trn) |> 
  step_normalize(all_numeric_predictors())
  
rec_prep <- rec |> 
  prep(data_trn)

feat_trn <- rec_prep |> 
  bake(NULL)

feat_test <- rec_prep |> 
  bake(data_test)
```

```{r}
data_trn
```


```{r}
rerun_setting <- F
```

```{r}
set.seed(42)
splits_boot <- data_trn |> 
   bootstraps(times = 5, strata = "rt") 
```

```{r}
grid_glmnet <- expand_grid(penalty = exp(seq(-2, 4, length.out = 50)),
                           mixture = seq(0, 1, length.out = 6))
```

```{r}
fits_glmnet <- xfun::cache_rds(
  expr = {
    linear_reg(penalty = tune(), 
               mixture = tune()) |> 
    set_engine("glmnet") |> 
    tune_grid(preprocessor = rec, 
              resamples = splits_boot, 
              grid = grid_glmnet, 
              metrics = metric_set(rmse))
  
  },
  rerun = rerun_setting,
  dir = here::here("lab/cache//"),
  file = "brain_demo")
```

```{r}
plot_hyperparameters(fits_glmnet, hp1 = "penalty", hp2 = "mixture", metric = "rmse")
```

```{r}
select_best(fits_glmnet, metric = "rmse")
```

It's a pure LASSO model. Then we fit the best model with the configuration above.

```{r}
fit_glmnet <- linear_reg(penalty = select_best(fits_glmnet)$penalty,
                         mixture = select_best(fits_glmnet)$mixture) |> 
  set_engine("glmnet") |> 
  fit(rt ~ ., data = feat_trn) |> 
  suppressWarnings()
```

```{r}
fit_glmnet |> 
  tidy() |> 
  print(n = Inf)
```

```{r}
# number of non-zero feature coeffients in the end (-1, for the intercept)
fit_glmnet |> 
  tidy() |> 
  pull(estimate) |> 
  as.logical() |> 
  sum()
```

The selected features are less than the real predictive features. So, penalty (especially LASSO) could also cause oversparsity (in data with correlated features). We should be mindful with this trade-off if you indeed want correlated features at some point).


We could have a look on the performance:

```{r}
rmse_vec(truth = feat_trn$rt,
         estimate = predict(fit_glmnet, feat_trn)$.pred) |> 
  suppressWarnings()
```

```{r}
rmse_vec(truth = feat_test$rt,
         estimate = predict(fit_glmnet, feat_test)$.pred) |> 
  suppressWarnings()
```


### 4. Example of nested CV

We have a nice example in the slides. In inner loops we use bootstrapping CV and k-fold for the outer loop.

### 5. Other questions

- How does future_map() handle parallelization differently from foreach()? Can you provide examples where one might be preferred over the other, particularly in the context of resampling techniques like cross-validation?
  - Nothing further is needed to use `foreach()` or tidymodels functions.
  - For `future_map()`, you need to set up a multisession plan with this code chunk


- In Bayesian, I believe **LOO-CV** is frequently used. Do the same limits apply to **Bayesian, or is it just ML** that prefers K-fold over LOO-CV?
  - I will deem this as a decision based on sample size and compute capability. Because LOOCV is very computationally expsensive, we rarely use it in real world cases. If in your example of Bayesian analysis it’s not computationally expensive, LOOCV will tremendousely decrease the bias on estimation with acceptable variance. So, it will be helpful.

- If we use Grouped K-Fold, how do we verify that individuals (or units of analysis) are never split across training and testing folds?
  - The parameter `group` will group all the row-wise cases with same value in this column as a group, which works similarly as `group_by()` in logic. Maybe go through the doc of this function will help. (?group_vfold_cv in console.)

- I hope to have more explanation on hyper-parameter tuning (bootstrap tuning part).
  - We use multiple resampled data to carry out cross-validation, together with the hyperparameter tuning grid we could have estimate of the performance for each configuration in the grid and then decide which to use.
  - Bootstrpping is mainly different becuase we draw samples with repetition from the original data.
  
-  I would like to know some naming conventions for our models and fits considering we are going to have many of them in the future work and poor naming could cause potential chaos in our codes. 
  - Good convention of naming is important, for both personal tracking and team collaboration.
  - In the case of our course, use `fit_x` format should be fine. 
  - In general use, we may want and are able to use data structure to hold a set of models or other complex objects. In R, one useful type is **list**. It's formated as `list[[i]] <- obj`. And we could also dynamically add more objects (like models) to a list. This will merit more complicated oprations, like looping though the models.

