PFMF23

rm(list = ls())

2 Pre Processing

Loading the packages….

invisible({capture.output({
# Packages for stats 
library(lme4)
library(afex)
library(broom.mixed)
library(metan)
library(nlme)
library(stats)
library(purrr)
library(BayesFactor)
library(rmcorr) 


# Package to perform the Anova tablet contrast the results
library(car)
library(emmeans)
  library(easystats)
  library(dunn.test)
  library(coin)

#Package for clustering 
library(cluster)
library(factoextra)
library(stats)
library(clusterCrit)
  
#Package for load online data
library(osfr)    #To access to OSF
library("httr")  #To load url links

#For the figures and tables
library(Rmisc)
library(readxl)
library(PupillometryR)
library(smplot2)
library(ggeasy)
library(pheatmap)
library(reshape2)
library(kableExtra)
library(tidyverse)
library(writexl)
require(ggpubr)
library(dplyr)
library(plotly)
library(cowplot) 
library(RColorBrewer) 
library(ggrepel)
})})

3 Manipulation check

MVC

We assess if the intervention was effective to increase physical fatigue objectively by means of a MVC test before and after the running exercise or the walking control condition; and by means of the rate or perceived exertion at subjective level.

# Standardize size of the text across all figures. We will use this for all figures.
size = theme(
  title = element_text(size = 12),
  axis.text.x = element_text(size = 12),
  axis.title.y = element_text(size = 12),
  axis.title.x = element_text(size = 12))


invisible({capture.output({
  
  # Read data from OSF.
  
  url <- 'https://osf.io/r42gz//?action=download'
  filename <- "data_PFMF23.xlsx"
  GET(url, write_disk(filename, overwrite = TRUE))
  MVC <- readxl::read_xlsx(filename, sheet= "MVC")
  RPE <- readxl::read_xlsx(filename, sheet= "RPE")
  MT <- readxl::read_xlsx(filename, sheet= "MT")
  
})})


# Calculate ratios for the MVC
MVC <- MVC %>%
  mutate(
    FAT = (MVC_POS_FAT - MVC_PRE_FAT) / (MVC_POS_FAT + MVC_PRE_FAT))

MVC <- MVC %>%
  mutate(
    CON = (MVC_POS_CON - MVC_PRE_CON) / (MVC_POS_CON + MVC_PRE_CON))  

MVC <- MVC[, c(1,6:7)]
  
MVC_long <- pivot_longer(MVC, cols =2:3, names_to = c("Condition"),
                           values_to = "Value", values_drop_na = TRUE)


MVCRain <- ggplot(data = MVC_long, mapping = aes(x = Condition, y = Value)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Condition, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  geom_point(aes(fill = Condition, group=Participant), size=2,shape=21)+
  geom_line(aes(group=Participant),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Normalized score (AU)", x = element_blank(),
        title ="MVC change",
        caption = "n = 29")+
  scale_colour_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

MVCRain

bf_result_MVC <- ttestBF(x = MVC$FAT, y = MVC$CON, paired = TRUE)
bf_result_MVC
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 59.87906 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = MVC$FAT, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = MVC$CON, nullInterval = NULL)  # Against 0

# Extract posterior samples for each mean
posterior_exp <- posterior(bf_exp, iterations = 10000)
posterior_con <- posterior(bf_con, iterations = 10000)

# Convert to data frames
posterior_exp_df <- as.data.frame(posterior_exp)
posterior_con_df <- as.data.frame(posterior_con)

# Calculate the 95% credible intervals for each mean
credible_interval_exp <- quantile(posterior_exp_df$mu, probs = c(0.025, 0.975))
credible_interval_con <- quantile(posterior_con_df$mu, probs = c(0.025, 0.975))

summary_stats <- MVC_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value)
  )

