blocks <- c("B1", "B2", "B3", "B4", "B5")
blocks <- rep(blocks, each=104)

sequence = c("B1", "B2", "B3", "B4", "B5")
Period=rep(sequence, each = 104) 

flankerData$Block <- rep(Period, times=60)

flankerData$Block <- factor(flankerData$Block)
flankerData$Stim <- factor(flankerData$Stim)
flankerData<-flankerData %>% mutate(Stim= recode(Stim, `A` = "Anodal", `B` = "Sham", `C` = "Cathodal"))
flankerData$Subject <-  factor(flankerData$Subject)

flankerData <- filter(flankerData, RT > 150) #All stimuli 

flankerData_filtered <- filter(flankerData, Error == 1)
# Simplify dataset. Calculate one RT per participant and block

flankerData_averaged <- flankerData_filtered %>% 
    group_by(Subject, Stim, Block, Comp) %>%
  summarize(RT = mean(RT, na.rm = TRUE))


# Full model with all fixed effects and interactions
model_full <- lmer(RT ~ Stim * Block * Comp  + (1 | Subject), data = flankerData_averaged)

# Model without interaction terms
model_main_effects <- lmer(RT ~ Stim + Block +  (1 | Subject), data = flankerData_averaged)

# Null model (no fixed effects)
model_null <- lmer(RT ~ 1 + (1 | Subject), data = flankerData_averaged)

anova(model_null, model_main_effects,  model_full)

# anova(model_full)
anova_results <- anova(model_full)
anova_table <- as.data.frame(anova_results)
print(anova_table)

# eta_squared(model_full)

Comp <- flankerData_averaged %>% 
    group_by(Comp) %>%
  summarize(RT = mean(RT, na.rm = TRUE))

ggplot(flankerData_averaged, aes(x = Block, y = RT, fill= Comp)) +
  geom_boxplot() +
  labs(title = "Distribution of RT across blocks and stimuli type",
       x = "Block",
       y = "Reaction times (ms)") +
  theme_minimal()

#For ACC
# Calculate the % of correct answer for each block


ACC <- flankerData %>%
  group_by(Subject, Block, Stim, Comp) %>%
  summarize(ACC = mean(Error == 1) * 100)


# Full model with all fixed effects and interactions
model_full <- lmer(ACC ~ Stim * Block *Comp  + (1 | Subject), data = ACC)

# Model without interaction terms
model_main_effects <- lmer(ACC ~ Stim + Block + Comp +  (1 | Subject), data = ACC)

# Null model (no fixed effects)
model_null <- lmer(ACC ~ 1 + (1 | Subject), data = ACC)

anova(model_null, model_main_effects,  model_full)

# eta_squared(model_full)
anova_results <- anova(model_full)
anova_table <- as.data.frame(anova_results)
print(anova_table)

Comp_ACC <- ACC %>% 
    group_by(Comp) %>%
  summarize(ACC = mean(ACC, na.rm = TRUE))

ggplot(ACC, aes(x = Block, y = ACC, fill= Comp)) +
  geom_boxplot() +
  labs(title = "Distribution of correct respons across blocks and stimuli type",
       x = "Block",
       y = "ACC (%)") +
  theme_minimal()

# Create a raincloud plot

  # Resting
delta_plot <-  ggplot(delta, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = " Normalized score (AU)", title = "Delta") +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())


theta_plot <-  ggplot(theta, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "Normalized score (AU)", title = "Theta") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
  


alpha_plot <-  ggplot(alpha, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = " Normalized score (AU)", title = "Alpha") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
  


beta_plot <-  ggplot(beta, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(y = " Normalized score (AU)", title = "Beta") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())


gamma_plot <-  ggplot(gamma, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(y = " Normalized score (AU)", title = "Gamma") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())


p1 <- ggarrange(delta_plot, theta_plot, alpha_plot, beta_plot, gamma_plot,
          # labels = c("delta", "theta", "alpha", "beta", "gamma"),
          ncol = 3, nrow = 2, 
          # common.legend = T,
          # legend = "right",
          align = "h")

check_outliers(
  delta$normalized, 
  method = c("ci", "zscore_robust", "iqr"),
  threshold = NULL, #Default threshold like this. If you want to specify, use: list('zscore' = 2, 'iqr' = 1.97))
  ID = NULL,
  verbose = TRUE
)

# Task

delta_plot_task <-  ggplot(delta_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "Normalized score (AU)", title = "Delta") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())


