# Load required packages
library(lme4)
library(emmeans)
library(dplyr)
library(lmerTest)

# Get the current directory of the script (or use getwd() for the working directory)
code_dir <- dirname(rstudioapi::getSourceEditorContext()$path)

# Construct the path to the CSV file
tablePath <- file.path(code_dir, "..", "output", "Quantitative_outputs", "Quantitative_results_from_dwi_processing.csv")

# Read the CSV file
data <- read.csv(tablePath)

# Ensure the relevant variables are treated as factors
data$merged_timepoint <- as.factor(data$merged_timepoint)  # Time points
data$Group <- as.factor(data$Group)  # Stroke vs Control
data$mask_name <- as.factor(data$mask_name)  # Brain regions
data$subjectID <- as.factor(data$subjectID)  # Mouse ID
data$dialation_amount <- as.factor(data$dialation_amount)  # Ensure dialation_amount is a factor

# Filter the data for FA type only
fa_data <- subset(data, Qtype == "fa")

# Check how many rows are left after filtering
cat("Number of rows after filtering for FA data:", nrow(fa_data), "\n")

# Check for missing data in the relevant columns
missing_data <- colSums(is.na(fa_data[c("Value", "merged_timepoint", "Group", "mask_name", "dialation_amount", "subjectID")]))
cat("Missing values in each relevant column:\n")
print(missing_data)

# If there are missing values, remove rows with missing data
fa_data_clean <- na.omit(fa_data[c("Value", "merged_timepoint", "Group", "mask_name", "dialation_amount", "subjectID")])

# Check how many rows remain after removing missing data
cat("Number of rows after removing missing values:", nrow(fa_data_clean), "\n")

# Check the levels of merged_timepoint to ensure reference is correct
print(levels(fa_data_clean$merged_timepoint))

# Optionally, set the reference level for merged_timepoint if needed
fa_data_clean$merged_timepoint <- relevel(fa_data_clean$merged_timepoint, ref = "0")

# Optionally, set the reference level for Group if needed
fa_data_clean$Group <- relevel(fa_data_clean$Group, ref = "Sham")

# Fit the linear mixed-effects model using the cleaned dataset
# Simplify the random effects structure to reduce boundary singular fit issues
fa_model <- lmer(Value ~ merged_timepoint * Group + mask_name + dialation_amount + 
                   (1 | subjectID) + (1 | mask_name), 
                 data = fa_data_clean)

# Summarize the model
cat("Model Summary:\n")
summary(fa_model)

# Set emm_options to allow degrees of freedom calculations for large datasets and increase rg.limit to avoid issues with reference grid
emm_options(pbkrtest.limit = 342163, lmerTest.limit = 342163, rg.limit = 12000)

# Post-hoc comparisons for significant differences between Stroke and Control within each time point for each mask and dilation amount
posthoc_group <- emmeans(fa_model, pairwise ~ Group | merged_timepoint * mask_name * dialation_amount)

# Post-hoc comparisons for significant differences between time points within each group, mask, and dilation amount
posthoc_time <- emmeans(fa_model, pairwise ~ merged_timepoint | Group * mask_name * dialation_amount)

# Construct the path to the output text file (same location as CSV file)
output_file <- file.path(dirname(tablePath), "posthoc_results_FA_model2.txt")

# Open a connection to the file
sink(output_file)

# Display the results of the pairwise comparisons

# Question 1: Are there significant differences between Stroke and Control within each time point?
cat("Post-hoc comparison: Stroke vs Control within each time point, mask, and dilation amount\n")
print(summary(posthoc_group))

# Question 2: Are there significant differences between time points within each group?
cat("\nPost-hoc comparison: Time point differences within each group, mask, and dilation amount\n")
print(summary(posthoc_time))

# Print pairwise results for Stroke vs Control within each time point, mask, and dilation amount
cat("\nPairwise results for Stroke vs Control within each time point, mask, and dilation amount:\n")
print(posthoc_group$contrasts)

# Print pairwise results for time point differences within each group, mask, and dilation amount
cat("\nPairwise results for time point differences within each group, mask, and dilation amount:\n")
print(posthoc_time$contrasts)

# Close the connection to the file
sink()

# Notify user that the results have been saved
cat("Results have been saved to", output_file, "\n")