summary_stats
## # A tibble: 2 × 3
##   Condition     mean     sd
##   <chr>        <dbl>  <dbl>
## 1 CON       -0.00290 0.0250
## 2 FAT       -0.0396  0.0476
# Print the credible intervals
print("Credible Interval for Fixation Exp:")
## [1] "Credible Interval for Fixation Exp:"
print(credible_interval_exp)
##        2.5%       97.5% 
## -0.05536238 -0.01898197
print("Credible Interval for Fixation Con:")
## [1] "Credible Interval for Fixation Con:"
print(credible_interval_con)
##         2.5%        97.5% 
## -0.011381120  0.006104872
t.test(x = MVC$FAT, y = MVC$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  MVC$FAT and MVC$CON
## t = -3.922, df = 28, p-value = 0.0005178
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.05584968 -0.01752620
## sample estimates:
## mean difference 
##     -0.03668794

RPE

RPE_long <- pivot_longer(RPE, cols =2:3, names_to = c("Condition"),
                         values_to = "Value", values_drop_na = TRUE)


RPERain <- ggplot(data = RPE_long, mapping = aes(x = Condition, y = Value)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Condition, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  scale_colour_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+
  geom_point(aes(fill = Condition, group=Participant), size=2,shape=21)+
  geom_line(aes(group=Participant),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( title = "RPE",
    y = "RPE 1-10", x = element_blank(),
        caption = "n = 29")+
  size

RPERain

bf_result_rpe <- ttestBF(x = RPE$FAT, y = RPE$CON, paired = TRUE)

bf_result_rpe
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 3.16384e+21 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
bf_exp <- ttestBF(x = RPE$FAT, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = RPE$CON, nullInterval = NULL)  # Against 0

# Extract posterior samples for each mean
posterior_exp <- posterior(bf_exp, iterations = 10000)
posterior_con <- posterior(bf_con, iterations = 10000)

# Convert to data frames
posterior_exp_df <- as.data.frame(posterior_exp)
posterior_con_df <- as.data.frame(posterior_con)

# Calculate the 95% credible intervals for each mean
credible_interval_exp <- quantile(posterior_exp_df$mu, probs = c(0.025, 0.975))
credible_interval_con <- quantile(posterior_con_df$mu, probs = c(0.025, 0.975))

# Print the credible intervals
print("Credible Interval for RPE EXP:")
## [1] "Credible Interval for RPE EXP:"
print(credible_interval_exp)
##     2.5%    97.5% 
## 7.802398 8.490682
print("Credible Interval for RPE CON:")
## [1] "Credible Interval for RPE CON:"
print(credible_interval_con)
##      2.5%     97.5% 
## 0.5100812 1.0641797
summary_stats <- RPE_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value)
  )

summary_stats
## # A tibble: 2 × 3
##   Condition  mean    sd
##   <chr>     <dbl> <dbl>
## 1 CON       0.828 0.723
## 2 FAT       8.16  0.877
t.test(x = RPE$FAT, y = RPE$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  RPE$FAT and RPE$CON
## t = 35.478, df = 28, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  6.904514 7.750658
## sample estimates:
## mean difference 
##        7.327586
ManipulationPlot <- ggarrange(MVCRain, RPERain,
                      labels = c("A", "B"),
                      ncol = 2, nrow = 1, 
                      # common.legend = T,
                      # legend = "right",
                      align = "none"
)

ManipulationPlot

4 TET

set.seed(123)
#Read data

invisible({capture.output({

  # Read data from OSF.
url <- 'https://osf.io/ud2zp//?action=download'
filename <- 'PMF23_TET_combined_data_FAT.csv'
GET(url, write_disk(filename, overwrite = TRUE))
TET_fat <- read.csv(filename)


url <- 'https://osf.io/u5xb6//?action=download'
filename <- 'PMF23_TET_combined_data_CON.csv'
GET(url, write_disk(filename, overwrite = TRUE))
TET_con <- read.csv(filename)

})})

TET_fat$Condition <- rep("Fat", 16793) #add condition based on the total number of observations
TET_fat$Dimension <- as.factor(TET_fat$Dimension)
TET_fat$Participant <- factor(TET_fat$Participant)
TET_fat$Condition <- factor(TET_fat$Condition)


TET_fat <- TET_fat %>%
  group_by(Participant) %>%
  mutate(Epoch = rep(1:(n() / length(unique(Dimension))), each = length(unique(Dimension)))) %>%
  ungroup()
TET_fat$Epoch <- factor(TET_fat$Epoch)

# transform in wide format for the clustering algorithm 
TET_fat <- TET_fat %>% 
  pivot_wider(
    names_from  = c(Dimension), 
    values_from = c(Value)
  )


# Control condition
TET_con$Condition <- rep("Con", 17143) #add condition based on the total number of observations
TET_con$Dimension <- as.factor(TET_con$Dimension)
TET_con$Participant <- factor(TET_con$Participant)
TET_con$Condition <- factor(TET_con$Condition)


# To add a order for each data point. It is needed for wide format
TET_con <- TET_con %>%
  group_by(Participant) %>%
  mutate(Epoch = rep(1:(n() / length(unique(Dimension))), each = length(unique(Dimension)))) %>%
  ungroup()
TET_con$Epoch <- factor(TET_con$Epoch)

# transform in wide format for the clustering algorithm 
TET_con <- TET_con %>% 
  pivot_wider(
    names_from  = c(Dimension), 
    values_from = c(Value)
  )

# Merge both experimental and control conditions
TET_merged <- rbind(TET_fat, TET_con)

names(TET_merged)[6] <- "Mind-wandering"

Finding the opptimal number of clusters. We will use the the Elbow, Silhouette methods and Caliński-Harabasz score.

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

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

# 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 = TET_merged[, 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")
}
## k = 2 : CH Score = 4437.253 
## k = 3 : CH Score = 3476.765 
## k = 4 : CH Score = 2751.299 
## k = 5 : CH Score = 2354.846 
## k = 6 : CH Score = 2028.48 
## k = 7 : CH Score = 1872.045 
## k = 8 : CH Score = 1757.752 
## k = 9 : CH Score = 1651.08 
## k = 10 : CH Score = 1541.467
# 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")

It seems that two clusters is the optimal number of clusters for this dataset

km.out <- kmeans(TET_merged[, 4:10], centers=2,nstart=100)
km.out$centers
##     Fatigue     Fedup Mind-wandering Execution Concentration   Boredom
## 1 0.4017923 0.7508808      0.2658710 0.7231078     0.2730512 0.1956047
## 2 0.7552468 0.3085921      0.6401075 0.3352308     0.6969019 0.6805452
##      Effort
## 1 0.2582507
## 2 0.6464405
TET_merged$Cluster <-km.out$cluster

The weighted average intensity or centroid is represented there, and it can be visualized here:

Cluster_long<- pivot_longer(TET_merged, cols = c(Fatigue, Fedup, 'Mind-wandering', Execution, Concentration, Boredom, Effort), names_to = "Dimension",
                        values_to = "Value", values_drop_na = TRUE)



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


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

# # Define your custom labels for the clusters
custom_labels <- c("1" = "Cluster 1 (low demand)", "2" = "Cluster 2 (high demand)")
custom_colors <- c("1" = "#999999", "2" = "#F781BF")  # Example Set2 colors

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_boxplot(aes(x = Dimension, y = Value, fill = Cluster),outlier.shape = NA, alpha = .5, width = .3, colour = "black")+
  geom_point(data = sumrepdatClust, aes(x = Dimension, y = Value, group = Cluster, colour = Cluster), shape = 18, size = 3) +
 scale_colour_manual(values = custom_colors, labels = custom_labels) +   # Custom labels and colors for color
  scale_fill_manual(values = custom_colors, labels = custom_labels) +
  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))




# 
# # Modify your ggplot code
# 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_boxplot(aes(x = Dimension, y = Value, fill = Cluster), outlier.shape = NA, alpha = .5, width = .3, colour = "black") +
#   geom_point(data = sumrepdatClust, aes(x = Dimension, y = Value, group = Cluster, colour = Cluster), shape = 18, size = 3) +
#   scale_colour_manual(values = custom_colors, labels = custom_labels) +   # Custom labels and colors for color
#   scale_fill_manual(values = custom_colors, labels = custom_labels) +     # Custom labels and colors for fill
#   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

Is there evidence of a difference between the two clusters?

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

data_relative_cluster <- data_relative_cluster[, c(1:6)]