theta_plot_task <-  ggplot(theta_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "Normalized score (AU)", title = "Theta") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())



alpha_plot_task <-  ggplot(alpha_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "Normalized score (AU)", title = "Alpha") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())



beta_plot_task <-  ggplot(beta_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "postNormalized score (AU)", title = "Beta") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())


gamma_plot_task <-  ggplot(gamma_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  labs(x = "Condition", y = "Normalized score (AU)", title = "Gamma") +
 theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5))+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

p2 <- ggarrange(delta_plot_task,theta_plot_task, alpha_plot_task, beta_plot_task, gamma_plot_task, 
          # labels = c("delta", "theta", "alpha", "beta", "gamma"),
          ncol = 3, nrow = 2, 
          common.legend = T,
          legend = "right",
          align = "h")

# Add labels to the combined plots
labeled_combined_plot1 <- plot_grid(p1, labels = "Resting", label_size = 20, ncol = 1)
labeled_combined_plot2 <- plot_grid(p2, labels = "Task", label_size = 20, ncol = 1)

final_combined_plot <- plot_grid(labeled_combined_plot1, labeled_combined_plot2, ncol = 1)


# grid.arrange(p1, p2, ncol = 1)

lzc_resting[, 3:300] <- lapply(lzc_resting[, 3:300], as.numeric) #transform the data point in numeric
lzc_resting$Condition <- factor(lzc_resting$Condition)

lzc_resting$pre <- rowMeans(lzc_resting[, 3:151], na.rm = TRUE)
lzc_resting$post <- rowMeans(lzc_resting[, 152:300],na.rm = TRUE)
lzc_resting <- lzc_resting[,c(1:2, 301:302)] # Simplify the data set for clarity

lzsum_resting[, 3:300] <- lapply(lzsum_resting[, 3:300], as.numeric)
lzsum_resting$Condition <- factor(lzsum_resting$Condition)
lzsum_resting$pre <- rowMeans(lzsum_resting[, 3:151], na.rm = TRUE)
lzsum_resting$post <- rowMeans(lzsum_resting[, 152:300],na.rm = TRUE)
lzsum_resting <- lzsum_resting[,c(1:2, 301:302)] # Simplify the data set for clarity


# Calculate ratios
lzc_resting <- lzc_resting %>%
  mutate(post_ratio = post / (pre + post))
lzc_resting <- lzc_resting %>%
  mutate(diff = (post - pre) / pre*100)
lzc_resting <- lzc_resting %>%
  mutate(normalized = (post - pre) / (post + pre))

lzsum_resting <- lzsum_resting %>%
  mutate(post_ratio = post / (pre + post))
lzsum_resting <- lzsum_resting %>%
  mutate(diff = (post - pre) / pre*100)
lzsum_resting <- lzsum_resting %>%
  mutate(normalized = (post - pre) / (post + pre))

