README.Rmd

---
title: "VEESA R Package"
output: rmarkdown::github_document
always_allow_html: yes
---

<!-- badges: start -->
[![R-CMD-check](https://github.com/sandialabs/veesa/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/sandialabs/veesa/actions/workflows/R-CMD-check.yaml)
[![Codecov test coverage](https://codecov.io/gh/sandialabs/veesa/graph/badge.svg)](https://app.codecov.io/gh/sandialabs/veesa)
<!-- badges: end -->
  
```{r setup, include = FALSE}
# Rmarkdown options
knitr::opts_chunk$set(echo = TRUE, message = FALSE, dpi = 300)
```

`veesa` is an R package for implementing the VEESA pipeline for an explainable approach to training machine learning models with functional data inputs. See a preprint manuscript describing the approach on [arXiv](https://arxiv.org/abs/2501.07602). Installing `veesa` can be implemented using either of the commands below.

```{r}
#| eval: false

# CRAN
install.packages("veesa")

# Development version from GitHub
remotes::install_github("sandialabs/veesa")
```

Keep reading for an example using `veesa` to implement the VEESA pipeline.

## Demonstration

#### Set Up and Data Generation

```{r}
#| message: false
# Load R packages
library(cowplot)
library(dplyr)
library(ggplot2)
library(purrr)
library(randomForest)
library(tidyr)
library(veesa)

# Specify a color palette
color_pal = wesanderson::wes_palette("Zissou1", 5, type = "continuous")

# Specify colors for PC direction plots
col_plus1 = "#784D8C"
col_plus2 = "#A289AE"
col_minus1 = "#EA9B44"
col_minus2 = "#EBBC88"
col_pcdir_1sd = c(col_plus1, "black", col_minus1)
col_pcdir_2sd = c(col_plus2, col_plus1, "black", col_minus1, col_minus2)
```

Simulate data:

```{r}
sim_data = simulate_functions(M = 100, N = 75, seed = 20211130)
```

Separate data into training/testing:

```{r}
set.seed(20211130)
id = unique(sim_data$id)
M_test = length(id) * 0.25
id_test = sample(x = id, size = M_test, replace = FALSE)
sim_data = sim_data %>% mutate(data = ifelse(id %in% id_test, "test", "train"))
```

Simulated functions colored by covariates: 

```{r echo = FALSE}
#| fig-height: 5
#| fig-width: 27
plot_sim <- function(cov) {
  sim_data %>%
    ggplot(aes(x = t, y = y, color = get(cov), group = id)) + 
    geom_line(alpha = 0.75) + 
    scale_color_gradientn(colours =  color_pal) +
    theme_bw(base_size = 20) + 
    labs(color = cov)
}
sim_plot = map(.x = c("x1", "x2", "x3"), .f = plot_sim)
plot_grid(plotlist = sim_plot, ncol = 3, byrow = FALSE)
```

Prepare matrices from the data frames:

```{r}
prep_matrix <- function(df, train_test) {
  df %>%
    filter(data == train_test) %>%
    select(id, t, y) %>%
    ungroup() %>%
    pivot_wider(id_cols = t,
                names_from = id,
                values_from = y) %>%
    select(-t) %>%
    as.matrix()
}
sim_train_matrix = prep_matrix(df = sim_data, train_test = "train")
sim_test_matrix = prep_matrix(df = sim_data, train_test = "test")
```

Create a vector of times:

```{r}
times = sim_data$t %>% unique()
```

#### Alignment and fPCA

Prepare train data

```{r}
train_transformed_jfpca <-
  prep_training_data(
    f = sim_train_matrix,
    time = times, 
    fpca_method = "jfpca",
    optim_method = "DPo"
  )
```

Prepare test data:

```{r}
test_transformed_jfpca <-
  prep_testing_data(
    f = sim_test_matrix,
    time = times,
    train_prep = train_transformed_jfpca,
    optim_method = "DPo"
  )
```

Plot several PCs:

```{r}
#| fig-height: 5
#| fig-width: 20
#| echo: false
plot_pc_directions(
  fpcs = 1:3,
  fdasrvf = train_transformed_jfpca$fpca_res,
  fpca_method = "jfpca",
  time = times,
  linesizes = rep(1, 5)
) + 
  scale_color_manual(values = col_pcdir_2sd) + 
  theme_bw(base_size = 20) + 
  labs(
    x = "Time",
    y = "Intensity"
  )
```

Compare jfPCA coefficients from train and test data:

```{r}
#| fig-height: 4
#| fig-width: 10
#| out-width: "75%"
#| fig-align: "center"
#| echo: false

train_plot_df_jfpca <- 
  train_transformed_jfpca$fpca_res$coef %>%
  t() %>%
  data.frame() %>%
  mutate(pc = 1:n()) %>%
  pivot_longer(cols = -pc, names_to = "id") %>%
  mutate(data = "train")

test_plot_df_jfpca <- 
  test_transformed_jfpca$coef %>%
  data.frame() %>%
  mutate(pc = 1:n()) %>%
  pivot_longer(cols = -pc, names_to = "id") %>%
  mutate(data = "test")

bind_rows(train_plot_df_jfpca, test_plot_df_jfpca) %>%
  ggplot(aes(
    x = pc,
    y = value,
    group = factor(id):factor(data),
    color = data
  )) +
  geom_line(alpha = 0.5) +
  labs(
    title = "Joint fPCA", 
    color = "Data",
    x = "PC number",
    y = "Coefficient") +
  theme_bw()
```

#### Models

Create response variable: 

```{r}
x1_train <- 
  sim_data %>% filter(data == "train") %>%
  select(id, x1) %>%
  distinct() %>% 
  pull(x1)
```

Create data frame with PCs and response for random forest:

```{r}
rf_jfpca_df <- 
  train_transformed_jfpca$fpca_res$coef %>%
  data.frame() %>%
  rename_all(.funs = function(x) stringr::str_replace(x, "X", "pc")) %>%
  mutate(x1 = x1_train) %>%
  select(x1, everything())
```

Fit random forest:

```{r}
set.seed(20211130)
rf_jfpca = randomForest(x1 ~ ., data = rf_jfpca_df)
```

#### PFI

Compute PFI:

```{r}
set.seed(20211130)
pfi_jfpca <- compute_pfi(
  x = rf_jfpca_df %>% select(-x1),
  y = rf_jfpca_df$x1,
  f = rf_jfpca,
  K = 10,
  metric = "nmse"
)
```

PFI results (mean of reps):

```{r}
#| fig-height: 4.5
#| fig-width: 12
#| out-width: "75%"
#| fig-align: "center"
#| echo: false
data.frame(pfi = pfi_jfpca$pfi) %>%
  mutate(pc = 1:n()) %>%
  ggplot(aes(x = pc, y = pfi)) +
  geom_point() +
  geom_segment(aes(yend = 0, xend = pc)) +
  theme_bw(base_size = 14) +
  labs(
    x = "Principal Component",
    y = "Permutation Feature Importance"
  )
```

PFI results (variability across reps):

```{r}
#| fig-height: 4.5
#| fig-width: 12 
#| out-width: "75%"
#| fig-align: "center"
#| echo: false
pfi_jfpca$pfi_single_reps %>%
  data.frame() %>%
  mutate(rep = 1:n()) %>%
  pivot_longer(cols = -rep,
               names_to = "pc",
               values_to = "pfi") %>%
  mutate(pc = stringr::str_remove(pc, "X")) %>%
  mutate(pc = as.numeric(pc)) %>%
  ggplot(aes(x = pc, y = pfi, group = pc)) +
  geom_boxplot() +
  theme_bw(base_size = 14) +
  labs(
    x = "Principal Component",
    y = "Permutation Feature Importance"
  )
```

Identify the top PC for each elastic fPCA method:

```{r}
top_pc_jfpca <- 
  data.frame(pfi = pfi_jfpca$pfi) %>%
  mutate(pc = 1:n()) %>%
  arrange(desc(pfi)) %>%
  slice(1) %>%
  pull(pc)
```

Principal directions of top PC for each jfPCA method:

```{r}
#| fig-height: 4
#| fig-width: 8 
#| out-width: "60%"
#| fig-align: "center"
#| echo: false
plot_pc_directions(
  fpcs = top_pc_jfpca,
  fdasrvf = train_transformed_jfpca$fpca_res,
  fpca_method = "jfpca",
  nrow = 2, 
  linesizes = rep(1, 5)
) +
  scale_color_manual(values = col_pcdir_2sd) +
  labs(title = "Top PC for Joint fPCA")
```