# Define the dimensions to loop over
dimensions <- c("Effort", "Mind-wandering", "Boredom", "Concentration", "Fatigue", "Fedup", "Execution")

# Initialize a list to store Bayes Factor results
bayes_factors <- list()

for (dim in dimensions) {
    # Subset the data for the current dimension
    data_subset <- subset(data_relative_cluster, Dimension == dim)
    
    # Ensure that 'Cluster' is a factor (if not already)
    data_subset$Cluster <- as.factor(data_subset$Cluster)
    
    # Perform the Bayesian independent t-test
    bf_test <- ttestBF(formula = Value ~ Cluster, data = data_subset)
    
    # Extract the Bayes Factor
    bf_value <- as.numeric(extractBF(bf_test)$bf)
    
    # Store the result with the corresponding dimension
    bayes_factors[[dim]] <- data.frame(Dimension = dim, BayesFactor = bf_value)
}

# Combine all Bayes Factors into a single data frame
combined_bf_results <- do.call(rbind, bayes_factors)

# # Print the combined Bayes Factors
# print(combined_bf_results)

# Create a table for model results
kable(combined_bf_results, format = "html", table.attr = "class='table table-striped'") %>%
  kable_styling(full_width = F, position = "left", bootstrap_options = c("striped", "hover", "condensed"))
Dimension BayesFactor
Effort Effort 2.904331e+10
Mind-wandering Mind-wandering 1.557482e+14
Boredom Boredom 2.939886e+10
Concentration Concentration 2.544063e+12
Fatigue Fatigue 2.681505e+07
Fedup Fedup 3.070581e+14
Execution Execution 2.645532e+16

We can also see the average intensity further descriptions.

# 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"))
Cluster count mean sd se ci_lower ci_upper Dimension
1 48 0.2723227 0.1737210 0.0250745 0.2218794 0.3227660 Effort
2 51 0.6116423 0.2176544 0.0304777 0.5504260 0.6728586 Effort
1 48 0.2645126 0.1392841 0.0201039 0.2240687 0.3049565 Mind-wandering
2 51 0.6261187 0.2011230 0.0281629 0.5695519 0.6826854 Mind-wandering
1 48 0.2363964 0.1451311 0.0209479 0.1942547 0.2785381 Boredom
2 51 0.6113452 0.2694953 0.0377369 0.5355485 0.6871420 Boredom
1 48 0.3321259 0.2220498 0.0320501 0.2676493 0.3966024 Concentration
2 51 0.6945362 0.1541550 0.0215860 0.6511794 0.7378930 Concentration
1 48 0.4252060 0.2415324 0.0348622 0.3550723 0.4953397 Fatigue
2 51 0.7314042 0.1905112 0.0266769 0.6778221 0.7849864 Fatigue
1 48 0.7249202 0.1905229 0.0274996 0.6695982 0.7802423 Fedup
2 51 0.3285599 0.1857294 0.0260073 0.2763227 0.3807971 Fedup
1 48 0.7167098 0.1564772 0.0225855 0.6712736 0.7621460 Execution
2 51 0.3382278 0.1724593 0.0241491 0.2897229 0.3867328 Execution

Cluster representation depends on the condition?

# Obtain the number of data points within each cluster for all participants at every epoch
stacked_data <- summarySE(Cluster_long, measurevar = "Value",
                       groupvars=c("Epoch", "Condition", "Cluster"))


# Filter the conditions for the stacked plots
stacked_exp <- stacked_data %>% filter(Condition == "Fat")
  stacked_exp <- stacked_exp[,c(1,3:4)]

  stacked_con <- stacked_data %>% filter(Condition == "Con")
  stacked_con <- stacked_con[,c(1,3:4)]
  

  # obtain relative frequency of each cluster at each epoch
exp_relative <- stacked_exp %>%
  group_by(Epoch) %>%
  mutate(Frequency = N / sum(N)) %>%
  ungroup()


exp_relative$Epoch <- as.numeric(exp_relative$Epoch)

# To add vertical lines for the duration of each participant completing the task

num_bins <- TET_merged %>%
  group_by(Participant, Condition) %>%  # Group by Participant and Condition
  summarise(Max_Epoch = max(as.numeric(as.character(Epoch)))) %>%  # Convert Epoch to numeric and get max
  ungroup()  # Ungroup to remove grouping structure


num_bins_exp <-num_bins %>% filter(Condition == "Fat")
num_bins_con <-num_bins %>% filter(Condition == "Con")

# # Define your custom labels for the clusters
custom_labels <- c("1" = "Cluster 1 (low demand)", "2" = "Cluster 2 (high demand)")
custom_colors <- c("1" = "#999999", "2" = "#F781BF")  # Example Set2 colors


plot_stacked_exp <- ggplot(exp_relative, aes(x = as.numeric(as.character(Epoch)), y = Frequency, fill = Cluster)) +
  geom_area(position = 'fill') +
   sm_classic(legends = FALSE)+
  easy_remove_x_axis(what = "text")+
  labs(title = "Experimental",
       x = "Task execution",
       y = "Relative Frequency of the cluster") +
scale_colour_manual(values = custom_colors, labels = custom_labels) +   # Custom labels and colors for color
  scale_fill_manual(values = custom_colors, labels = custom_labels) +    # Custom labels and colors for fill
   # Add vertical lines for each Participant's Max_Epoch
  geom_vline(data = num_bins_exp, aes(xintercept = Max_Epoch, color = Participant), linetype = "dashed", size = 1)+
  geom_label_repel(data = num_bins_exp, aes(x = Max_Epoch, y = 1, label = Participant), 
                   angle = 90, size = 3, color = "black", inherit.aes = FALSE,
                   nudge_x = 0.2, nudge_y = 0.1, # Nudging to avoid overlap
                  direction = "y", # Keep the labels aligned in the y-direction
                  max.overlaps = Inf) 
  

plot_stacked_exp

con_relative <- stacked_con %>%
  group_by(Epoch) %>%
  mutate(Frequency = N / sum(N)) %>%
  ungroup()

con_relative$Epoch <- as.numeric(con_relative$Epoch)