lzc_plot_resting <-  ggplot(lzc_resting, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  scale_y_continuous(limits = c(-0.2, 0.4))+ 
  labs(x = FALSE, y =  "Normalized score (AU)") +
  theme(legend.position = "none")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# lzc_plot_resting


lzsum_plot_resting <-  ggplot(lzsum_resting, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  scale_y_continuous(limits = c(-0.2, 0.3))+ 
  labs(y = "Normalized score (AU)") +
  theme(legend.position = "none")+
  easy_remove_x_axis(what = "text")+
   theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# lzsum_plot_resting

ggarrange(lzc_plot_resting, lzsum_plot_resting, 
          labels = c("Global LZc", "Global LZ sum"),
          ncol = 2, nrow = 1, 
          common.legend = T,
          legend = "right",
          align = "h")

lzc_task[, 3:360] <- lapply(lzc_task[, 3:360], as.numeric) #transform the data point in numeric
lzc_task$Condition <- factor(lzc_task$Condition)

lzc_task$pre <- rowMeans(lzc_task[, 3:181], na.rm = TRUE)
lzc_task$post <- rowMeans(lzc_task[, 182:360],na.rm = TRUE)
lzc_task <- lzc_task[,c(1:2, 361:362)] # Simplify the data set for clarity

lzsum_task[, 3:360] <- lapply(lzsum_task[, 3:360], as.numeric)
lzsum_task$Condition <- factor(lzsum_task$Condition)
lzsum_task$pre <- rowMeans(lzsum_task[, 3:181], na.rm = TRUE)
lzsum_task$post <- rowMeans(lzsum_task[, 182:360],na.rm = TRUE)
lzsum_task <- lzsum_task[,c(1:2, 361:362)] # Simplify the data set for clarity


# Calculate ratios
lzc_task <- lzc_task %>%
  mutate(post_ratio = post / (pre + post))
lzc_task <- lzc_task %>%
  mutate(diff = (post - pre) / pre*100)
lzc_task <- lzc_task %>%
  mutate(normalized = (post - pre) / (post+pre))

lzsum_task <- lzsum_task %>%
  mutate(post_ratio = post / (pre + post))
lzsum_task <- lzsum_task %>%
  mutate(diff = (post - pre) / pre*100)
lzsum_task <- lzsum_task %>%
  mutate(normalized = (post - pre) / (post+pre))

lzc_plot_task <-  ggplot(lzc_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  scale_y_continuous(limits = c(-0.2, 0.3))+ 
  labs(x = "Condition", y = "Normalized Score (AU)") +
  theme(legend.position = "none")+
   theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# lzc_plot_task

lzsum_plot_task <-  ggplot(lzsum_task, aes(x = Condition, y = normalized, fill = Condition)) +
  geom_flat_violin(position = position_nudge(x = 0.2), alpha = 0.7) +
  geom_point(position = position_jitter(width = 0.2, height = 0), alpha = 0.3, size = 2) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.5) +
  theme_minimal() +
  scale_y_continuous(limits = c(-0.2, 0.3))+ 
  labs(x = "Condition", y = "Normalized Score (AU)") +
  theme(legend.position = "none")+
   theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# lzsum_plot_task

ggarrange(lzc_plot_task, lzsum_plot_task, 
          labels = c("Global LZc", "Global LZ sum"),
          ncol = 2, nrow = 1, 
          common.legend = T,
          legend = "right",
          align = "h")

data_wide <- data %>% 
  pivot_wider(
    names_from  = c(Dimension), # Can accommodate more variables, if needed.
    values_from = c(Value)
  )

knitr::kable(head(data_wide[, 1:10]), "pipe", caption = "Data organization for the clustering")

data_cluster<- data_wide[,c(1:5, 7:11)] #Remove the Discomfort and Stimulation variables which are not related 


fviz_nbclust(data_cluster[, 4:10], kmeans, method = "wss")+
  labs(subtitle="Elbow Method")

fviz_nbclust(data_cluster[, 4:10], # data
             kmeans, # clustering algorithm
             method = "silhouette")+ # silhouette
  labs(subtitle="Silhouette Method")

# fviz_nbclust(data_wide[, 4:9], # data, 
#              kmeans ,
#              nstart = 25, 
#              method = "gap_stat")+
#     labs(subtitle="Gap Methods")

# Load necessary libraries
library(stats)
library(clusterCrit)



# Function to calculate the CH score for a given number of clusters
calculate_ch_score <- function(data, k) {
  set.seed(123)
  kmeans_result <- kmeans(data, centers = k)
  clustering_vector <- kmeans_result$cluster
  ch_score <- intCriteria(as.matrix(data), clustering_vector, "Calinski_Harabasz")
  return(ch_score$calinski_harabasz)
}

# Range of clusters to evaluate
k_values <- 2:10

# Calculate CH scores for each number of clusters
ch_scores <- sapply(k_values, calculate_ch_score, data = data_cluster[, 4:10])

# Print CH scores for each k
for (i in 1:length(k_values)) {
  cat("k =", k_values[i], ": CH Score =", ch_scores[i], "\n")
}

# Plot the CH scores to visualize the optimal number of clusters
plot(k_values, ch_scores, type = "b", pch = 19, frame = FALSE,
     xlab = "Number of Clusters (k)", ylab = "Caliński-Harabasz Score",
     main = "Caliński-Harabasz Score vs. Number of Clusters")

# Optional: Highlight the optimal k
optimal_k <- k_values[which.max(ch_scores)]
points(optimal_k, max(ch_scores), col = "red", pch = 19)
text(optimal_k, max(ch_scores), labels = paste("Optimal k =", optimal_k), pos = 4, col = "red")

km.out <- kmeans(data_cluster[, 4:10], centers=2,nstart=100)
km.out$centers

#  Total within-cluster sum of squares
# km.out$tot.withinss
# km.clusters<-km.out$cluster
data_cluster$Cluster <-km.out$cluster


# fviz_cluster(list(data=data_wide[, 4:9], cluster = km.clusters))
# 
# cluster_stats <- cluster.stats(data_wide[, 4:9], clustering = km.clusters)
# 
# ch_index <- cluster.stats(data_wide[, 4:9], clustering = km.clusters, criterion = "ch")

#Change Epoch time into 3 different phases of (pre-stim, stimulation and post-stim)
segment_lengths <- c(12, 36, 12)

# Create a vector to specify the period for each observation
period <- rep(1:length(segment_lengths), times = segment_lengths)
data_cluster$Phase <- rep(period, 60)
data_cluster<-data_cluster %>% mutate(Phase= recode(Phase, `1` = "Pre", `2` = "During", `3` = "Post"))
data_cluster$Phase <- factor(data_cluster$Phase, levels=c("Pre", "During", "Post"))

# Transform data to long format with the cluster assignation. We are going to create data frame with Participant, Epoch, Dimension, Cluster, Value(intensity)

Cluster_long<- pivot_longer(data_cluster, cols = c(Alertness, Boredom, Effort, Execution, Fedup, Mindwandering, Stimulation), names_to = "Dimension",
                        values_to = "Value", values_drop_na = TRUE)

Cluster_long$Cluster <- as.factor(Cluster_long$Cluster)


sumrepdatClust <- summarySE(Cluster_long, measurevar = "Value",
                       groupvars=c("Dimension", "Cluster"))

ClusterRain <- ggplot(Cluster_long, aes(x = Dimension, y = Value, fill = Cluster)) +
  geom_flat_violin(aes(fill = Cluster),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  # geom_point(aes(x = Dimension, y = Value, colour = Cluster),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Dimension, y = Value, fill = Cluster),outlier.shape = NA, alpha = .5, width = .3, colour = "black")+
  # geom_line(data = sumrepdatClust, aes(x = Dimension, y = Value, group = Cluster, colour = Cluster), linetype = 3)+
  geom_point(data = sumrepdatClust, aes(x = Dimension, y = Value, group = Cluster, colour = Cluster), shape = 18, size = 3) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "Intensity", x = element_blank(),
        title ="Intensity of each dimension for each cluster",
        )+
    theme(legend.position="right")+
     scale_x_discrete(guide = guide_axis(n.dodge=3))


ClusterRain

# If we want to divide by conditions

# ggplot(Cluster_long, aes(x = Dimension, y = Value, fill = Cluster)) +
#   # geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = 0.8) +
#   # geom_point(position = position_jitter(width = .15), size = 1, color = "black", alpha = 0.6) +
#   geom_boxplot(width = 0.2, outlier.shape = NA, alpha = 0.5) +
#   scale_fill_manual(values = c("#00AFBB", "#E7B800")) +  # Choose your colors
#   # theme_minimal() +
#   labs(x = "Dimension", y = "Value", fill = "Cluster") +
#   theme(legend.position = "right") + facet_grid(~Condition)

data_relative_cluster <- summarySE(Cluster_long, measurevar = "Value",
                       groupvars=c("Participant", "Dimension", "Cluster"))

data_relative_cluster <- data_relative_cluster[, c(1:3,5)]


# Define the dimensions to loop over
dimensions <- c("Alertness", "Boredom", "Effort", "Execution", "Fedup", "Mindwandering", "Stimulation")

# Initialize a list to store model results
model_results <- list()

for (dim in dimensions) {
    # Subset the data for the current dimension
    data_subset <- subset(data_relative_cluster, Dimension == dim)
    
    # Fit the linear mixed-effects model
    model <- lmer(Value ~ Cluster + (1 | Participant), data = data_subset)
    
    # Extract and tidy the fixed effects summary
    model_summary <- tidy(model, effects = "fixed") %>%
        mutate(Dimension = dim) %>% 
      mutate(across(c(estimate, std.error, statistic, df, p.value), ~ round(., 2)))
    
    # Append the results to the list
    model_results[[dim]] <- model_summary
}

# Combine all model results into a single data frame
model_results_df <- bind_rows(model_results)

# Combine all model results into a single data frame
model_results_df <- bind_rows(model_results)

# Create a table for model results
kable(model_results_df, format = "html", table.attr = "class='table table-striped'") %>%
  kable_styling(full_width = F, position = "left", bootstrap_options = c("striped", "hover", "condensed"))

# Descriptive 
descriptive_stats <- list()

for (dim in dimensions) {
    # Subset the data for the current dimension
    data_subset <- subset(data_relative_cluster, Dimension == dim)
    
    # Calculate descriptive statistics by Cluster
    descriptive_by_cluster <- data_subset %>%
        group_by(Cluster) %>%
        summarise(
            count = n(),
            mean = mean(Value, na.rm = TRUE),
            sd = sd(Value, na.rm = TRUE),
            se = sd / sqrt(count),
            ci_lower = mean - qt(1 - (0.05 / 2), count - 1) * se,
            ci_upper = mean + qt(1 - (0.05 / 2), count - 1) * se
        ) %>%
        mutate(Dimension = dim)
    
    # Append the results to the list
    descriptive_stats[[dim]] <- descriptive_by_cluster
}

# Combine all descriptive statistics into a single data frame
descriptive_stats_df <- bind_rows(descriptive_stats)

kable(descriptive_stats_df, format = "html", table.attr = "class='table table-striped'") %>%
  kable_styling(full_width = F, position = "left", bootstrap_options = c("striped", "hover", "condensed"))

# for t-test

# Initialize a list to store t-test results
t_test_results <- list()

for (dim in dimensions) {
    # Subset the data for the current dimension
    data_subset <- subset(data_relative_cluster, Dimension == dim)
    
    # Perform independent t-test between clusters
    t_test <- t.test(Value ~ Cluster, data = data_subset)
    
    # Create a summary data frame for t-test results
    t_test_summary <- data.frame(
        Dimension = dim,
        t_statistic = t_test$statistic,
        p_value = t_test$p.value,
        conf_int_lower = t_test$conf.int[1],
        conf_int_upper = t_test$conf.int[2],
        mean_diff = t_test$estimate[1] - t_test$estimate[2]
    )
    
    # Append the results to the list
    t_test_results[[dim]] <- t_test_summary
}

# Combine all t-test results into a single data frame
t_test_results_df <- bind_rows(t_test_results)

# Create a table for t-test results
kable(t_test_results_df, format = "html", table.attr = "class='table table-striped'") %>%
  kable_styling(full_width = F, position = "left", bootstrap_options = c("striped", "hover", "condensed"))

# Split data into participant condition phase and cluster to obtain the relative frequency 

relative_complete_phases <-  summarySE(Cluster_long, measurevar = "Value",
                       groupvars=c("Participant", "Condition", "Phase", "Cluster"))
relative_complete_phases <- relative_complete_phases[,1:5]



# To fill the NA with 0 to include phases with no representation of clusters
relative_complete_phases <- relative_complete_phases %>%
  complete(Participant, Cluster, Phase, fill = list(N = 0))

# After including the 0 in the N, it does not fill in the conditions, so
# thus we need to fill Condition based on the closest match:
relative_complete_phases <- relative_complete_phases %>%
  group_by(Participant) %>%
  fill(Condition, .direction = "downup") %>%
  ungroup()
relative_complete_phases$Cluster <- as.factor(relative_complete_phases$Cluster)

# Calculate the relative frequency for each cluster by participant and phase
relative_complete_phases <- relative_complete_phases %>%
  group_by(Participant, Phase) %>%
  mutate(Frequency = N / sum(N)) %>%
  ungroup()

# filter each condition to plot the relative frequency

relative_anodal <- relative_complete_phases %>% filter(Condition == "Anodal") 
relative_cathodal <- relative_complete_phases %>% filter(Condition == "Cathodal")
relative_sham <- relative_complete_phases %>% filter(Condition == "Sham")

# To represent only cluster 1

# relative_anodal <- relative_complete_phases %>% filter(Condition == "Anodal" & Cluster == "1") 
# relative_cathodal <- relative_complete_phases %>% filter(Condition == "Cathodal" & Cluster == "1")
# relative_sham <- relative_complete_phases %>% filter(Condition == "Sham"& Cluster == "1")



 relat_anodal <- ggplot(relative_anodal, aes(x = Phase, y = Frequency, fill = Cluster)) +
  geom_flat_violin(aes(fill = Cluster),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_point(aes(x = as.numeric(Phase)-.15, y = Frequency, colour = Cluster),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Phase, y = Frequency, fill = Cluster),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  # geom_line(data = sumrepdatAVG20, aes(x = as.numeric(Phase)+.1, y = Value, group = Cluster, colour = Condition), linetype = 3)+
  # geom_point(data = sumrepdatAVG20, aes(x = as.numeric(Phase)+.1, y = Value, group = Cluster, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
sm_classic(legends = FALSE)+
      labs( y = "Relative frequency", x = element_blank(),
        title ="Anodal")

 
 relat_cathodal <- ggplot(relative_cathodal, aes(x = Phase, y = Frequency, fill = Cluster)) +
  geom_flat_violin(aes(fill = Cluster),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_point(aes(x = as.numeric(Phase)-.15, y = Frequency, colour = Cluster),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Phase, y = Frequency, fill = Cluster),outlier.shape = NA, alpha = .5, width = .6, colour = "black")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
sm_classic(legends = FALSE)+
  # easy_remove_x_axis(what = "text")+
   labs( y = element_blank(), x = element_blank(),
        title ="Cathodal")

 
relat_sham <- ggplot(relative_sham, aes(x = Phase, y = Frequency, fill = Cluster)) +
  geom_flat_violin(aes(fill = Cluster),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_point(aes(x = as.numeric(Phase)-.15, y = Frequency, colour = Cluster),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Phase, y = Frequency, fill = Cluster),outlier.shape = NA, alpha = .5, width = .6, colour = "black")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic(legends = FALSE)+
  # easy_remove_x_axis(what = "text")+
  labs( y = element_blank(), x = element_blank(),
        title ="Sham")


Relat_plots <- ggarrange(relat_anodal, relat_cathodal, relat_sham,
                      ncol =3 ,labels = c("A", "B", "C"),
                      nrow=1,
                      common.legend = T,
                      legend = "right",
                      align = "h")

Relat_plots

descriptive_relative <- summarySE(relative_complete_phases, measurevar = "Frequency", groupvars=c("Condition", "Phase", "Cluster"))


 descriptive_relative_2 <- relative_complete_phases %>% 
  group_by(Condition, Phase, Cluster) %>%
  summarise(average = median(Frequency))


kable(descriptive_relative, format = "html", table.attr = "class='table table-striped'") %>%
  kable_styling(full_width = F, position = "left", bootstrap_options = c("striped", "hover", "condensed"))

# Calculate the Pre value of each point
relative_complete_phases <- relative_complete_phases %>%
  group_by(Participant, Cluster, Condition) %>%
  mutate(Pre_value = Frequency[Phase == "Pre"]) %>%
  ungroup()


# Full model with all fixed effects and interactions
model_full <- lmer(Frequency ~ Cluster * Phase * Condition + Pre_value + (1 | Participant), data = relative_complete_phases)

# Model without interaction terms
model_no_interactions <- lmer(Frequency ~ Cluster + Phase + Condition + Pre_value + (1 | Participant), data = relative_complete_phases)

# Model without pre_value
model_no_prevalue <- lmer(Frequency ~ Cluster * Phase * Condition + (1 | Participant), data = relative_complete_phases)

# Model with only main effects
model_main_effects <- lmer(Frequency ~ Cluster + Phase + Condition + (1 | Participant), data = relative_complete_phases)

# Null model (no fixed effects)
model_null <- lmer(Frequency ~ 1 + (1 | Participant), data = relative_complete_phases)

anova(model_null, model_main_effects, model_no_prevalue, model_no_interactions, model_full)

AIC(model_null, model_main_effects, model_no_prevalue, model_no_interactions, model_full)

BIC(model_null, model_main_effects, model_no_prevalue, model_no_interactions, model_full)

anova_results <- anova(model_full)
anova_table <- as.data.frame(anova_results)
print(anova_table)

emm <- emmeans(model_full, ~ Cluster * Phase * Condition)

# Pairwise comparisons with adjustment for multiple testing (Tukey method)
pairwise_comparisons <- pairs(emm, adjust = "Tukey")

# Displaying the results in a clean format
summary(pairwise_comparisons)

# considering only cluster 1
cluster1 <- data_filter(relative_complete_phases, Cluster== "1")
cluster2 <- data_filter(relative_complete_phases, Cluster== "2")

# Full model with all fixed effects and interactions
model_full <- lmer(Frequency ~ Phase * Condition + Pre_value + (1 | Participant), data = cluster2)

# Model without interaction terms
model_no_interactions <- lmer(Frequency ~  Phase + Condition + Pre_value + (1 | Participant), data = cluster2)

# Model without pre_value
model_no_prevalue <- lmer(Frequency ~  Phase * Condition + (1 | Participant), data = cluster2)

# Model with only main effects
model_main_effects <- lmer(Frequency ~ Phase + Condition + (1 | Participant), data = cluster2)

# Null model (no fixed effects)
model_null <- lmer(Frequency ~ 1 + (1 | Participant), data = cluster2)

anova(model_null, model_main_effects, model_no_prevalue, model_no_interactions, model_full)