---
title: "Results for elastic net regression models"
author: "Matt Hawrilenko, PhD"
date: "`r Sys.Date()`"
output: 
  html_document:
    toc: yes
    toc_float: yes
    number_sections: yes
    code_folding: hide
    theme: lumen
  word_document:
      toc: yes
  pdf_document:
      toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      warning = FALSE,
                      message = FALSE,
                      fig.path = 'figures/')
options(DT.options = list(pageLength = 100, language = list(search = 'Filter:')))

knitr::opts_knit$set(root.dir = 'X:/Workspace/Schiz_Repo/Methods_Paper/')
```


# Setup
``` {r load-functions}
# Load packages
library(rio)
library(knitr)
library(tidyverse)
library(kableExtra)
library(table1)
library(caret)
library(ggplot2)
library(see)
library(ggforce)


# Load functions
source("MJH/02_models/01_CrossTrial_pipeline.R")
```


# Prep
``` {r trial-names}
study_labels <- data.frame(
  study_name = c(
    'SCH-3015',
    'SCH-305',
    'SCH-303',
    'SCH-302',
    'PSZ-3001',
    'All others'
  ),
  study_label = c(
    'Adults FE',
    'Chronic #1',
    'Chronic #2',
    'Older Adults',
    'Teens',
    'All others'
  )
)
```


# Table 1
``` {r table-1}
df_pascore <- readRDS("Methods_Paper/$Pre_submission_reviewed_files/baseline_predictors_survivors_dropped_imputed.tu.RDS") %>%
  select(DUSUBJID, PASCORETOTBASE) %>%
  mutate(DUSUBJID = str_replace(DUSUBJID, "R076477-", ""))

df_table1 <- df_analyze %>%
  select(-PASCORETOTBASE) %>% # Drop the standardized value
  left_join(., df_pascore) # pull in raw value

table1(~ resp_25 + resp_50 + rswg_remission + PASORETOTPERCENTCHANGE + PASCORETOTBASE | STUDYID, 
       overall = TRUE,
       data = df_table1)
```


# Elastic net modeling results

## RSWG remission

### Table of results
``` {r enet-rswg}
table_enet_cv_rswg <- import(file = "MJH/output/table_enet_cv_rswg.csv")
table_enet_xtrial_rswg <- import(file = "MJH/output/table_enet_xtrial_rswg.csv")
table_enet_loto_rswg <- import(file = "MJH/output/table_enet_loto_rswg.csv")

enet_rswg <- rbind(
  table_enet_cv_rswg %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_enet_xtrial_rswg %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_enet_loto_rswg %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
enet_rswg <- enet_rswg %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no")) 

# Make table
enet_rswg %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for elastic net models: RSWG outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")

# Make overall summary
enet_rswg %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-enet-rswg, fig.height=8, fig.width=14}