plot_stacked_con <- ggplot(con_relative, aes(x = Epoch, y = Frequency, fill = Cluster)) +
  geom_area(position = 'stack') +
   sm_classic(legends = FALSE)+
  easy_remove_x_axis(what = "text")+
  labs(title = "Control",
       x = "Task duration",
       y = "Relative Frequency of the cluster") +
  theme(axis.title.y=element_blank(),
        axis.text.x=element_blank(),
        )+
scale_colour_manual(values = custom_colors, labels = custom_labels) +   # Custom labels and colors for color
  scale_fill_manual(values = custom_colors, labels = custom_labels)+    # Custom labels and colors for fill
   geom_vline(data = num_bins_exp, aes(xintercept = Max_Epoch, color = Participant), linetype = "dashed", size = 1)+
  geom_label_repel(data = num_bins_con, aes(x = Max_Epoch, y = 1, label = Participant), 
                   angle = 90, size = 3, color = "black", inherit.aes = FALSE,
                   nudge_x = 0.2, nudge_y = 0.1, # Nudging to avoid overlap
                  direction = "y", # Keep the labels aligned in the y-direction
                  max.overlaps = Inf) 


plot_stacked_con

For the article, both together

Stacked_plots <- ggarrange(plot_stacked_exp, plot_stacked_con,
                      ncol =2 ,labels = c("A", "B"),
                      nrow=1,
                      common.legend = T,
                      legend = "right",
                      align = "h")

Stacked_plots

In any case, at a glance, one cluster is more likely to occur than another within each condition. We can calculate the odds ratios for this.

# Create a contingency table

contingency_table <- table(Cluster_long$Condition, Cluster_long$Cluster)
print(contingency_table)
##      
##          1    2
##   Fat 9121 7672
##   Con 9408 7735
# Extract counts

exp_1 <- contingency_table["Fat", "1"]
exp_2 <- contingency_table["Fat", "2"]
con_1 <- contingency_table["Con", "1"]
con_2 <- contingency_table["Con", "2"]

# Calculate the odds
odds_exp <- exp_2 / exp_1
odds_con <- con_2 / con_1

# Calculate the Odds Ratio
odds_ratio <- odds_exp / odds_con

cat("Odds for experimental condition (Cluster 2/Cluster 1):", odds_exp, "\n")
## Odds for experimental condition (Cluster 2/Cluster 1): 0.8411358
cat("Odds for control condition (Cluster 2/Cluster 1):", odds_con, "\n")
## Odds for control condition (Cluster 2/Cluster 1): 0.8221726
# cat("Odds Ratio (Control/Experimental):", odds_ratio, "\n")

# Fisher's Exact Test
fisher_test <- fisher.test(contingency_table)
print(fisher_test)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  contingency_table
## p-value = 0.3003
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.9363449 1.0204019
## sample estimates:
## odds ratio 
##  0.9774485
# Chi-Square Test
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_table
## X-squared = 1.0707, df = 1, p-value = 0.3008

5 Lempel Ziv 2

invisible({capture.output({

url <- 'https://osf.io/pmq73//?action=download'
filename <- "combined_LZ2.xlsx"
GET(url, write_disk(filename, overwrite = TRUE))
LZ_c <- readxl::read_xlsx(filename, sheet= "Feuil1")

url <- 'https://osf.io/pmq73//?action=download'
filename <- "combined_LZ2.xlsx"
GET(url, write_disk(filename, overwrite = TRUE))
LZ_sum <- readxl::read_xlsx(filename, sheet= "Feuil2")


})})

LZ_c <- LZ_c %>%
  separate(Participant, into = c("Participant", "Condition"), sep = "_")

LZ_sum <- LZ_sum %>%
  separate(Participant, into = c("Participant", "Condition"), sep = "_")

LZc

LZ_c <- LZ_c %>% mutate(across(3:3660, as.numeric))
LZ_c$Participant <- as.numeric(LZ_c$Participant)
LZ_c$Participant <- paste0("P", LZ_c$Participant)


LZ_c_long <- LZ_c %>% 
 pivot_longer(
    cols = 3:3660,         # Select columns 3 to 3656 (these are the time point columns)
    names_to = "Time",      # Create a 'Time' column to store the time points (column names)
    values_to = "Value"     # Create a 'Value' column to store the actual values
  ) %>%
  mutate(Time = as.numeric(Time))

LZ_c_long <- LZ_c_long %>% drop_na(Value)


LZ_c_long <- LZ_c_long %>%
  left_join(num_bins, by = c("Participant", "Condition"))


# Function to calculate binned averages of LZc
calculate_bin_averages <- function(LZc_data, num_bins) {
  # Calculate the bin index for each data point
  bin_indices <- cut(seq_along(LZc_data), breaks = num_bins, labels = FALSE)
  
  # Calculate the average for each bin
  bin_means <- tapply(LZc_data, bin_indices, mean)
  
  return(bin_means)
}

#To replace the Na in Max epoch of those participants which we did not process the TET
LZ_c_long <- LZ_c_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P4" & Condition == "Con" & is.na(Max_Epoch), 170, Max_Epoch)
  )

LZ_c_long <- LZ_c_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P6" & Condition == "Con" & is.na(Max_Epoch), 44, Max_Epoch)
  )

LZ_c_long <- LZ_c_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P8" & Condition == "Con" & is.na(Max_Epoch), 149, Max_Epoch)
  )

LZ_c_long <- LZ_c_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P10" & Condition == "Fat" & is.na(Max_Epoch), 112, Max_Epoch)
  )



binned_LzC <- LZ_c_long %>%
  group_by(Participant, Condition) %>%
  summarise(
    Epoch = 1:unique(Max_Epoch),  # Create bin index
    LZc = calculate_bin_averages(Value, unique(Max_Epoch))  # Calculate the bin averages
  ) %>%
  ungroup()
# Calculate average value across time to include in the figure.
average_data_fat <- binned_LzC %>%
  group_by(Epoch, Condition) %>%
  filter(Condition =="Fat") %>% 
  summarise(average_change = mean(LZc, na.rm = TRUE))

