# Ensure necessary packages are loaded
library(dplyr)
library(ggplot2)
library(gridExtra)

# Prompt the user to select a file (LongHFScores.csv <- available in the Supplementary Materials)
selected_file <- file.choose()

# Load the dataframe file (LongHFScores.csv <- available in the Supplementary Materials)
data <- read.csv(selected_file)

# Function to check normality, apply the appropriate test, and plot histograms with density curves
perform_test_and_plot <- function(data, var, condition) {
  # Filter the data for each level of SOCIAL_ISOLATION
  before <- data %>% filter(SOCIAL_ISOLATION == 0) %>% pull(var)
  during <- data %>% filter(SOCIAL_ISOLATION == 1) %>% pull(var)
  
  # Check normality
  shapiro_before <- shapiro.test(before)$p.value
  shapiro_during <- shapiro.test(during)$p.value
  
  # Print normality test results
  cat("Shapiro-Wilk normality test for", var, "(Before isolation): p-value =", shapiro_before, "\n")
  cat("Shapiro-Wilk normality test for", var, "(During isolation): p-value =", shapiro_during, "\n")
  
  # Select appropriate test
  if (shapiro_before > 0.05 & shapiro_during > 0.05) {
    # If both groups are normal, use the paired Student's t-test
    t_test <- t.test(before, during, paired = TRUE)
    cat("Paired t-test for", var, ": p-value =", t_test$p.value, "\n")
  } else {
    # Otherwise, use the paired Wilcoxon test
    wilcox_test <- wilcox.test(before, during, paired = TRUE)
    cat("Paired Wilcoxon test for", var, ": p-value =", wilcox_test$p.value, "\n")
  }
  
  # Set number of bins for histograms
  num_bins <- 8  # Adjust for between 7 to 10 intervals

  # Plot histograms with density curve
  plot_before <- ggplot(data.frame(value = before), aes(x = value)) +
    geom_histogram(aes(y = after_stat(density)), bins = num_bins, fill = "blue", alpha = 0.4) +
    geom_density(alpha = 0.5, fill = "blue") +
    ggtitle(paste("Before -", var, "(Condition =", condition, ")")) +
    xlab(var) +
    ylab("Density")

  plot_during <- ggplot(data.frame(value = during), aes(x = value)) +
    geom_histogram(aes(y = after_stat(density)), bins = num_bins, fill = "red", alpha = 0.4) +
    geom_density(alpha = 0.5, fill = "red") +
    ggtitle(paste("During -", var, "(Condition =", condition, ")")) +
    xlab(var) +
    ylab("Density")
  
  return(list(plot_before, plot_during))
}

# Filtered data for each condition
data_adhd <- data %>% filter(CONDITION == 1)
data_asd <- data %>% filter(CONDITION == 2)
data_adhd_asd <- data %>% filter(CONDITION == 3)

# Create a list to store all plots
all_plots <- list()

# 1. SPECIFIC_SCENARIOS
cat("\nComparison for SPECIFIC_SCENARIOS:\n")
cat("Condition 1 (ADHD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd, "SPECIFIC_SCENARIOS", 1))

cat("Condition 2 (ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_asd, "SPECIFIC_SCENARIOS", 2))

cat("Condition 3 (ADHD & ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd_asd, "SPECIFIC_SCENARIOS", 3))

# 2. SCREEN_TIME
cat("\nComparison for SCREEN_TIME:\n")
cat("Condition 1 (ADHD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd, "SCREEN_TIME", 1))

cat("Condition 2 (ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_asd, "SCREEN_TIME", 2))

cat("Condition 3 (ADHD & ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd_asd, "SCREEN_TIME", 3))

# 3. TOTAL_SCORE
cat("\nComparison for TOTAL_SCORE:\n")
cat("Condition 1 (ADHD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd, "TOTAL_SCORE", 1))

cat("Condition 2 (ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_asd, "TOTAL_SCORE", 2))

cat("Condition 3 (ADHD & ASD):\n")
all_plots <- c(all_plots, perform_test_and_plot(data_adhd_asd, "TOTAL_SCORE", 3))

# Display all plots in one window
do.call(grid.arrange, c(all_plots, ncol = 3))