ggplot(enet_rswg, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
        geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25), 
          panel.border = element_rect(color="grey35", fill = NA),
          plot.caption = element_text(size=11)
        ) +
        labs(title ="Balanced accuracy across all study conditions",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance") 

ggsave(filename = "X:/Workspace/Schiz_Repo/Methods_Paper/MJH/02_models/Figure_1_vectorized.pdf", width=14, height=8, units = "in", device = cairo_pdf)
```



### Variable importances for all trials
``` {r plot-importance-rswg, fig.width = 12, fig.height = 10}
enet_xtrial_rswg <- readRDS(file = "MJH/ML objects/enet_xtrial_rswg.RData") # load data
#NOTE: To save memory space, we later saved the plotting df in ml objects/importances_enet_rswg.RData
# Preserving code as was originally run.

# Extract importances into df
df_imp <- list(varImp(enet_xtrial_rswg[[1]])$importance,
               varImp(enet_xtrial_rswg[[2]])$importance,
               varImp(enet_xtrial_rswg[[3]])$importance,
               varImp(enet_xtrial_rswg[[4]])$importance,
               varImp(enet_xtrial_rswg[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

  saveRDS(df_imp, file = "MJH/ML objects/importances_enet_rswg.RData")
```



## 25% improvement

### Table of results
``` {r enet-resp_25}
table_enet_cv_resp_25 <- import(file = "MJH/output/table_enet_cv_resp_25.csv")
table_enet_xtrial_resp_25 <- import(file = "MJH/output/table_enet_xtrial_resp_25.csv")
table_enet_loto_resp_25 <- import(file = "MJH/output/table_enet_loto_resp_25.csv")

enet_resp_25 <- rbind(
  table_enet_cv_resp_25 %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_enet_xtrial_resp_25 %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_enet_loto_resp_25 %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
enet_resp_25 <- enet_resp_25 %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no"))

# Make table
enet_resp_25 %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for elastic net models: 25% improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
enet_resp_25 %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-enet-resp_25, fig.height=8, fig.width=14}

ggplot(enet_resp_25, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
   geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25), 
          panel.border = element_rect(color="grey35", fill = NA),
          plot.caption = element_text(size=11)
        ) +
        labs(title ="Balanced accuracy across all study conditions: Elastic net modeling, 25% improvement",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance")
```



### Variable importances
``` {r plot-importance-resp_25, fig.width = 12, fig.height = 10}
enet_xtrial_resp_25 <- readRDS(file = "MJH/ML objects/enet_xtrial_resp_25.RData") # load data

# Extract importances into df
df_imp <- list(varImp(enet_xtrial_resp_25[[1]])$importance,
               varImp(enet_xtrial_resp_25[[2]])$importance,
               varImp(enet_xtrial_resp_25[[3]])$importance,
               varImp(enet_xtrial_resp_25[[4]])$importance,
               varImp(enet_xtrial_resp_25[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

 saveRDS(df_imp, file = "MJH/ML objects/importances_enet_resp25.RData")
```



## 50% improvement

### Table of results
``` {r enet-resp_50}
table_enet_cv_resp_50 <- import(file = "MJH/output/table_enet_cv_resp_50.csv")
table_enet_xtrial_resp_50 <- import(file = "MJH/output/table_enet_xtrial_resp_50.csv")
table_enet_loto_resp_50 <- import(file = "MJH/output/table_enet_loto_resp_50.csv")

enet_resp_50 <- rbind(
  table_enet_cv_resp_50 %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_enet_xtrial_resp_50 %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_enet_loto_resp_50 %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
enet_resp_50 <- enet_resp_50 %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no"))

# Make table
enet_resp_50 %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for elastic net models: 50% improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
enet_resp_50 %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-enet-resp_50, fig.height=8, fig.width=14}

ggplot(enet_resp_50, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
   geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25), 
          panel.border = element_rect(color="grey35", fill = NA),
          plot.caption = element_text(size=11)
        ) +
        labs(title ="Balanced accuracy across all study conditions: Elastic net modeling, 50% improvement",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance")
```



### Variable importances
``` {r plot-importance-resp_50, fig.width = 12, fig.height = 10}
rm(enet_xtrial_resp_50, enet_xtrial_resp_25) # drop these because server runs out of memory

enet_xtrial_resp_50 <- readRDS(file = "MJH/ML objects/enet_xtrial_resp_50.RData") # load data

# Extract importances into df
df_imp <- list(varImp(enet_xtrial_resp_50[[1]])$importance,
               varImp(enet_xtrial_resp_50[[2]])$importance,
               varImp(enet_xtrial_resp_50[[3]])$importance,
               varImp(enet_xtrial_resp_50[[4]])$importance,
               varImp(enet_xtrial_resp_50[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

 saveRDS(df_imp, file = "MJH/ML objects/importances_enet_resp50.RData")
 
```


## Percent total improvement

### Table of results
``` {r enet-percent_change}
table_enet_cv_pctchg <- import(file = "MJH/output/table_enet_cv_pctchg.csv")
table_enet_xtrial_pctchg <- import(file = "MJH/output/table_enet_xtrial_pctchg.csv")
table_enet_loto_pctchg <- import(file = "MJH/output/table_enet_loto_pctchg.csv")

enet_pctchg <- rbind(
  table_enet_cv_pctchg %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_enet_xtrial_pctchg %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_enet_loto_pctchg %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Make table
enet_pctchg %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, r2, lb, ub, RMSE) %>%
  mutate(across('r2':'ub', ~ round(.x, 3))) %>%
  mutate(RMSE = round(RMSE, 1)) %>%
  kbl(
    caption = "Classification quality for elastic net models: Percent total improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "RMSE")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  add_header_above(c(" " = 2, "R-squared" = 3, " " = 1))  %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
enet_pctchg %>%
  group_by(condition) %>%
  summarize(across(c(r2, RMSE), ~ mean(.x))) %>%
  mutate(across('r2':'RMSE', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "R-squared", "RMSE")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: R-squared
``` {r plot-enet-pctchg, fig.height=8, fig.width=14}

ggplot(enet_pctchg, aes(x = study_pair, y = r2)) +
   geom_errorbar(aes(ymin=lb, ymax=ub), color = "light grey", width=0.3) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25), 
          panel.border = element_rect(color="grey35", fill = NA),
          plot.caption = element_text(size=11)
        ) +
        labs(title ="R-squared across all study conditions: Elastic net modeling, Percent total improvement",
             x = "Study pair",
             y = "R-squared",
             caption = "FE = first episode")
```



### Variable importances
``` {r plot-importance-percent_change, fig.width = 12, fig.height = 10}
enet_xtrial_pctchg <- readRDS(file = "MJH/ML objects/enet_xtrial_pctchg.RData") # load data

# Extract importances into df
df_imp <- list(varImp(enet_xtrial_pctchg[[1]])$importance,
               varImp(enet_xtrial_pctchg[[2]])$importance,
               varImp(enet_xtrial_pctchg[[3]])$importance,
               varImp(enet_xtrial_pctchg[[4]])$importance,
               varImp(enet_xtrial_pctchg[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))
 
 saveRDS(df_imp, file = "MJH/ML objects/importances_enet_r2.RData")
```



# random forest modeling results

## RSWG remission

### Table of results
``` {r rf-rswg}
table_rf_cv_rswg <- import(file = "MJH/output/table_rf_cv_rswg.csv")
table_rf_xtrial_rswg <- import(file = "MJH/output/table_rf_xtrial_rswg.csv")
table_rf_loto_rswg <- import(file = "MJH/output/table_rf_loto_rswg.csv")

rf_rswg <- rbind(
  table_rf_cv_rswg %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_rf_xtrial_rswg %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_rf_loto_rswg %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
rf_rswg <- rf_rswg %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no"))

# Make table
rf_rswg %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for random forest models: RSWG outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
rf_rswg %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-rf-rswg, fig.height=8, fig.width=14}

ggplot(rf_rswg, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
   geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25)
        ) +
        labs(title ="Balanced accuracy across all study conditions",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance")
```



### Variable importances for all trials
``` {r plot-importance-rswg, fig.width = 12, fig.height = 10}
rf_xtrial_rswg <- readRDS(file = "MJH/ML objects/rf_xtrial_rswg.RData") # load data

# Extract importances into df
df_imp <- list(varImp(rf_xtrial_rswg[[1]])$importance,
               varImp(rf_xtrial_rswg[[2]])$importance,
               varImp(rf_xtrial_rswg[[3]])$importance,
               varImp(rf_xtrial_rswg[[4]])$importance,
               varImp(rf_xtrial_rswg[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

 saveRDS(df_imp, file = "MJH/ML objects/importances_rf_rswg.RData")
```



## 25% improvement

### Table of results
``` {r rf-resp_25}
table_rf_cv_resp_25 <- import(file = "MJH/output/table_rf_cv_resp_25.csv")
table_rf_xtrial_resp_25 <- import(file = "MJH/output/table_rf_xtrial_resp_25.csv")
table_rf_loto_resp_25 <- import(file = "MJH/output/table_rf_loto_resp_25.csv")

rf_resp_25 <- rbind(
  table_rf_cv_resp_25 %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_rf_xtrial_resp_25 %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_rf_loto_resp_25 %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
rf_resp_25 <- rf_resp_25 %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no"))

# Make table
rf_resp_25 %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for random forest models: 25% improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
rf_resp_25 %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-rf-resp_25, fig.height=8, fig.width=14}

ggplot(rf_resp_25, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
   geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25)
        ) +
        labs(title ="Balanced accuracy across all study conditions: random forest modeling, 25% improvement",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance")
```



### Variable importances
``` {r plot-importance-resp_25, fig.width = 12, fig.height = 10}
rf_xtrial_resp_25 <- readRDS(file = "MJH/ML objects/rf_xtrial_resp_25.RData") # load data

# Extract importances into df
df_imp <- list(varImp(rf_xtrial_resp_25[[1]])$importance,
               varImp(rf_xtrial_resp_25[[2]])$importance,
               varImp(rf_xtrial_resp_25[[3]])$importance,
               varImp(rf_xtrial_resp_25[[4]])$importance,
               varImp(rf_xtrial_resp_25[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))
 
 saveRDS(df_imp, file = "MJH/ML objects/importances_rf_resp25.RData")
```



## 50% improvement

### Table of results
``` {r rf-resp_50}
table_rf_cv_resp_50 <- import(file = "MJH/output/table_rf_cv_resp_50.csv")
table_rf_xtrial_resp_50 <- import(file = "MJH/output/table_rf_xtrial_resp_50.csv")
table_rf_loto_resp_50 <- import(file = "MJH/output/table_rf_loto_resp_50.csv")

rf_resp_50 <- rbind(
  table_rf_cv_resp_50 %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_rf_xtrial_resp_50 %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_rf_loto_resp_50 %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Adjust p-values
rf_resp_50 <- rf_resp_50 %>%
  mutate( p_adj = p.adjust(ba_p)) %>%
  mutate(sig_vs_noinfo = ifelse(p_adj < .05, "yes", "no"))

# Make table
rf_resp_50 %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, balanced_accuracy, ba_lb, ba_ub, p_adj, sensitivity, specificity, accuracy, no_info_rate, auc) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  mutate(p_adj = case_when(
    p_adj < .001 ~ "< .001",
    p_adj > .999 ~ "> .999",
    TRUE ~ as.character(p_adj))) %>%
  mutate(p_adj = str_replace(p_adj, "0[.]", ".")) %>%
  kbl(
    caption = "Classification quality for random forest models: 50% improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "Adjusted p-value", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  column_spec(6, width = "5em") %>%
  add_header_above(c(" " = 2, "Balanced accuracy" = 4, "Other metrics" = 5)) %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")

# Make overall summary
rf_resp_50 %>%
  group_by(condition) %>%
  summarize(across(c(balanced_accuracy, sensitivity, specificity, accuracy, no_info_rate, auc), ~ mean(.x))) %>%
  mutate(across('balanced_accuracy':'auc', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "Balanced accuracy", "Sensitivity", "Specificity", "Accuracy", "No info rate", "AUC")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: Balanced accuracy
``` {r plot-rf-resp_50, fig.height=8, fig.width=14}

ggplot(rf_resp_50, aes(x = study_pair, y = balanced_accuracy, color = sig_vs_noinfo)) +
   geom_errorbar(aes(ymin=ba_lb, ymax=ba_ub), color = "light grey", width=0.3) +
        geom_hline(yintercept=0.5, linetype="dashed", color="black", size=.7) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        scale_y_continuous(limits=c(.43, 1), 
                           breaks = seq(0.5, 1, by=0.1)) +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        scale_color_manual(values = c("dark grey", "red")) +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25)
        ) +
        labs(title ="Balanced accuracy across all study conditions: random forest modeling, 50% improvement",
             x = "Study pair",
             y = "Balanced accuracy",
             caption = "FE = first episode", 
             color = "Balanced accuracy is higher than chance")
```



### Variable importances
``` {r plot-importance-resp_50, fig.width = 12, fig.height = 10}

rf_xtrial_resp_50 <- readRDS(file = "MJH/ML objects/rf_xtrial_resp_50.RData") # load data

# Extract importances into df
df_imp <- list(varImp(rf_xtrial_resp_50[[1]])$importance,
               varImp(rf_xtrial_resp_50[[2]])$importance,
               varImp(rf_xtrial_resp_50[[3]])$importance,
               varImp(rf_xtrial_resp_50[[4]])$importance,
               varImp(rf_xtrial_resp_50[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

 saveRDS(df_imp, file = "MJH/ML objects/importances_rf_resp50.RData")
 
```


## Percent total improvement

### Table of results
``` {r rf-percent_change}
table_rf_cv_pctchg <- import(file = "MJH/output/table_rf_cv_pctchg.csv")
table_rf_xtrial_pctchg <- import(file = "MJH/output/table_rf_xtrial_pctchg.csv")
table_rf_loto_pctchg <- import(file = "MJH/output/table_rf_loto_pctchg.csv")

rf_pctchg <- rbind(
  table_rf_cv_pctchg %>%
    mutate(condition = "Within-trial CV")
  ,
  
  table_rf_xtrial_pctchg %>%
    mutate(  condition = case_when(
    predictor == target ~ "Within-trial no validation",
    predictor != target ~ "Paired-trial"
  ))
  , 
  
  table_rf_loto_pctchg %>%
    mutate(condition = "Leave-one-trial-out")
  ) %>%
  
  # Add shorthand labels
  left_join(study_labels, by = c("predictor" = "study_name")) %>% # predictors
  rename(predictor_label = study_label) %>%
  
  left_join(study_labels, by = c("target" = "study_name")) %>% # targets
  rename(target_label = study_label) %>%
  
  mutate(
    study_pair = paste0(predictor_label, " \U2192 ", target_label),
    condition = fct_relevel(condition, "Within-trial no validation",  "Within-trial CV", "Paired-trial", "Leave-one-trial-out")
  )

# Make table
rf_pctchg %>%
  arrange(condition, predictor_label, target_label) %>%
  select(predictor_label, target_label, r2, lb, ub, RMSE) %>%
  mutate(across('r2':'ub', ~ round(.x, 3))) %>%
  mutate(RMSE = round(RMSE, 1)) %>%
  kbl(
    caption = "Classification quality for random forest models: Percent total improvement outcome",
    col.names = c("Predictor", "Target", "Estimate", "Lower", "Upper", "RMSE")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  pack_rows("Within-trial no validation", 1, 5) %>%
  pack_rows("Within-trial cross-validation", 6, 10) %>%
  pack_rows("Paired-trial", 11, 30) %>%
  pack_rows("Leave-one-trial-out", 31, 35) %>%
  add_header_above(c(" " = 2, "R-squared" = 3, " " = 1))  %>%
  footnote(general = "Lower/Upper = 95% confidence interval lower and upper bounds.")


# Make overall summary
rf_pctchg %>%
  group_by(condition) %>%
  summarize(across(c(r2, RMSE), ~ mean(.x))) %>%
  mutate(across('r2':'RMSE', ~ round(.x, 3))) %>%
  kbl(
    caption = "Summary of classification quality metrics across all prediction scenarios",
    col.names = c("Condition", "R-squared", "RMSE")
    ) %>%
  kable_classic(full_width = F, html_font = "Cambria")

```

### Figure: R-squared
``` {r plot-rf-pctchg, fig.height=8, fig.width=14}

ggplot(rf_pctchg, aes(x = study_pair, y = r2)) +
   geom_errorbar(aes(ymin=lb, ymax=ub), color = "light grey", width=0.3) +
        geom_point(size=3, shape="triangle") +
        theme_modern() +
        facet_grid(.~ condition, scales = "free_x", space = "free_x") +
        theme(
          axis.text.x = element_text(angle=60, vjust=1, hjust=1),
          legend.position = "bottom",
          axis.ticks.x.bottom = element_line(),
          strip.text.x = element_text(size = 11),
          panel.grid.major.y = element_line(colour="light grey", size=0.25)
        ) +
        labs(title ="R-squared across all study conditions: random forest modeling, Percent total improvement",
             x = "Study pair",
             y = "R-squared",
             caption = "FE = first episode")
```



### Variable importances
``` {r plot-importance-percent_change, fig.width = 12, fig.height = 10}
rf_xtrial_pctchg <- readRDS(file = "MJH/ML objects/rf_xtrial_pctchg.RData") # load data

# Extract importances into df
df_imp <- list(varImp(rf_xtrial_pctchg[[1]])$importance,
               varImp(rf_xtrial_pctchg[[2]])$importance,
               varImp(rf_xtrial_pctchg[[3]])$importance,
               varImp(rf_xtrial_pctchg[[4]])$importance,
               varImp(rf_xtrial_pctchg[[5]])$importance) %>%
  map_df(rownames_to_column, "predictor", .id = "study") %>%
  rename("importance" = "Overall") %>%
  group_by(predictor) %>%
  mutate(
    mean_importance = mean(importance),
    predictor = str_to_title(predictor)
  ) %>%
  ungroup() %>%
  mutate(rank = dense_rank(-mean_importance)) %>%
  arrange(desc(mean_importance), predictor) %>%
  mutate(study_name = case_when(
    study == 1 ~ study_names[1],
    study == 2 ~ study_names[2],
    study == 3 ~ study_names[3],
    study == 4 ~ study_names[4],
    study == 5 ~ study_names[5]
  )) %>%
  left_join(study_labels) %>%
  filter(rank <= 50)

# Plot
ggplot(data = df_imp, aes(x = reorder(predictor, mean_importance), y = importance, fill = study_label)) +
  geom_bar(stat = "identity") +
  coord_flip() + 
  facet_wrap(~ study_label, ncol = 5) +
  theme_minimal(base_family = "serif", base_size = 12) +
  ggtitle("Variable importance across all studies") +
  xlab("") +
  ylab("Variable importance") +
  theme(
    legend.position = "none",
    axis.text = element_text(size=12))

 saveRDS(df_imp, file = "MJH/ML objects/importances_rf_r2.RData")
 
```