average_data_con <- binned_LzC %>%
  group_by(Epoch, Condition) %>%
  filter(Condition =="Con") %>% 
  summarise(average_change = mean(LZc, na.rm = TRUE))


LzC_TOT <- ggplot(binned_LzC, aes(x = Epoch, y = LZc), fill = Condition) +
    geom_line(aes(color = Participant), alpha = 0.5) + 
    geom_line(data = average_data_con, aes(x = Epoch, y = average_change), color = "red", size = 1.5) +
    geom_line(data = average_data_fat, aes(x = Epoch, y = average_change), color = "blue", size = 1.5) +
  labs(title = "Time course of Lempel Ziv(c) across the task",
       x = "Epoch",
       y = "Arbinatry units ") +
  sm_minimal()


LzC_TOT

Calculate average complexity in both conditions

average_LZc <- binned_LzC %>%
  group_by(Participant, Condition) %>%
  summarise(average_change = mean(LZc, na.rm = TRUE))

LZcPlot <- ggplot(data = average_LZc, mapping = aes(x = Condition, y = average_change)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Condition, y = average_change, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  geom_point(aes(fill = Condition, group=Participant), size=2,shape=21)+
  geom_line(aes(group=Participant),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Lempel Ziv complenxity (AU)", x = element_blank(),
        title ="Average LZc complexity")+
  scale_colour_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

LZcPlot 

average_LZc_wide <- pivot_wider(average_LZc,  names_from = "Condition",
    values_from = "average_change")

average_LZc_wide <- drop_na(average_LZc_wide)

Stats:

ttestBF(x = average_LZc_wide$Con, y = average_LZc_wide$Fat, paired = TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 7.874107 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
t.test(x = average_LZc_wide$Con, y = average_LZc_wide$Fat, paired = TRUE)
## 
##  Paired t-test
## 
## data:  average_LZc_wide$Con and average_LZc_wide$Fat
## t = 3.0308, df = 27, p-value = 0.005327
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.01258012 0.06531239
## sample estimates:
## mean difference 
##      0.03894626
bf_exp <- ttestBF(x = average_LZc_wide$Fat, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = average_LZc_wide$Con, nullInterval = NULL)  # Against 0

# Extract posterior samples for each mean
posterior_exp <- posterior(bf_exp, iterations = 10000)
posterior_con <- posterior(bf_con, iterations = 10000)

# Convert to data frames
posterior_exp_df <- as.data.frame(posterior_exp)
posterior_con_df <- as.data.frame(posterior_con)

# Calculate the 95% credible intervals for each mean
credible_interval_exp <- quantile(posterior_exp_df$mu, probs = c(0.025, 0.975))
credible_interval_con <- quantile(posterior_con_df$mu, probs = c(0.025, 0.975))


# Print the credible intervals
print("Credible Interval for LZc EXP:")
## [1] "Credible Interval for LZc EXP:"
print(credible_interval_exp)
##      2.5%     97.5% 
## 0.6957501 0.7402365
print("Credible Interval for LZc CON:")
## [1] "Credible Interval for LZc CON:"
print(credible_interval_con)
##      2.5%     97.5% 
## 0.7392714 0.7744453
summary_stats <- average_LZc %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(average_change)
  )

summary_stats
## # A tibble: 2 × 2
##   Condition  mean
##   <chr>     <dbl>
## 1 Con       0.757
## 2 Fat       0.715

LZ Sum

LZ_sum <- LZ_sum %>% mutate(across(3:3660, as.numeric))
LZ_sum$Participant <- as.numeric(LZ_sum$Participant)
LZ_sum$Participant <- paste0("P", LZ_sum$Participant)


LZ_sum_long <- LZ_sum %>% 
 pivot_longer(
    cols = 3:3660,         # Select columns 3 to 3656 (these are the time point columns)
    names_to = "Time",      # Create a 'Time' column to store the time points (column names)
    values_to = "Value"     # Create a 'Value' column to store the actual values
  ) %>%
  mutate(Time = as.numeric(Time))

LZ_sum_long <- LZ_sum_long %>% drop_na(Value)


LZ_sum_long <- LZ_sum_long %>%
  left_join(num_bins, by = c("Participant", "Condition"))


# Function to calculate binned averages of LZc
calculate_bin_averages <- function(LZc_data, num_bins) {
  # Calculate the bin index for each data point
  bin_indices <- cut(seq_along(LZc_data), breaks = num_bins, labels = FALSE)
  
  # Calculate the average for each bin
  bin_means <- tapply(LZc_data, bin_indices, mean)
  
  return(bin_means)
}

#To replace the Na in Max epoch of those participants which we did not process the TET
LZ_sum_long <- LZ_sum_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P4" & Condition == "Con" & is.na(Max_Epoch), 170, Max_Epoch)
  )

LZ_sum_long <- LZ_sum_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P6" & Condition == "Con" & is.na(Max_Epoch), 44, Max_Epoch)
  )

LZ_sum_long <- LZ_sum_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P8" & Condition == "Con" & is.na(Max_Epoch), 149, Max_Epoch)
  )


LZ_sum_long <- LZ_sum_long %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P10" & Condition == "Fat" & is.na(Max_Epoch), 112, Max_Epoch)
  )



binned_LzSum <- LZ_sum_long %>%
  group_by(Participant, Condition) %>%
  summarise(
    Epoch = 1:unique(Max_Epoch),  # Create bin index
    LZc = calculate_bin_averages(Value, unique(Max_Epoch))  # Calculate the bin averages
  ) %>%
  ungroup()
# Calculate average value across time to include in the figure.
average_data_fat <- binned_LzSum %>%
  group_by(Epoch, Condition) %>%
  filter(Condition =="Fat") %>% 
  summarise(average_change = mean(LZc, na.rm = TRUE))

average_data_con <- binned_LzSum %>%
  group_by(Epoch, Condition) %>%
  filter(Condition =="Con") %>% 
  summarise(average_change = mean(LZc, na.rm = TRUE))


