library(dplyr)
library(tidyr)
library(cluster)
library(ggrepel)
library(FactoMineR)
library(ggplot2)  
library(factoextra) 
library(ggforce) 
library(psych)
library(plotly)
library(purrr) 
library(stringr)
library(lmerTest)
library(reshape2)
library(tidyverse)
library(pls)
library(corrplot)
library(glmnet)
library(vip)
library(readxl)


# ============================================================================
# FUNCTION 1: Process individual metric data with flexible subject naming
# ============================================================================

process_metric_data <- function(
    data_file_33,      # CSV file for timepoint 33
    data_file_40,      # CSV file for timepoint 40
    metric_name,       # e.g., "BOLD", "MK", "RK"
    output_dir         # Where to save results
) {
  
  cat("\n========================================\n")
  cat("Processing:", metric_name, "\n")
  cat("========================================\n")
  
  # Create output directory if it doesn't exist
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  # Read data
  data_33 <- read.csv(data_file_33)
  data_40 <- read.csv(data_file_40)
  
  # Rename first column to "subject"
  names(data_33)[1] <- "subject"
  names(data_40)[1] <- "subject"
  
  # DIAGNOSTIC: Print original subject IDs
  cat("\nOriginal subject IDs in 33 weeks file (first 5):\n")
  print(head(data_33$subject, 5))
  cat("\nOriginal subject IDs in 40 weeks file (first 5):\n")
  print(head(data_40$subject, 5))
  
  # Clean subject labels - extract numbers after underscore OR after "Subject_"
  # This works for: "T0_04", "T1_04", "Subject_002", or just "004"
  data_33 <- data_33 %>%
    mutate(subject_original = subject,  # Keep original for debugging
           subject = case_when(
             str_detect(subject, "_") ~ str_extract(subject, "(?<=_)[0-9]+"),  # After underscore
             str_detect(subject, "^[0-9]+$") ~ subject,  # Already just numbers
             TRUE ~ str_extract(subject, "[0-9]+")  # Any numbers
           ))
  
  data_40 <- data_40 %>%
    mutate(subject_original = subject,
           subject = case_when(
             str_detect(subject, "_") ~ str_extract(subject, "(?<=_)[0-9]+"),  # After underscore
             str_detect(subject, "^[0-9]+$") ~ subject,  # Already just numbers
             TRUE ~ str_extract(subject, "[0-9]+")  # Any numbers
           ))
  
  # Remove empty rows
  data_33 <- data_33 %>% filter(!is.na(subject) & subject != "")
  data_40 <- data_40 %>% filter(!is.na(subject) & subject != "")
  
  # DIAGNOSTIC: Print cleaned subject IDs
  cat("\nCleaned subject IDs in 33 weeks (first 5):\n")
  print(head(data_33$subject, 5))
  cat("\nCleaned subject IDs in 40 weeks (first 5):\n")
  print(head(data_40$subject, 5))
  
  # DIAGNOSTIC: Check for matches
  cat("\nNumber of subjects in 33 weeks:", nrow(data_33), "\n")
  cat("Number of subjects in 40 weeks:", nrow(data_40), "\n")
  cat("Number of matching subjects:", sum(data_33$subject %in% data_40$subject), "\n")
  
  if (sum(data_33$subject %in% data_40$subject) == 0) {
    cat("\n⚠️ WARNING: NO MATCHING SUBJECTS FOUND!\n")
    cat("This likely means the subject ID format is different between files.\n")
    cat("\nShowing unmatched subjects from 33 weeks:\n")
    print(setdiff(data_33$subject, data_40$subject))
    cat("\nShowing unmatched subjects from 40 weeks:\n")
    print(setdiff(data_40$subject, data_33$subject))
    stop("Cannot proceed - no matching subjects between timepoints")
  }
  
  # Remove the original column before merging
  data_33 <- data_33 %>% select(-subject_original)
  data_40 <- data_40 %>% select(-subject_original)
  
  # Define regions
  regions <- c("limbic", "paralimbic", "thalamus", "PFC", "visual",
               "auditory", "PCC", "SSM", "PCUN")
  
  # Rename columns
  names(data_33)[2:10] <- paste0(regions, "_33")
  names(data_40)[2:10] <- paste0(regions, "_40")
  
  # Merge datasets
  data_merged <- data_33 %>%
    inner_join(data_40, by = "subject")
  
  # Compute deltas
  for (r in regions) {
    data_merged[[paste0(r, "_delta")]] <- 
      data_merged[[paste0(r, "_40")]] - data_merged[[paste0(r, "_33")]]
  }
  
  # Save merged data
  write.csv(data_merged, 
            file.path(output_dir, paste0(metric_name, "_merged.csv")), 
            row.names = FALSE)
  
  data_merged <- data_merged %>%
    mutate(subject = as.numeric(as.character(subject))) %>%
    filter(!subject %in% c(05, 06, 07, 49))
  
  # Save merged data
  write.csv(data_merged, 
            file.path(output_dir, paste0(metric_name, "_merged_NoMotion.csv")), 
            row.names = FALSE)
  
  # ==================== STATISTICAL ANALYSIS (WITH OUTLIERS) ====================
  
  results <- data.frame(
    Region = regions,
    mean_delta = NA,
    sd_delta = NA,
    CI_low = NA,
    CI_high = NA,
    t_value = NA,
    p_value = NA
  )
  
  for (r in regions) {
    t_res <- t.test(
      data_merged[[paste0(r, "_40")]], 
      data_merged[[paste0(r, "_33")]], 
      paired = TRUE
    )
    
    delta_vals <- data_merged[[paste0(r, "_40")]] - data_merged[[paste0(r, "_33")]]
    
    results[results$Region == r, "mean_delta"] <- t_res$estimate
    results[results$Region == r, "sd_delta"] <- sd(delta_vals, na.rm = TRUE)
    results[results$Region == r, "CI_low"] <- t_res$conf.int[1]
    results[results$Region == r, "CI_high"] <- t_res$conf.int[2]
    results[results$Region == r, "t_value"] <- t_res$statistic
    results[results$Region == r, "p_value"] <- t_res$p.value
  }
  
  # FDR correction
  results$p_FDR <- p.adjust(results$p_value, method = "fdr")
  
  # Significance stars
  results$significance <- cut(
    results$p_FDR,
    breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
    labels = c("***", "**", "*", "ns")
  )
  
  write.csv(results, 
            file.path(output_dir, paste0(metric_name, "_results.csv")), 
            row.names = FALSE)
  
  # ==================== OUTLIER DETECTION ====================
  
  # Convert to long format
  data_long <- data_merged %>%
    pivot_longer(
      cols = -subject,
      names_to = c("region", "timepoint"),
      names_sep = "_",
      values_to = "value"
    ) %>%
    mutate(
      timepoint = factor(timepoint, levels = c("33", "40", "delta")),
      region = factor(region),
      ID = subject
    )
  
  write.csv(data_long, 
            file.path(output_dir, paste0(metric_name, "_long.csv")), 
            row.names = FALSE)
  
  # Flag outliers
  data_flagged <- data_long %>%
    group_by(region, timepoint) %>%
    mutate(
      Q1 = quantile(value, 0.25, na.rm = TRUE),
      Q3 = quantile(value, 0.75, na.rm = TRUE),
      IQR_val = IQR(value, na.rm = TRUE),
      lower = Q1 - 1.5 * IQR_val,
      upper = Q3 + 1.5 * IQR_val,
      outlier = value < lower | value > upper
    ) %>%
    ungroup()
  
  # Save outlier tables
  outlier_table <- data_flagged %>% filter(outlier == TRUE)
  write.csv(outlier_table, 
            file.path(output_dir, paste0(metric_name, "_all_outliers.csv")), 
            row.names = FALSE)
  
  # Outliers by timepoint and region
  outliers_33 <- data_flagged %>% filter(timepoint == "33", outlier == TRUE)
  outliers_40 <- data_flagged %>% filter(timepoint == "40", outlier == TRUE)
  outliers_delta <- data_flagged %>% filter(timepoint == "delta", outlier == TRUE)
  
  # By region
  outliers_33_region <- outliers_33 %>%
    group_by(region) %>%
    summarise(n_outliers = n(), IDs = paste(unique(ID), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  outliers_40_region <- outliers_40 %>%
    group_by(region) %>%
    summarise(n_outliers = n(), IDs = paste(unique(ID), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  outliers_delta_region <- outliers_delta %>%
    group_by(region) %>%
    summarise(n_outliers = n(), IDs = paste(unique(ID), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  write.csv(outliers_33_region, 
            file.path(output_dir, paste0(metric_name, "_outliers_33_region.csv")), 
            row.names = FALSE)
  write.csv(outliers_40_region, 
            file.path(output_dir, paste0(metric_name, "_outliers_40_region.csv")), 
            row.names = FALSE)
  write.csv(outliers_delta_region, 
            file.path(output_dir, paste0(metric_name, "_outliers_delta_region.csv")), 
            row.names = FALSE)
  
  # By ID
  outliers_33_ID <- outliers_33 %>%
    group_by(ID) %>%
    summarise(n_outliers = n(), regions = paste(unique(region), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  outliers_40_ID <- outliers_40 %>%
    group_by(ID) %>%
    summarise(n_outliers = n(), regions = paste(unique(region), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  outliers_delta_ID <- outliers_delta %>%
    group_by(ID) %>%
    summarise(n_outliers = n(), regions = paste(unique(region), collapse = ", ")) %>%
    arrange(desc(n_outliers))
  
  write.csv(outliers_33_ID, 
            file.path(output_dir, paste0(metric_name, "_outliers_33_ID.csv")), 
            row.names = FALSE)
  write.csv(outliers_40_ID, 
            file.path(output_dir, paste0(metric_name, "_outliers_40_ID.csv")), 
            row.names = FALSE)
  write.csv(outliers_delta_ID, 
            file.path(output_dir, paste0(metric_name, "_outliers_delta_ID.csv")), 
            row.names = FALSE)
  
  # ==================== VISUALIZATIONS ====================
  
  # Boxplot for 33 and 40 weeks
  data_plot <- data_flagged %>% filter(timepoint %in% c("33","40"))
  
  data_plot <- data_plot %>%
    mutate(
      region_num = as.numeric(region),
      region_dodge = region_num + ifelse(timepoint == "33", -0.2, 0.2)
    )
  
  p1 <- ggplot(data_plot, aes(x = region_num, y = value, fill = timepoint)) +
    geom_boxplot(
      aes(group = interaction(region, timepoint)),
      outlier.shape = NA,
      width = 0.4
    ) +
    geom_point(
      data = data_plot %>% filter(outlier == TRUE),
      aes(x = region_dodge, y = value),
      color = "red",
      size = 2
    ) +
    geom_text_repel(
      data = data_plot %>% filter(outlier == TRUE),
      aes(x = region_dodge, y = value, label = ID),
      size = 3,
      nudge_y = 0.1,
      segment.color = 'gray50'
    ) +
    scale_x_continuous(
      breaks = 1:length(levels(data_plot$region)),
      labels = levels(data_plot$region)
    ) +
    scale_fill_manual(values = c("33" = "skyblue", "40" = "orange")) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
      title = paste0(metric_name, " - Outliers at 33 and 40 weeks"),
      x = "Region",
      y = "Value",
      fill = "Timepoint"
    )
  
  ggsave(
    file.path(output_dir, paste0(metric_name, "_boxplot_outliers_33_40.png")), 
    p1, 
    width = 12, 
    height = 6, 
    dpi = 300
  )
  
  # Boxplot for delta with outliers
  data_delta <- data_flagged %>% filter(timepoint == "delta")
  
  region_order <- data_delta %>%
    group_by(region) %>%
    summarise(mean_delta = mean(value, na.rm = TRUE)) %>%
    arrange(desc(mean_delta)) %>%
    pull(region)
  
  data_delta <- data_delta %>%
    mutate(region = factor(region, levels = region_order))
  
  p2 <- ggplot(data_delta, aes(x = region, y = value)) +
    geom_boxplot(
      outlier.shape = NA,
      fill = "orange",
      width = 0.6
    ) +
    geom_point(
      data = data_delta %>% filter(outlier == TRUE),
      aes(x = region, y = value),
      color = "red",
      size = 2
    ) +
    geom_text_repel(
      data = data_delta %>% filter(outlier == TRUE),
      aes(x = region, y = value, label = ID),
      size = 3,
      nudge_y = 0.1,
      segment.color = "gray50"
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
      title = paste0(metric_name, " - Delta (40-33) with Outliers Flagged"),
      x = "Region",
      y = "Delta Value"
    )
  
  ggsave(
    file.path(output_dir, paste0(metric_name, "_boxplot_outliers_delta.png")), 
    p2, 
    width = 12, 
    height = 6, 
    dpi = 300
  )
  
  cat("\n✓ Preprocessing complete for", metric_name, "\n")
  cat("✓ Check the output folder for outlier tables\n")
  cat("✓ Review outliers before running outlier removal\n\n")
  
  return(list(
    results = results,
    data_merged = data_merged,
    data_long = data_long,
    outliers = outlier_table
  ))
}

# ============================================================================
# FUNCTION 2: Remove outliers and reanalyze
# ============================================================================

remove_outliers_and_reanalyze <- function(
    metric_name,
    output_dir,
    outlier_list_delta = NULL   # Named list: region = c(IDs)
) {
  
  cat("\n========================================\n")
  cat("Removing outliers and reanalyzing:", metric_name, "\n")
  cat("========================================\n")
  
  # Read merged data
  data_merged <- read.csv(file.path(output_dir, paste0(metric_name, "_merged_NoMotion.csv")))
  
  regions <- c("limbic", "paralimbic", "thalamus", "PFC", "visual",
               "auditory", "PCC", "SSM", "PCUN")
  
  # Create copy for outlier removal
  data_noOutliers <- data_merged
  
  # Remove delta outliers
  if (!is.null(outlier_list_delta)) {
    cat("Removing delta outliers...\n")
    for (r in names(outlier_list_delta)) {
      outlier_IDs <- outlier_list_delta[[r]]
      cols_to_remove <- c(paste0(r, "_33"), paste0(r, "_40"), paste0(r, "_delta"))
      data_noOutliers[data_noOutliers$subject %in% outlier_IDs, cols_to_remove] <- NA
    }
  }
  
  write.csv(data_noOutliers, 
            file.path(output_dir, paste0(metric_name, "_noOutliers.csv")), 
            row.names = FALSE)
  
  # ==================== STATISTICAL ANALYSIS (WITHOUT OUTLIERS) ====================
  
  results_noOUT <- data.frame(
    Region = regions,
    mean_delta = NA,
    sd_delta = NA,
    CI_low = NA,
    CI_high = NA,
    t_value = NA,
    p_value = NA
  )
  
  for (r in regions) {
    t_res <- t.test(
      data_noOutliers[[paste0(r, "_40")]], 
      data_noOutliers[[paste0(r, "_33")]], 
      paired = TRUE
    )
    
    delta_vals <- data_noOutliers[[paste0(r, "_40")]] - data_noOutliers[[paste0(r, "_33")]]
    
    results_noOUT[results_noOUT$Region == r, "mean_delta"] <- t_res$estimate
    results_noOUT[results_noOUT$Region == r, "sd_delta"] <- sd(delta_vals, na.rm = TRUE)
    results_noOUT[results_noOUT$Region == r, "CI_low"] <- t_res$conf.int[1]
    results_noOUT[results_noOUT$Region == r, "CI_high"] <- t_res$conf.int[2]
    results_noOUT[results_noOUT$Region == r, "t_value"] <- t_res$statistic
    results_noOUT[results_noOUT$Region == r, "p_value"] <- t_res$p.value
  }
  
  results_noOUT$p_FDR <- p.adjust(results_noOUT$p_value, method = "fdr")
  
  results_noOUT$significance <- cut(
    results_noOUT$p_FDR,
    breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
    labels = c("***", "**", "*", "ns")
  )
  
  write.csv(results_noOUT, 
            file.path(output_dir, paste0(metric_name, "_results_noOutliers.csv")), 
            row.names = FALSE)
  
  # ==================== LONG FORMAT ====================
  
  data_long_noOUT <- data_noOutliers %>%
    pivot_longer(
      cols = -subject,
      names_to = c("region", "timepoint"),
      names_sep = "_",
      values_to = "value"
    ) %>%
    mutate(
      timepoint = factor(timepoint, levels = c("33","40","delta")),
      region = factor(region),
      ID = subject
    )
  
  write.csv(data_long_noOUT, 
            file.path(output_dir, paste0(metric_name, "_long_noOutliers.csv")), 
            row.names = FALSE)
  
  # ==================== VISUALIZATION ====================
  
  data_delta_noOUT <- data_long_noOUT %>% filter(timepoint == "delta")
  
  region_order <- data_delta_noOUT %>%
    group_by(region) %>%
    summarise(mean_delta = mean(value, na.rm = TRUE)) %>%
    arrange(desc(mean_delta)) %>%
    pull(region)
  
  data_delta_noOUT <- data_delta_noOUT %>%
    mutate(region = factor(region, levels = region_order))
  
  # With remaining outliers labeled
  p3 <- ggplot(data_delta_noOUT, aes(x = region, y = value)) +
    geom_boxplot(fill = "orange", width = 0.6) +
    geom_text_repel(
      data = data_delta_noOUT %>%
        group_by(region) %>%
        filter(value < quantile(value, 0.25, na.rm = TRUE) - 1.5*IQR(value, na.rm = TRUE) |
                 value > quantile(value, 0.75, na.rm = TRUE) + 1.5*IQR(value, na.rm = TRUE)),
      aes(x = region, y = value, label = ID),
      size = 3,
      nudge_y = 0.1,
      segment.color = "gray50"
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
      title = paste0(metric_name, " - Delta after removing outliers"),
      x = "Region",
      y = "Delta Value"
    )
  
  ggsave(
    file.path(output_dir, paste0(metric_name, "_boxplot_noOutliers_delta.png")), 
    p3, 
    width = 12, 
    height = 6, 
    dpi = 300
  )
  
  # Without outliers displayed
  p3b <- ggplot(data_delta_noOUT, aes(x = region, y = value)) +
    geom_boxplot(fill = "orange", width = 0.6, outlier.shape = NA) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    labs(
      title = paste0(metric_name, " - Delta after removing outliers"),
      x = "Region",
      y = "Delta Value"
    )
  
  ggsave(
    file.path(output_dir, paste0(metric_name, "_boxplot_noOutliers_delta_clean.png")), 
    p3b, 
    width = 12, 
    height = 6, 
    dpi = 300
  )
  
  cat("\n✓ Outlier removal complete for", metric_name, "\n")
  cat("✓ New results saved with '_noOutliers' suffix\n\n")
  
  return(list(
    results_noOutliers = results_noOUT,
    data_noOutliers = data_noOutliers,
    data_long_noOutliers = data_long_noOUT
  ))
}

# ============================================================================
# EXAMPLE USAGE
# ============================================================================

# STEP 1: Run preprocessing for ALL metrics
# ------------------------------------------

setwd("")

metrics <- list(
list(name = "MK", file_33 = "MK_33.csv", 
     file_40 = "MK_40.csv", color = "forestgreen"),
list(name = "RK", file_33 = "RK_33.csv", 
     file_40 = "RK_40.csv", color = "#6A5ACD"),
list(name = "AK", file_33 = "AK_33.csv", 
     file_40 = "AK_40.csv", color = "#003366")
)

metrics <- list(
list(name = "BOLD", file_33 = "ICA_BOLD_SD_PT_33_all_regions.csv", 
     file_40 = "ICA_BOLD_SD_PT_40_regions_longitudinal.csv", color = "#c21e56")
)

metrics <- list(
list(name = "MD", file_33 = "DKI_MD_33.csv", 
     file_40 = "DKI_MD_40.csv", color = "#007a7a"),
list(name = "FA", file_33 = "FA-dki_33.csv", 
     file_40 = "FA-dki_40.csv", color = "#00CED1"),
list(name = "intra", file_33 = "intra_33.csv", 
     file_40 = "intra_40.csv", color = "#9BCC5F"),
list(name = "diff", file_33 = "diff_33.csv", 
     file_40 = "diff_40.csv", color = "#6D3F1F"),
list(name = "extraMD", file_33 = "extramd_33.csv", 
     file_40 = "extramd_40.csv", color = "#F08080"),
list(name = "extraTrans", file_33 = "extratrans_33.csv", 
     file_40 = "extratrans_40.csv", color = "#FFC107")
)

# Define metrics to process
metrics <- list(
  list(name = "BOLD", file_33 = "ICA_BOLD_SD_PT_33_all_regions.csv", 
       file_40 = "ICA_BOLD_SD_PT_40_regions_longitudinal.csv", color = "#c21e56"),
  list(name = "MK", file_33 = "MK_33.csv", 
       file_40 = "MK_40.csv", color = "forestgreen"),
  list(name = "RK", file_33 = "RK_33.csv", 
       file_40 = "RK_40.csv", color = "#6A5ACD"),
  list(name = "AK", file_33 = "AK_33.csv", 
       file_40 = "AK_40.csv", color = "#003366"),
  list(name = "MD", file_33 = "DKI_MD_33.csv", 
       file_40 = "DKI_MD_40.csv", color = "#007a7a"),
  list(name = "FA", file_33 = "FA-dki_33.csv", 
       file_40 = "FA-dki_40.csv", color = "#00CED1"),
  list(name = "intra", file_33 = "intra_33.csv", 
       file_40 = "intra_40.csv", color = "#9BCC5F"),
  list(name = "diff", file_33 = "diff_33.csv", 
       file_40 = "diff_40.csv", color = "#6D3F1F"),
  list(name = "extraMD", file_33 = "extramd_33.csv", 
       file_40 = "extramd_40.csv", color = "#F08080"),
  list(name = "extraTrans", file_33 = "extratrans_33.csv", 
       file_40 = "extratrans_40.csv", color = "#FFC107")
)

# Process each metric
for (metric in metrics) {
  output_dir <- metric$name  # Results will be saved in the same folder
  
  process_metric_data(
    data_file_33 = file.path(metric$name, metric$file_33),
    data_file_40 = file.path(metric$name, metric$file_40),
    metric_name = metric$name,
    output_dir = file.path(output_dir, "Long_noMotion")
  )
}

cat("\n========================================\n")
cat("✓ All metrics processed!\n")
cat("✓ Review the outlier CSV files in each folder\n")
cat("✓ Then proceed to STEP 2 to remove outliers\n")
cat("========================================\n\n")


# STEP 2: Review outlier tables and define which to remove
# ---------------------------------------------------------
# After reviewing the outlier CSV files in each metric folder define your outlier lists for each metric - defined below as calculated before

BOLD_outliers_delta <- list(
  PFC      = c("24", "39", "42", "43", "44", "65"),
  limbic   = c("44", "51", "65"),
  visual   = c("49", "50", "63", "69"),
  PCC      = c("50", "63", "69"),
  thalamus = c("50", "51")
)

MK_outliers_delta <- list(
  limbic   = c("35"),
  SSM      = c("39", "41"),
  PFC      = c("35"),
  thalamus = c("4"),
  PCC      = c("41"),
  visual   = c("41")
)

RK_outliers_delta <- list(
  PCUN     = c("4", "35", "38"),
  SSM      = c("4","39", "41"),
  auditory   = c("4"),
  limbic = c("35"),
  PCC      = c("41"),
  PFC      = c("35"),
  thalamus = c("4"),
  visual   = c("41")
)

AK_outliers_delta <- list(
  SSM      = c("39", "41"),
  auditory   = c("41"),
  PCC      = c("41"),
  PCUN     = c("35"),
  visual   = c("41"),
  PFC      = c("35")
)

MD_outliers_delta <- list(
  auditory   = c("13", "28", "42"),
  paralimbic = c("32", "33", "42"),
  visual   = c("41", "44"),
  limbic = c("17"),
  PCC      = c("32"),
  SSM      = c("41"),
  thalamus = c("4")
)

FA_outliers_delta <- list(
  paralimbic = c("4", "18", "37", "51"),
  PCUN     = c("4", "37", "38", "51"),
  limbic   = c("4", "37", "42"),
  PCC      = c("4", "38", "41"),
  SSM      = c("4", "41", "51"),
  visual   = c("4", "41", "51"),
  auditory   = c("4", "37"),
  PFC      = c("4"),
  thalamus = c("4")
)

intra_outliers_delta <- list(
  auditory = c("4", "17", "41"),
  PCC      = c("4", "30", "45"),
  SSM      = c("4", "39", "41"),
  limbic   = c("4", "41"),
  paralimbic = c("4"),
  PCUN = c("4"),
  thalamus = c("4"),
  visual   = c("4")
)

diff_outliers_delta <- list(
  limbic      = c("4", "28", "33", "35"),
  PCC         = c("4", "32", "45"),
  visual      = c("4", "32", "45" ),
  paralimbic  = c("4", "35"),
  PFC         = c("4", "28"),
  auditory    = c("4"),
  PCUN        = c("4"),
  SSM         = c("4"),
  thalamus = c("4")
)

extraMD_outliers_delta <- list(
  limbic      = c("4", "24", "28", "33","35"),
  paralimbic  = c("4", "35"),
  visual      = c("4", "32", "44", "45"),
  PCC         = c("4", "32"),
  auditory    = c("4"),
  PCUN        = c("4"),
  PFC         = c("4"),
  SSM         = c("4"),
  thalamus = c("4")
)

extraTrans_outliers_delta <- list(
  limbic      = c("4", "28", "33"),
  visual      = c("32", "44", "45"),
  paralimbic  = c("4", "33"),
  PFC = c("4"),
  thalamus = c("4")
)



# STEP 3: Remove outliers and reanalyze for each metric
# ------------------------------------------------------

# BOLD
result_BOLD_clean <- remove_outliers_and_reanalyze(
  metric_name = "BOLD",
  output_dir = "BOLD",
  outlier_list_delta = BOLD_outliers_delta
)
# DKI
# MK
result_MK_clean <- remove_outliers_and_reanalyze(
  metric_name = "MK",
  output_dir = file.path("MK", "Long_noMotion"),
  outlier_list_delta = MK_outliers_delta
)
# RK
result_RK_clean <- remove_outliers_and_reanalyze(
  metric_name = "RK",
  output_dir = file.path("RK", "Long_noMotion"),
  outlier_list_delta = RK_outliers_delta
)
# AK
result_AK_clean <- remove_outliers_and_reanalyze(
  metric_name = "AK",
  output_dir = file.path("AK", "Long_noMotion"),
  outlier_list_delta = AK_outliers_delta
)
# DTI
# MD
result_MD_clean <- remove_outliers_and_reanalyze(
  metric_name = "MD",
  output_dir = file.path("MD", "Long_noMotion"),
  outlier_list_delta = MD_outliers_delta
)
# FA
result_FA_clean <- remove_outliers_and_reanalyze(
  metric_name = "FA",
  output_dir = file.path("FA", "Long_noMotion"),
  outlier_list_delta = FA_outliers_delta
)
# SMT
# intra
result_intra_clean <- remove_outliers_and_reanalyze(
  metric_name = "intra",
  output_dir = file.path("intra", "Long_noMotion"),
  outlier_list_delta = intra_outliers_delta
)
# diff
result_diff_clean <- remove_outliers_and_reanalyze(
  metric_name = "diff",
  output_dir = file.path("diff", "Long_noMotion"),
  outlier_list_delta = diff_outliers_delta
)
# extraMD
result_extraMD_clean <- remove_outliers_and_reanalyze(
  metric_name = "extraMD",
  output_dir = file.path("extraMD", "Long_noMotion"),
  outlier_list_delta = extraMD_outliers_delta
)
# extraTrans
result_extraTrans_clean <- remove_outliers_and_reanalyze(
  metric_name = "extraTrans",
  output_dir = file.path("extraTrans", "Long_noMotion"),
  outlier_list_delta = extraTrans_outliers_delta
)

cat("\n========================================\n")
cat("✓ All outlier removal complete!\n")
cat("✓ Check '_noOutliers' files in each folder\n")
cat("========================================\n\n")