LzSum_TOT <- ggplot(binned_LzSum, aes(x = Epoch, y = LZc), fill = Condition) +
    geom_line(aes(color = Participant), alpha = 0.5) + 
    geom_line(data = average_data_con, aes(x = Epoch, y = average_change), color = "red", size = 1.5) +
    geom_line(data = average_data_fat, aes(x = Epoch, y = average_change), color = "blue", size = 1.5) +
  labs(title = "Time course of Lempel Ziv(Sum) across the task",
       x = "Epoch",
       y = "Arbinatry units ") +
  sm_minimal()

LzSum_TOT

Calculate average complexity in both conditions

average_LZSum <- binned_LzSum %>%
  group_by(Participant, Condition) %>%
  summarise(average_change = mean(LZc, na.rm = TRUE))

LZsumPlot <- ggplot(data = average_LZSum, mapping = aes(x = Condition, y = average_change)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Condition, y = average_change, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  geom_point(aes(fill = Condition, group=Participant), size=2,shape=21)+
  geom_line(aes(group=Participant),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Lempel Ziv complenxity (AU)", x = element_blank(),
        title ="Average LZSum complexity")+
  scale_colour_brewer(palette = "Set1")+
  scale_fill_brewer(palette = "Set1")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

LZsumPlot 

average_LZsum_wide <- pivot_wider(average_LZSum,  names_from = "Condition",
    values_from = "average_change")

average_LZsum_wide <- drop_na(average_LZsum_wide)

ttestBF(x = average_LZsum_wide$Con, y = average_LZsum_wide$Fat, paired = TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 1452340 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
t.test(x = average_LZsum_wide$Con, y = average_LZsum_wide$Fat, paired = TRUE)
## 
##  Paired t-test
## 
## data:  average_LZsum_wide$Con and average_LZsum_wide$Fat
## t = 8.1745, df = 27, p-value = 8.858e-09
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.05353630 0.08941852
## sample estimates:
## mean difference 
##      0.07147741

stats:

bf_exp <- ttestBF(x = average_LZsum_wide$Fat, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = average_LZsum_wide$Con, nullInterval = NULL)  # Against 0

# Extract posterior samples for each mean
posterior_exp <- posterior(bf_exp, iterations = 10000)
posterior_con <- posterior(bf_con, iterations = 10000)

# Convert to data frames
posterior_exp_df <- as.data.frame(posterior_exp)
posterior_con_df <- as.data.frame(posterior_con)

# Calculate the 95% credible intervals for each mean
credible_interval_exp <- quantile(posterior_exp_df$mu, probs = c(0.025, 0.975))
credible_interval_con <- quantile(posterior_con_df$mu, probs = c(0.025, 0.975))


# Print the credible intervals
print("Credible Interval for LZ sum EXP:")
## [1] "Credible Interval for LZ sum EXP:"
print(credible_interval_exp)
##      2.5%     97.5% 
## 0.4766975 0.5046144
print("Credible Interval for LZ sum CON:")
## [1] "Credible Interval for LZ sum CON:"
print(credible_interval_con)
##      2.5%     97.5% 
## 0.5460380 0.5781498
summary_stats <- average_LZSum %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(average_change)
  )

summary_stats
## # A tibble: 2 × 2
##   Condition  mean
##   <chr>     <dbl>
## 1 Con       0.562
## 2 Fat       0.491

Combine plots Lempel:

LempelPlots <- ggarrange(LZcPlot, LzC_TOT, LZsumPlot,LzSum_TOT,
                      ncol =2 ,labels = c("A", "", "B"),
                      nrow=2,
                      # common.legend = T,
                      # legend = "right",
                      align = "h")


LempelPlots

Cognitive task until failure

invisible({capture.output({
  
  # Read data from OSF.
  
  url <- 'https://osf.io/r42gz//?action=download'
  filename <- "data_PFMF23.xlsx"
  GET(url, write_disk(filename, overwrite = TRUE))
  TaskTTE <- readxl::read_xlsx(filename, sheet= "Task")
  
  url <- 'https://osf.io/q94zd//?action=download'
  filename <- "combined_data_task_PFMF.txt"
  GET(url, write_disk(filename, overwrite = TRUE))
  Task_perf <- read.delim(filename)
  
  
  
})})


# Task_perf <- read.delim("~/Library/CloudStorage/OneDrive-UniversitédeLausanne/PFMF/Codes/combined_data_task_PFMF.txt")

Task_perf <-  Task_perf %>% 
separate(Subject, into = c("Participant", "Condition"), sep = "_")
Task_perf <- Task_perf[, c(1,2,6)]

Task_perf$Participant <- as.numeric(Task_perf$Participant)
Task_perf$Participant <- paste0("P", Task_perf$Participant)
Task_perf$Condition <- as.factor(Task_perf$Condition)


Task_perf <- Task_perf %>%
  left_join(num_bins, by = c("Participant", "Condition"))



# Function to calculate binned averages of Task
calculate_bin_averages <- function(Task_data, num_bins) {
  # Calculate the bin index for each data point
  bin_indices <- cut(seq_along(Task_data), breaks = num_bins, labels = FALSE)
  
  # Calculate the average for each bin
  bin_means <- tapply(Task_data, bin_indices, mean)
  
  return(bin_means)
}



#To replace the Na in Max epoch of those participants which we did not process the TET
Task_perf <- Task_perf %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P4" & Condition == "Con" & is.na(Max_Epoch), 170, Max_Epoch)
  )

Task_perf <- Task_perf %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P6" & Condition == "Con" & is.na(Max_Epoch), 44, Max_Epoch)
  )

Task_perf <- Task_perf %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P8" & Condition == "Con" & is.na(Max_Epoch), 149, Max_Epoch)
  )

Task_perf <- Task_perf %>%
  mutate(
    Max_Epoch = ifelse(Participant == "P10" & Condition == "Fat" & is.na(Max_Epoch), 112, Max_Epoch)
  )


binned_Task_perf <- Task_perf %>%
  group_by(Participant, Condition) %>%
  summarise(
    Epoch = 1:unique(Max_Epoch),  # Create bin index
    Perf = calculate_bin_averages(PERFORMANCE, unique(Max_Epoch))  # Calculate the bin averages
  ) %>%
  ungroup()

The total time to complete the task in both conditions:

bf_result_task <- ttestBF(x = TaskTTE$Fat, y = TaskTTE$Con, paired = TRUE)

bf_result_task
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.3296233 ±0.03%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = TaskTTE$Fat, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = TaskTTE$Con, nullInterval = NULL)  # Against 0

# Extract posterior samples for each mean
posterior_exp <- posterior(bf_exp, iterations = 10000)
posterior_con <- posterior(bf_con, iterations = 10000)

# Convert to data frames
posterior_exp_df <- as.data.frame(posterior_exp)
posterior_con_df <- as.data.frame(posterior_con)

# Calculate the 95% credible intervals for each mean
credible_interval_exp <- quantile(posterior_exp_df$mu, probs = c(0.025, 0.975))
credible_interval_con <- quantile(posterior_con_df$mu, probs = c(0.025, 0.975))

# Print the credible intervals
print("Credible Interval for Perfomance EXP:")
## [1] "Credible Interval for Perfomance EXP:"
print(credible_interval_exp)
##     2.5%    97.5% 
## 3328.883 5109.181
print("Credible Interval for Perfomance CON:")
## [1] "Credible Interval for Perfomance CON:"
print(credible_interval_con)
##     2.5%    97.5% 
## 3895.648 5478.124
TaskTTELong <- pivot_longer(TaskTTE, cols = 2:3, names_to = "Condition")
TaskTTELong$Condition <- factor(TaskTTELong$Condition)
TaskTTELong <- TaskTTELong %>% mutate(Condition = recode(Condition, Con = "Control", Fat = "Fatigue"))



ggplot(TaskTTELong, aes(x = factor(Participant), y = value, fill = Condition)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7) +  # Adjust dodge and bar width
  scale_fill_brewer(palette = "Set1") +  # Set color palette for conditions
  theme_minimal() +
  labs(title = "Time on Task for Each Participant",
       x = "Participant",
       y = "Time (s)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(TaskTTELong, aes(x = factor(Participant), y = value, fill = Condition)) +
  geom_bar(stat = "identity", position = "dodge") +  # Bars are still side-by-side, but facets split by condition
  scale_fill_brewer(palette = "Set1") +  # Set color palette for conditions
  theme_minimal() +
  labs(title = "Time on Task for Each Participant (Faceted by Condition)",
       x = "Participant",
       y = "Time (s)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~ Condition)

TET and cognitive task

TET_perf <- merge(TET_merged, binned_Task_perf, by = c("Participant","Condition" , "Epoch"))

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

# To check the colors for the plots
# display.brewer.pal(n = 8, name = 'Set1')
# 
# brewer.pal(n = 8, name = 'Set1')


custom_x_labels <- c("1" = "Cluster 1 (low demand)", "2" = "Cluster 2 (high demand)")
custom_colors <- c("Fat" = "#377EB8", "Con" = "#E41A1C")  # Replace with the colors you prefercustom_labels <- c("Con" = "Control", "Fat" = "Fatigue")


#Divide by Cluster
PerfTET <- ggplot(data = TET_perf, mapping = aes(x = Cluster, y = Perf, fill = Condition)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Cluster, y = Perf, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  labs(y = "Performance (%)", x = element_blank(),
         fill = "Condition")+
  scale_fill_manual(values = custom_colors)+
  theme(plot.background = element_rect(colour = NA)) +
  scale_x_discrete(labels = custom_x_labels) + 
  theme(legend.position = "right") +
# Custom x-axis labels
  size

PerfTET

#Divide by Condition
PerfTET <- ggplot(data = TET_perf, mapping = aes(x = Condition, y = Perf, 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_boxplot(aes(x = Condition, y = Perf, fill = Cluster),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  guides(fill = "none") +
  labs(y = "Performance (%)", x = element_blank(),
        )+
  scale_colour_brewer(palette = "Set1", labels = custom_labels)+
  scale_fill_brewer(palette = "Set1", labels = custom_labels)+ 
  theme(plot.background = element_rect(colour = NA)) +
  scale_x_discrete(labels = custom_x_labels) +   # Custom x-axis labels
  size

PerfTET

Is performance different according to the cluster? We need to obtain the performance associate to each cluster and condition for each participant. Note that for some participants, they might have not contributed to both clusters, so the structure is unbalanced.

stats_perf_cluster <- TET_perf %>%
  group_by(Participant, Cluster, Condition) %>%
  summarise(
    Perf = mean(Perf)
  )

# stats_perf_cluster


# write.csv(stats_perf_cluster, file = "stats_perf_cluster.csv", row.names = FALSE)


# model <- lmer( Perf ~ Cluster * Condition + (1 | Participant), data = stats_perf_cluster)


bf_result <- anovaBF(Perf ~ Cluster + Condition + Participant, 
                     data = stats_perf_cluster,
                     whichRandom = "Participant")  # Random factor


bf_values <- c(8.44e+14, 0.2260224, 1.72437e+14, 5.7809e+13)

# Compute Bayes factors relative to the best model (Model 1)
relative_bf <- bf_values / max(bf_values)

# Print the relative Bayes factors
relative_bf
## [1] 1.000000e+00 2.677991e-16 2.043092e-01 6.849408e-02
# Fit the model with Cluster and Participant to obtain the "main effect"

bf_model <- anovaBF(Perf ~ Cluster + Participant, 
                    data = stats_perf_cluster, 
                    whichRandom = "Participant")
posterior_samples <- posterior(bf_model, iterations = 10000)


# # Inspect the posterior samples
# summary(posterior_samples)
# 
# # Extract the cluster effect from the posterior samples
# cluster_effect <- posterior_samples[, "Cluster"]
# 
# # Compute the posterior mean and credible intervals for the cluster effect
# mean_cluster_effect <- mean(cluster_effect)
# credible_interval <- quantile(cluster_effect, probs = c(0.025, 0.975))
# 
# # Print the results
# mean_cluster_effect
# credible_interval

TET and complexiity

For LZc:

TET_lzc <- merge(TET_merged, binned_LzC, by = c("Participant","Condition" , "Epoch"))

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

TET_lzc$Conditon <- factor(TET_lzc$Condition, levels=c("Fat" , "Con"))

# To check the colors for the plots
# display.brewer.pal(n = 8, name = 'Set1')
# 
# brewer.pal(n = 8, name = 'Set1')

custom_colors <- c("Fat" = "#377EB8", "Con" = "#E41A1C")  # Replace with the colors you prefer


#Divide by Cluster
LzcTET <- ggplot(data = TET_lzc, mapping = aes(x = Cluster, y = LZc, fill = Condition)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Cluster, y = LZc, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  labs(y = "Complexity (AU)", x = element_blank(),
         fill = "Condition")+
  scale_fill_manual(values = custom_colors)+
  theme(plot.background = element_rect(colour = NA)) +
  scale_x_discrete(labels = custom_x_labels) +   # Custom x-axis labels
    # theme(legend.position = "right") +

  size

LzcTET

stats_lzc_cluster <- TET_lzc %>%
  group_by(Participant, Cluster, Condition) %>%
  summarise(
    LZc = mean(LZc)
  )

# stats_lzc_cluster

# write.csv(stats_lzc_cluster, file = "stats_lzc_cluster.csv", row.names = FALSE)

bf_result <- anovaBF(LZc ~ Cluster + Condition + Participant, 
                     data = stats_lzc_cluster,
                     whichRandom = "Participant")  # Random factor

bf_result
## Bayes factor analysis
## --------------
## [1] Cluster + Participant                                 : 0.2183384 ±1.27%
## [2] Condition + Participant                               : 361.6438  ±1.51%
## [3] Cluster + Condition + Participant                     : 82.43831  ±1.57%
## [4] Cluster + Condition + Cluster:Condition + Participant : 22.91265  ±1.91%
## 
## Against denominator:
##   LZc ~ Participant 
## ---
## Bayes factor type: BFlinearModel, JZS
bf_values <- c(0.2211868, 360.4836, 80.70137, 25.15605)

# Compute Bayes factors relative to the best model (Model 1)
relative_bf <- bf_values / max(bf_values)

# Print the relative Bayes factors
relative_bf
## [1] 0.0006135835 1.0000000000 0.2238697405 0.0697841733
bf_model <- anovaBF(LZc ~ Condition + Participant, 
                    data = stats_lzc_cluster, 
                    whichRandom = "Participant")

bf_model
## Bayes factor analysis
## --------------
## [1] Condition + Participant : 352.271 ±0.83%
## 
## Against denominator:
##   LZc ~ Participant 
## ---
## Bayes factor type: BFlinearModel, JZS

For Lempel Ziv Sum:

TET_lzsum <- merge(TET_merged, binned_LzSum, by = c("Participant","Condition" , "Epoch"))

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

TET_lzsum$Conditon <- factor(TET_lzsum$Condition, levels=c("Fat" , "Con"))

# To check the colors for the plots
# display.brewer.pal(n = 8, name = 'Set1')
# 
# brewer.pal(n = 8, name = 'Set1')

custom_colors <- c("Fat" = "#377EB8", "Con" = "#E41A1C")  # Replace with the colors you prefer


#Divide by Cluster
LzsumTET <- ggplot(data = TET_lzsum, mapping = aes(x = Cluster, y = LZc, fill = Condition)) +
  geom_flat_violin(aes(fill = Condition),position = position_nudge(x = .1, y = 0), adjust = 1.5, trim = FALSE, alpha = .5, colour = NA)+
  geom_boxplot(aes(x = Cluster, y = LZc, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  labs(y = "Complexity (AU)", x = element_blank(),
         fill = "Condition")+
  scale_fill_manual(values = custom_colors)+
  theme(plot.background = element_rect(colour = NA)) +
  scale_x_discrete(labels = custom_x_labels) +   # Custom x-axis labels
    theme(legend.position = "right") +

  size

LzsumTET

stats_lzsum_cluster <- TET_lzsum %>%
  group_by(Participant, Cluster, Condition) %>%
  summarise(
    LZc = mean(LZc)
  )

# stats_lzsum_cluster

write.csv(stats_lzsum_cluster, file = "stats_lzsum_cluster.csv", row.names = FALSE)

bf_result <- anovaBF(LZc ~ Cluster + Condition + Participant, 
                     data = stats_lzsum_cluster,
                     whichRandom = "Participant")  # Random factor

bf_result
## Bayes factor analysis
## --------------
## [1] Cluster + Participant                                 : 0.2357635    ±1.17%
## [2] Condition + Participant                               : 1.637835e+16 ±1.07%
## [3] Cluster + Condition + Participant                     : 4.809103e+15 ±11.1%
## [4] Cluster + Condition + Cluster:Condition + Participant : 1.224882e+15 ±2.15%
## 
## Against denominator:
##   LZc ~ Participant 
## ---
## Bayes factor type: BFlinearModel, JZS
bf_values <- c(0.2372368, 1.636497e+16, 4.122328e+15, 1.188756e+15)

# Compute Bayes factors relative to the best model (Model 1)
relative_bf <- bf_values / max(bf_values)

# Print the relative Bayes factors
relative_bf
## [1] 1.449662e-17 1.000000e+00 2.518995e-01 7.264028e-02
bf_model <- anovaBF(LZc ~ Condition + Participant, 
                    data = stats_lzc_cluster, 
                    whichRandom = "Participant")

bf_model
## Bayes factor analysis
## --------------
## [1] Condition + Participant : 354.5148 ±1.15%
## 
## Against denominator:
##   LZc ~ Participant 
## ---
## Bayes factor type: BFlinearModel, JZS
TETLZPlot <- ggarrange(LzcTET, LzsumTET,
                      labels = c("A", "B"),
                      ncol = 2, nrow = 1, 
                      # common.legend = T,
                      # legend = "right",
                      align = "none"
)

TETLZPlot

Exploratory analysis

# Correlation between the Mental toughness questionnaire and performance in the running exercise
correlationBF(x=MT$MT, y= MT$TTF,
              rscale = "medium",
              nullInterval = NULL,
              posterior = FALSE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.333 : 0.4068249 ±0%
## 
## Against denominator:
##   Null, rho = 0 
## ---
## Bayes factor type: BFcorrelation, Jeffreys-beta*