IMF23

rm(list = ls())

1. Intro

The main goal of this study was to shed light on the impact of performing a mental exertion task on a subsequent physical exercise by a fully individualized protocol. The main hypotheses of this study were:

  1. Participants would reach mental failure at different duration and multidimensional subjective sensations would be involved at varying rates.

  2. The change in subjective sensations and the physiological impact would cause subsequent physical performance to worsen.

  3. More exploratory, psychophysiological data would be closely related to the dynamics of task performance and subjective experience.

The corresponding material to reproduce the document is available in OSF, but this document is self-sufficient.

A brief reminder of the study protocol: Protocol

In the experimental condition, participant completed the task until the performance was below 50% in the last consecutive 60 trails. Meanwhile, in the control condition, participant watched a documentary/serie for 70 minutes. 70 minutes were based on the average duration of a previous study completing the cognitive task until failure.

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)
  
#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) 
# #Packages for logistics
#   library(this.path) #To set working directory where the code file is located
# For HRV
library(RHRV)


  
  
# #To set working directory where the code is located. 
#   setwd(this.path::here())
  
})})

…and data

invisible({capture.output({
  
  # Read data from OSF.

url <- 'https://osf.io/zr5x9//?action=download'
filename <- "data_IMF23.xlsx"
GET(url, write_disk(filename, overwrite = TRUE))
perf <- readxl::read_xlsx(filename, sheet= "Perf")
VAS <- readxl::read_xlsx(filename, sheet= "VAS")  
RPE <- readxl::read_xlsx(filename, sheet= "RPE")
Task <- readxl::read_xlsx(filename, sheet= "TaskD")

})})

# 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))

Task duration

Given that not everyone has the same cognitive resources and we do not know which is the optimal duration to tax on these resources, we proposed this novel approach until task failure. Note that a reduction in performance is one of the classical indices of fatigue in the literature. We expected that whit this individualized cognitive task, participants would reach failure at individualized time, being “sure” that everyone is mentally “exhausted” or willing to continue with the task.

Task stopped when the performance over the last 60 trials was below 50% or when participants reached the maximal time allowed in this experiment. Also they had to complete a minimum number of trails (~ 400/6 min)) before this moving average started to count performance

We can see the variety of task performance here:

ggplot(Task, aes(x = factor(Participant), y = Task)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  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))

summary_stats <- Task %>%
  summarise(
    mean = mean(Task),
    sd = sd(Task),
    range(Task)
  )
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Manipulation Check

We have several manipulation check for this study, either the subjective (VAS and TET) and the psychophysiological data (Eye-tracking and HRV).

VAS questions

We asked participants the following questions before and after the task and documentary:

  1. “What is your perception of mental fatigue now?”
  2. “What is your perception of boredom now?”

We need calculate the normalized score. We opt for the formula post-test rating minus pre-test divided by post-test rating plus pre-test.

This formula can be interpreted as a way to normalize the change between the post and baseline values relative to their combined magnitude. This normalization method balances the difference with the total magnitude of the measures, providing a proportionate understanding of the change.

Range: The output will always be between -1 and 1. If the post value is much larger than the baseline, the result approaches 1. If the post value is much smaller than the baseline, the result approaches -1. If the post value is equal to the baseline, the result is 0.

Thus, the fist step is to calculate the normalized score for each question and condition

VAS <- VAS %>%
  mutate(fatigue_EXP = (FAT_POS_EXP - FAT_PRE_EXP) / (FAT_POS_EXP + FAT_PRE_EXP))

VAS <- VAS %>%
  mutate(fatigue_CON = (FAT_POS_CON -FAT_PRE_CON) / (FAT_POS_CON + FAT_PRE_CON))

VAS <- VAS %>%
  mutate(boredom_EXP = (BOR_POS_EXP - BOR_PRE_EXP) / (BOR_POS_EXP + BOR_PRE_EXP))

VAS <- VAS %>%
  mutate(boredom_CON = (BOR_POS_CON - BOR_PRE_CON) / (BOR_POS_CON + BOR_PRE_CON))

VAS <- VAS[,c(1,10:13)]

#It consider Na some values that should be 0
VAS[c(2, 11,16), 5] = 0


VAS_long <- pivot_longer(VAS, cols =2:5, names_to = c("Scale","Condition"),
                         names_sep ="_", values_to = "Value", values_drop_na = TRUE)

and then, the normalized score looks like this:

Fatigue <- filter(VAS_long, Scale == "fatigue")


FatigueRain <- ggplot(data = Fatigue, 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 ="Subjective level of mental fatigue",
        caption = "n = 21")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size



Boredom <- filter(VAS_long, Scale == "boredom")

BoredomRain <- ggplot(data = Boredom, 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 ="Subjective level of boredom",
        caption = "n = 20")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size


VASplots <- ggarrange(FatigueRain, BoredomRain,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1, 
                     # common.legend = T,
                     # legend = "right",
                     align = "none"
)

VASplots

What do the Bayes tests say about these scales?

For the mental fatigue question:

bf_result_fat <- ttestBF(x = VAS$fatigue_EXP, y = VAS$fatigue_CON, paired = TRUE)
bf_result_fat
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 16.10336 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = VAS$fatigue_EXP, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = VAS$fatigue_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 <- Fatigue %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value),
    n = n()
  )

summary_stats
## # A tibble: 2 × 4
##   Condition  mean    sd     n
##   <chr>     <dbl> <dbl> <int>
## 1 CON       0.189 0.379    21
## 2 EXP       0.495 0.307    21
# Print the credible intervals
print("Credible Interval for VAS Fatigue EXP:")
## [1] "Credible Interval for VAS Fatigue EXP:"
print(credible_interval_exp)
##      2.5%     97.5% 
## 0.3327238 0.6204797
print("Credible Interval for VAS Fatigue CON:")
## [1] "Credible Interval for VAS Fatigue CON:"
print(credible_interval_con)
##        2.5%       97.5% 
## 0.009497868 0.334281498
# Frequentist approach
t.test(x = VAS$fatigue_EXP, y = VAS$fatigue_CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  VAS$fatigue_EXP and VAS$fatigue_CON
## t = 3.4474, df = 20, p-value = 0.002547
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.1206847 0.4904924
## sample estimates:
## mean difference 
##       0.3055885

And for the boredom question:

VAS <- VAS %>% drop_na(boredom_CON)

bf_result_boredom <- ttestBF(x = VAS$boredom_EXP, y = VAS$boredom_CON, paired = TRUE)
bf_result_boredom
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 77.67717 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = VAS$boredom_EXP, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = VAS$boredom_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 <- Boredom %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value),
    n = n()
  )

summary_stats
## # A tibble: 2 × 4
##   Condition  mean    sd     n
##   <chr>     <dbl> <dbl> <int>
## 1 CON       0.212 0.574    20
## 2 EXP       0.704 0.391    21
# Print the credible intervals
print("Credible Interval for VAS Boredom EXP:")
## [1] "Credible Interval for VAS Boredom EXP:"
print(credible_interval_exp)
##      2.5%     97.5% 
## 0.4850422 0.8735760
print("Credible Interval for VAS Boredom CON:")
## [1] "Credible Interval for VAS Boredom CON:"
print(credible_interval_con)
##        2.5%       97.5% 
## -0.06170937  0.43894700
t.test(x = VAS$boredom_EXP, y = VAS$boredom_CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  VAS$boredom_EXP and VAS$boredom_CON
## t = 4.2607, df = 19, p-value = 0.0004226
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.2526419 0.7405339
## sample estimates:
## mean difference 
##       0.4965879

At subjective level with the VAS questions, the evidence is strong in favor of the alternative hypothesis that performing the mental effort task until task failure increase the level of mental fatigue and boredom. Then, If there were any effect in the subsequent physical performance, it could not be solely attribute to mental fatigue.

Temporal experience tracing

Intro

In this part, were are analyzing the subjective experience across the tasks (either the cognitive task until task failure or the documentary). Participants had to draw how they felt for 6 subjective dimensions (Perceived Effort, Mind-wandering, Boredom, Concentration, Mental Fatigue and Frustration). See an example:

TET

Pre-processing

I pre-processed the images in Matlab. All images were vectorized and double checked to ensure that the traces represented the draw of the participants. The data points for each participant in the experimental condition was individualized based on the total duration of the task. Each data point was considered as 30 seconds of the trace. That’s to say, if one participant performed the task for 1000 seconds, we divided this for 30, so it will have 33 data points. This assuming that participants can recall they subjective expericience with a precision of 30 seconds across time. The participants were not aware of the time the spent in the task.

For the control condition, we obtained the same number of data points, assuming a average duration of the documentary of 70 minutes.

First step is to organize the data from the matrix (Nrow x 6 Columns) obtained in Matlab after processing all the images. I did not have the variables “epoch” and “condition, thus I add this in this step.

set.seed(123)
#Read data

invisible({capture.output({

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


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

})})

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


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

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


# Control condition
TET_con$Condition <- rep("con", 17640) #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_exp, TET_con)

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

K-means cluster

A “cluster” refers to a group of data points or observations that are similar to each other and dissimilar to data points in other clusters. Clustering is a process of grouping a set of objects in such a way that objects in the same group (cluster) are more similar to each other than to those in other groups. Similarity in the context of clustering refers to how alike or comparable two data points are to each other.

The metric used for calculating the similarity is the Euclidean Distance. It calculates the straight-line distance between two data points in the feature space. It is defined as the square root of the sum of squared differences between corresponding elements of two vectors.

K-means minimizes the sum of squared distances within each cluster, effectively making it equivalent to minimizing the variance within clusters. This objective function is optimized iteratively by alternating between two steps:

Assignment Step (Expectation Step): Assign each data point to the nearest cluster center (centroid) based on the squared Euclidean distance.

Update Step (Maximization Step): Recalculate the cluster centroids by taking the mean of all data points assigned to each cluster.

The data needs to be organized as follows, in wide format:

knitr::kable(head(TET_merged[c(1:3, 500:503), 1:9]), "pipe", caption = "Data organization for the clustering")
Data organization for the clustering
Participant Condition Epoch Effort Mind-wandering Boredom Concentration Fatigue Frustration
P1 exp 1 0.0151363 0.8588650 0.5643682 0.2195388 0.1677642 0.0551477
P1 exp 2 0.0316939 0.8610704 0.5543321 0.2200673 0.1490859 0.0536380
P1 exp 3 0.0362245 0.8183531 0.5460332 0.2206294 0.1484836 0.0446818
P7 exp 76 0.4259803 0.2705053 0.3783038 0.3530884 0.4092780 0.1181638
P7 exp 77 0.4332258 0.2661568 0.3844858 0.3605145 0.4160527 0.1187221
P7 exp 78 0.4404091 0.2618829 0.3914838 0.3680246 0.4234085 0.1195441

We need to identify the optimal number of cluster that we want to use for the algorithm. The choice could be driven by Domain Knowledge: Utilize domain knowledge or context-specific information to decide the appropriate number of clusters. Or comparing the output of different techniques. As it is the first time we have used this task until failure, we do not have a true hypothesis of the number of cluster. In addition, we have two conditions with two clear difference in terms of cognitive load.

Independently of the number of cluster, we could expect two clears pattern in the data.

  1. If the control condition actually is able to keep a neutral emotional state, we should expect the representation of only one cluster across time, indicating that there are not fluctuation in the emotional state.

  2. In experimental condition, in the experimental condition, we should expect a transition from one cluster to another across time, characterized at the beginning with low fatigue, boredom, mind wandering, frustration and perceived effort, to change to into of a state of the contrary emotional state

Then, the main methods to identify the optimal number of clusters are:

Elbow Method Plot the within-cluster sum of squares (WCSS) or total within-cluster variation against the number of clusters. The “elbow” point, where the rate of decrease in WCSS slows down, can be considered as the optimal number of clusters.

Silhouette Method Calculate the average silhouette width for different numbers of clusters. The number of clusters with the highest average silhouette width is considered optimal.

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

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

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

We can also use the Caliński-Harabasz score to choose the optimal number of clusters:

# 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 = TET_merged[, 4:9])

# 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 = 4059.298 
## k = 3 : CH Score = 2967.086 
## k = 4 : CH Score = 2772.73 
## k = 5 : CH Score = 2366.189 
## k = 6 : CH Score = 2153.045 
## k = 7 : CH Score = 1959.114 
## k = 8 : CH Score = 1737.241 
## k = 9 : CH Score = 1679.404 
## k = 10 : CH Score = 1678.547
# 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")

In the end, everything seems like it could be two clusters the ideal scenario. The Caliński-Harabasz Score establish we should pick the cluster solution with the higher score.

Let’s check with 2 clusters:

km.out <- kmeans(TET_merged[, 4:9], centers=2,nstart=100)
km.out$centers
##      Effort Mind-wandering   Boredom Concentration   Fatigue Frustration
## 1 0.2028767      0.6439305 0.1747040     0.2045217 0.1260444   0.2186864
## 2 0.7306542      0.5134694 0.5995464     0.6350621 0.5015579   0.5177504
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(Effort, 'Mind-wandering', Boredom, Concentration, Fatigue, Frustration), 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"))

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_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))



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


# 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

Are these dimensions statistically different 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", "Frustration")

# 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.964393e+09
Mind-wandering Mind-wandering 1.871980e+00
Boredom Boredom 7.450041e+10
Concentration Concentration 3.888439e+13
Fatigue Fatigue 1.342122e+07
Frustration Frustration 6.526674e+04

This is how it looks like, independently of the conditions. We see that the intensity of each dimension is different.

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 38 0.2682513 0.1845228 0.0299336 0.2076001 0.3289024 Effort
2 27 0.6806289 0.1918180 0.0369154 0.6047482 0.7565096 Effort
1 38 0.6262567 0.2782615 0.0451400 0.5347945 0.7177190 Mind-wandering
2 27 0.4883874 0.2044461 0.0393457 0.4075112 0.5692636 Mind-wandering
1 38 0.2272029 0.1590593 0.0258028 0.1749214 0.2794844 Boredom
2 27 0.6198066 0.1677619 0.0322858 0.5534421 0.6861710 Boredom
1 38 0.2247825 0.1581614 0.0256572 0.1727961 0.2767688 Concentration
2 27 0.6520357 0.1389160 0.0267344 0.5970823 0.7069890 Concentration
1 38 0.1612134 0.1322631 0.0214559 0.1177396 0.2046871 Fatigue
2 27 0.5011429 0.2402706 0.0462401 0.4060950 0.5961908 Fatigue
1 38 0.2315332 0.1795249 0.0291228 0.1725248 0.2905416 Frustration
2 27 0.5088248 0.1986850 0.0382369 0.4302277 0.5874220 Frustration

*Note that counts means the number of data points identified in a cluster in the whole sample of participants. That to say, there are 27 counts for the cluster 2 for the 21 participants in the experimental conditions, and 6 participants presented at some point values considered as cluster 2, for example.

However, we can hypothesize that cluster 2 will be more present during the execution of the task, and cluster 1 would be more likely in the documentary.

# 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 == "exp")
  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)

num_bins <- Task %>% 
    mutate(Epochs = Task/30)

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


plot_stacked_exp <- ggplot(exp_relative, aes(x = 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
  geom_vline(data = num_bins, aes(xintercept = Epochs, y = 0), color = "black", size = 0.5, linetype="dotted") #add vertical line to point out the contribution of each participant
  

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 = "Duration of documentary",
       ) +
  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
  


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

Note that for the experimental condition, the length of the task is different for each participant, and the total contribution is reducing across time. The lined vertical lines delimits the contribution of each participant.

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
##   exp  3948  8940
##   con 14436  3204
# Extract counts

exp_1 <- contingency_table["exp", "1"]
exp_2 <- contingency_table["exp", "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_1 / con_2

# 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): 2.264438
cat("Odds for control condition (Cluster 1/Cluster 2):", odds_con, "\n")
## Odds for control condition (Cluster 1/Cluster 2): 4.505618
# 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 < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.09283757 0.10342500
## sample estimates:
## odds ratio 
## 0.09803481
# 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 = 8148.3, df = 1, p-value < 2.2e-16

In summary:

  1. these two conditions represent different emotional states

  2. In the experimental condition, performing a cognitive task until failure is more likely to be associated with an subjective state of high fatigue, boredom, frustration, sensation of effort, concentration and mind-wandering. However, this state is still dynamic as the contrary cluster of sensation is also present during the task.

  3. This approach of the self-chosen documentary seems effective to maintain an stable emotional state, as the odd ratio indicate that the cluster 1 is 4.5 times more likely to happens than the cluster 2 across the session.

Eye tracking

Here we have two levels of analysis, as in the control condition, participants completed only the task for one block before and after.

  1. We have the normalized score comparison between the two conditions. For being able to compare, we take the first minute and the last minute of the execution of the task for the experimental condition.

  2. We have the analysis of the time on task effect during the execution of the task. The duration of the task is different for each individual, but we split the continuous value according to the number of epochs generated for the subjective TET data.

We have there main parameters: pupil size, fixation (within the area of interest) and ¿blink rate?

Confirmatory Pre-Post score both conditions

Pupil Size

First, we need to calculate from the task file, the average from the first and last “block” of the task as the pre and post value. we need also to calculate the average from the control and then merge the two dataset together. Finally, we calculate the normalize score as before.

  # Read data from OSF.
invisible({capture.output({

# url <- 'https://osf.io/mgbn2//?action=download'
# filename <- "Pupil_exp.txt"
# GET(url, write_disk(filename, overwrite = TRUE))
# Pupil_exp <- read.delim(filename)


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


url <- 'https://osf.io/6mnaq//?action=download'
filename <- "Gaze_Pupil_20_08.txt"
GET(url, write_disk(filename, overwrite = TRUE))
gaze_pupil <- read.delim(filename)


})})


gaze_pupil$LEFT_PUPIL_SIZE <- as.numeric(gaze_pupil$LEFT_PUPIL_SIZE)
Pupil_exp <- filter(gaze_pupil, LEFT_PUPIL_SIZE !="Na")
Pupil_exp <- gaze_pupil[, c(1:2,5)]
Pupil_exp[Pupil_exp == 0] <- NA
Pupil_exp <- filter(Pupil_exp, LEFT_PUPIL_SIZE !="Na")


nove <- filter(Pupil_exp, Participant =="9")
  

Pupil_exp$Participant <- as.numeric(Pupil_exp$Participant)


Pupil_exp <- Pupil_exp %>% 
group_by(Participant) %>%
  mutate(Sequence = row_number()) %>%
  ungroup() 

n_select <- 30000


# # Group by Participant and select the first and last 30 000 rows. Approx first and las minute
 Pupil_exp <- Pupil_exp %>%
    group_by(Participant) %>%
   filter(Sequence <= n_select | Sequence > (max(Sequence) - n_select)) %>%
   ungroup() %>%  # Remove grouping to get a clean data frame
   select(-Sequence)  # Optionally, remove the Sequence column if not needed

sequence = c("Pre", "Pos")
Period=rep(sequence, each = 30000) #
Pupil_exp$Period <- rep(Period, 21) # 21 files with Fixations

pupil_confirmatory_exp <- Pupil_exp %>%
  group_by(Participant, Period) %>%
  summarise(
    EXP = mean(LEFT_PUPIL_SIZE, na.rm = TRUE))

pupil_confirmatory_exp <- pupil_confirmatory_exp %>% 
    spread(key = Period, value = EXP)


pupil_confirmatory_exp <- pupil_confirmatory_exp %>% 
  mutate(
  EXP = (Pos - Pre) / (Pos + Pre)
  )
pupil_confirmatory_exp <- pupil_confirmatory_exp[, c(1,4)]


# Pupil_exp$Pupil <- as.numeric(Pupil_exp$Pupil)

Pupil_con$Pupil <- as.numeric(Pupil_con$Pupil)
Pupil_con$Participant <- as.numeric(Pupil_con$Participant)
Pupil_con$Period <- as.factor(Pupil_con$Period)


# #Calculate pre and post for the pupil size for the experimental condition
# Pupil_exp <- filter(Pupil_exp, Participant != "13_Task")
# 
# 
# 
# pupil_confirmatory_exp <- Pupil_exp %>%
#   group_by(Participant) %>%
#   arrange(Participant, Trial) %>%  # Sort by Participant and trial
#   summarize(
#     Pre = mean(Pupil[1:3000], na.rm = TRUE),
#     Pos = mean(tail(Pupil, 3000), na.rm = TRUE)
#   )
# 
# pupil_confirmatory_exp <- pupil_confirmatory_exp %>%
#   mutate(EXP = (Pos - Pre) / (Pos + Pre))
# pupil_confirmatory_exp <- pupil_confirmatory_exp[, c(1,4)]


pupil_confirmatory_exp <- gaze_pupil %>%
  group_by(Participant) %>%
  summarise(
    Pre = mean(LEFT_PUPIL_SIZE[TRIAL_INDEX == min(TRIAL_INDEX)], na.rm = TRUE),
    Pos = mean(LEFT_PUPIL_SIZE[TRIAL_INDEX == max(TRIAL_INDEX)], na.rm = TRUE)
  )


pupil_confirmatory_exp <- pupil_confirmatory_exp %>%
  mutate(EXP = (Pos - Pre) / (Pos + Pre))
pupil_confirmatory_exp <- pupil_confirmatory_exp[, c(1,4)]

                          # -----------------------------------#

#Pupil Control condition

pupil_confirmatory_con <- Pupil_con %>%
 group_by(Participant, Period) %>%
  summarise(AveragePupil = mean(Pupil, na.rm = TRUE)) %>%
  ungroup()

pupil_confirmatory_con <- pupil_confirmatory_con %>%
  pivot_wider(names_from = Period, values_from = AveragePupil)

pupil_confirmatory_con <- pupil_confirmatory_con %>%
  mutate(CON = (Pos - Pre) / (Pos + Pre))
pupil_confirmatory_con <- pupil_confirmatory_con[, c(1,4)]


#Final data set pupil

Pupil <- merge(pupil_confirmatory_con, pupil_confirmatory_exp, by = "Participant")


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

The visualization is:

PupilRain <- ggplot(data = Pupil_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 ="Pupil size",
        caption = "n = 20")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

PupilRain

A reminder, a negative value means that the post value is lower than the pre.

The bayes factor shows an extreme evidence in favor of the alternative hypothesis:

Pupil <- Pupil %>% drop_na(CON)
bf_result_pupil <- ttestBF(x = Pupil$EXP, y = Pupil$CON, paired = TRUE)
bf_result_pupil
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 46.69359 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = Pupil$EXP, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = Pupil$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 <- Pupil_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value),
    n = n()
  )

summary_stats
## # A tibble: 2 × 4
##   Condition    mean     sd     n
##   <chr>       <dbl>  <dbl> <int>
## 1 CON       -0.0205 0.0315    20
## 2 EXP       -0.0657 0.0564    21
# Print the credible intervals
print("Credible Interval for Pupil Exp:")
## [1] "Credible Interval for Pupil Exp:"
print(credible_interval_exp)
##        2.5%       97.5% 
## -0.09132687 -0.03772185
print("Credible Interval for Pupil Con:")
## [1] "Credible Interval for Pupil Con:"
print(credible_interval_con)
##         2.5%        97.5% 
## -0.032983583 -0.004358865
t.test(x = Pupil$EXP, y = Pupil$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  Pupil$EXP and Pupil$CON
## t = -4.0061, df = 19, p-value = 0.0007555
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.07335876 -0.02301042
## sample estimates:
## mean difference 
##     -0.04818459

Fixation

Here, we focus weather the participants were looking at the stimulus when they appeared on the screen or not. We define a area of interest around the stimulus leaving a marge around this area.

We calculate the total numbers of fixation within this area and then the relative in relation to the total in a given period of time.

invisible({capture.output({


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

})})

# For the Control condition

Fixation_con$Fix_X <- as.numeric(Fixation_con$Fix_X)
Fixation_con$Fix_Y <- as.numeric(Fixation_con$Fix_Y)

#Create condition to establish if the current fixation is inside or outside of the area of interest

# Define the bounds according to the screen (with respect to the area of the simulus)
# x_min <- 960 - 74 / 2   # 923
# x_max <- 960 + 74 / 2   # 997
# y_min <- 600 - 93 / 2   # 553.5
# y_max <- 600 + 93 / 2   # 646.5

# Define the bounds according to the screen (leaving an extra marge of 14 pxs)
x_min <- 960 - 88 / 2   # 923
x_max <- 960 + 88 / 2   # 997
y_min <- 600 - 107 / 2   # 553.5
y_max <- 600 + 107 / 2   # 646.5


#Establish which are inside and which are outside
Fixation_con <- Fixation_con %>%
  mutate(fixation = ifelse(Fix_X >= x_min & Fix_X <= x_max & Fix_Y >= y_min & Fix_Y <= y_max, "Inside", "Outside"))


#  Count total fixations and "Inside" fixations for each Participant and Period
Fixation_count_con <- Fixation_con %>%
  group_by(Participant, Period) %>%
  summarise(
    TotalFixations = n(),
    OutsideFixations = sum(fixation == "Outside"),
    .groups = 'drop'
  )

# Calculate the relative number of "Inside" fixations
Fixation_count_con <- Fixation_count_con %>%
  mutate(fixation = OutsideFixations / TotalFixations)
Fixation_count_con <- Fixation_count_con[, c(1,2,5)]

# Calculate the normalized score

Fixation_count_con <- filter(Fixation_count_con, Participant != "10") #Missing Pre value for this participant.

# adapt to wide to calculate normalized score
Fixation_count_con <- Fixation_count_con %>%
  pivot_wider(names_from = Period, values_from = fixation)

# Normalized score
Fixation_count_con <- Fixation_count_con %>%
  mutate(CON = (Pos - Pre) / (Pos + Pre))

Fixation_count_con <- Fixation_count_con[, c(1,4)]

#It consider Na a value values that should be 0
Fixation_count_con[2, 2] = 0

                      ####### For the Experimental condition #####


# Obtain the pre and post for the experimental condition

# Fixation_exp <- Fixation_exp %>%
# group_by(Participant) %>%
#   mutate(Sequence = row_number()) %>%
#   ungroup()

gaze_pupil$LEFT_GAZE_X <- as.numeric(gaze_pupil$LEFT_GAZE_X)

Fixation_exp <- filter(gaze_pupil, LEFT_GAZE_X !="Na")


Fixation_exp <- Fixation_exp %>% 
group_by(Participant) %>%
  mutate(Sequence = row_number()) %>%
  ungroup() 

Fixation_exp$LEFT_GAZE_Y <- as.numeric(Fixation_exp$LEFT_GAZE_Y)
Fixation_exp$Participant <- as.numeric(Fixation_exp$Participant)
Fixation_exp <- Fixation_exp[, c(1:4,7)]


n_select <- 30000 #Approximately one block


# # Group by Participant and select the first and last 30 000 rows. Approx first and last minute
 Fixation_exp <- Fixation_exp %>%
    group_by(Participant) %>%
   filter(Sequence <= n_select | Sequence > (max(Sequence) - n_select)) %>%
   ungroup() %>%  # Remove grouping to get a clean data frame
   select(-Sequence)  # Optionally, remove the Sequence column if not needed

Fixation_exp <- Fixation_exp %>%
  mutate(fixation = ifelse(LEFT_GAZE_X >= x_min & LEFT_GAZE_X <= x_max & LEFT_GAZE_Y >= y_min & LEFT_GAZE_Y <= y_max, "Inside", "Outside"))

# It was an alternative to take the first and last trial of the sample
# However, the last trial could have much lower number of samples if the participant stopped the task.

# Fixation_exp <- Fixation_exp %>%
#   group_by(Participant) %>%
#   summarise(
#     Pre = mean(fixation[TRIAL_INDEX == min(TRIAL_INDEX)] == "Outside", na.rm = TRUE),
#     Pos = mean(fixation[TRIAL_INDEX == max(TRIAL_INDEX)] == "Outside", na.rm = TRUE)
#   )


# Introduce pre and pos label
sequence = c("Pre", "Pos")
Period=rep(sequence, each = 30000) #
Fixation_exp$Period <- rep(Period, 21) # 21 files with Fixations
 

 # Count total fixations and "Outise" fixations for each Participant and Period
Fixation_count_exp <- Fixation_exp %>%
  group_by(Participant, Period) %>%
  summarise(
    TotalFixations = n(),
    OutsideFixations = sum(fixation == "Outside"),
    .groups = 'drop'
  )


# Calculate the relative number of "Outside" fixations
Fixation_count_exp <- Fixation_count_exp %>%
  mutate(fixation = OutsideFixations / TotalFixations)
Fixation_count_exp <- Fixation_count_exp[, c(1,2,5)]

# Calculate the normalized score

# adapt to wide to calculate normalized score
Fixation_count_exp <- Fixation_count_exp %>%
  pivot_wider(names_from = Period, values_from = fixation)

# Normalized score
Fixation_count_exp <- Fixation_count_exp %>%
  mutate(EXP = (Pos - Pre) / (Pos + Pre))

Fixation_count_exp <- Fixation_count_exp[, c(1,4)]

# Merge the two dataset

Fixation_count_exp$Participant <-as_numeric(Fixation_count_exp$Participant) 
Fixations <- merge(Fixation_count_con, Fixation_count_exp, by = "Participant")


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

The visualization is:

FixationRain <- ggplot(data = Fixations_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 ="Fixations ",
        caption = "n = 20")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

FixationRain

The bayes factor shows an extreme evidence in favor of the alternative hypothesis:

bf_result_fixations <- ttestBF(x = Fixations$EXP, y = Fixations$CON, paired = TRUE)
bf_result_fixations
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 2.968834 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = Fixations$EXP, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = Fixations$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 <- Fixations_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value),
    sd = sd(Value),
    n = n()
  )

summary_stats
## # A tibble: 2 × 4
##   Condition    mean    sd     n
##   <chr>       <dbl> <dbl> <int>
## 1 CON       -0.0670 0.556    20
## 2 EXP        0.326  0.439    20
# 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.09710005 0.50094877
print("Credible Interval for Fixation Con:")
## [1] "Credible Interval for Fixation Con:"
print(credible_interval_con)
##       2.5%      97.5% 
## -0.3032550  0.1818677
t.test(x = Fixations$EXP, y = Fixations$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  Fixations$EXP and Fixations$CON
## t = 2.5546, df = 19, p-value = 0.01937
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.07097313 0.71461147
## sample estimates:
## mean difference 
##       0.3927923

Time on task and explorative analysis

The second step is to explore the time on task effect of the eye tracking measures. As we will later match number of data points from the subjective scales and the eye-tracking, we create a new variable (bins) in the data set to divide into blocks/epochs

Pupil size across time

I want to check if the pupil size is different if we only consider the samples with the fixation inside the area of interest

gaze_pupil$LEFT_GAZE_X <- as.numeric(gaze_pupil$LEFT_GAZE_X)
gaze_pupil$LEFT_GAZE_Y <- as.numeric(gaze_pupil$LEFT_GAZE_Y)
gaze_pupil$LEFT_PUPIL_SIZE <- as.numeric(gaze_pupil$LEFT_PUPIL_SIZE)
gaze_pupil$LEFT_IN_BLINK <- as.numeric(gaze_pupil$LEFT_IN_BLINK)


# Creating the bins from the total duration of each participant completing the task.
# each 30s, 1 epoch
num_bins <- Task %>% 
    mutate(Epochs = Task/30)


# Calculate the fixations inside or outside of the area of interest
gaze_pupil <- gaze_pupil %>%
  mutate(fixation = ifelse(LEFT_GAZE_X >= x_min & LEFT_GAZE_X <= x_max & LEFT_GAZE_Y >= y_min & LEFT_GAZE_Y <= y_max, "Inside", "Outside"))


pupil_inside <- filter(gaze_pupil, fixation =="Inside")

# to create the bins
pupil_inside <- pupil_inside %>%
  left_join(num_bins, by = "Participant")
pupil_inside$Participant <- as.numeric(pupil_inside$Participant)


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

# Create the epochs for each participant
binned_Pupil_inside <- pupil_inside %>%
  group_by(Participant) %>%
  summarise(
    Epoch = 1:unique(Epochs),  # Create bin index
    Pupil = calculate_bin_averages(LEFT_PUPIL_SIZE, unique(Epochs))  # Calculate the bin averages
  ) %>%
  ungroup()

The average pupil size across epochs for each participant

pupil_relative_change_inside <- binned_Pupil_inside %>%
group_by(Participant) %>%
  mutate(
    baseline = first(Pupil[Epoch == 1]),  # Get baseline value (Pupil size at Trial 1)
    relative_change = ifelse(Epoch == 1, 0, (Pupil - baseline) / baseline * 100)  # Set relative change to 0 for Trial 1
  ) %>%
  ungroup() # Remove grouping) 


# Calculate the average of each epoch to add to the final plot
average_data <- pupil_relative_change_inside %>%
  group_by(Epoch) %>%
  summarise(average_change = mean(relative_change, na.rm = TRUE))


relativeplot_inside <- ggplot(pupil_relative_change_inside, aes(x = Epoch, y = relative_change)) +
  geom_line(aes(color = Participant), alpha = 0.5) +  # Slightly transparent individual lines
  geom_line(data = average_data, aes(x = Epoch, y = average_change), color = "red", size = 1.5) +  # Thicker average line
  labs(
    title = "Relative change in pupil size over time",
    y = "Relative Change (%)"
  ) +
   sm_minimal()+
  easy_remove_x_axis(what = c("text"), teach = FALSE)+
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  geom_vline(data = num_bins, aes(xintercept = Epochs, y = 0), color = "black", size = 0.7, linetype="dotted")


relativeplot_inside

The vertical dotted lines point out the time each participant stopped the mental effort task to failure.

Fixations across time

# Add number of epochs
Fixation_clean <- gaze_pupil %>%
  left_join(num_bins, by = "Participant")


Fixation_clean <- Fixation_clean[, c(1,7,9)]

Fixation_clean <- filter(Fixation_clean, fixation !="Na")

# Function to calculate binned averages of "Outside" fixation
calculate_bin_averages <- function(fixation_data, num_bins) {
  # Determine if each fixation is "Outside"
  is_outside <- fixation_data == "Outside"
  
  # Calculate the bin index for each data point
  bin_indices <- cut(seq_along(is_outside), breaks = num_bins, labels = FALSE)
  
  # Calculate the average (proportion of "Outside") for each bin
  bin_means <- tapply(is_outside, bin_indices, mean)
  
  return(bin_means)
}

# Apply binning and averaging for each participant
binned_fixation_clean <- Fixation_clean %>%
  group_by(Participant) %>%
  summarise(
    Epoch = 1:unique(Epochs),  # Create bin index
    Relative_Outside = calculate_bin_averages(fixation, unique(Epochs))  # Calculate the bin averages
  ) %>%
  ungroup()

# Calculate average value across time to include in the figure.
average_data <- binned_fixation_clean %>%
  group_by(Epoch) %>%
  summarise(average_change = mean(Relative_Outside, na.rm = TRUE))



# Plot of 'Outside' fixation proportion over time for each participant
Fixation_TOT <- ggplot(binned_fixation_clean, aes(x = Epoch, y = Relative_Outside, )) +
    geom_line(aes(color = Participant), alpha = 0.5) +  # Slightly transparent individual lines
    geom_line(data = average_data, aes(x = Epoch, y = average_change), color = "blue", size = 1.5) +  # Thicker average line
  labs(title = " Proportion of 'Outside' fixations over time",
       x = "Epochs",
       y = "Relative outside fixation") +
  sm_minimal()+
    geom_vline(data = num_bins, aes(xintercept = Epochs, y = 0), color = "black", size = 0.7, linetype="dotted")

Fixation_TOT

blinks <- gaze_pupil[, c(1,6)]

blinks <- blinks %>%
  left_join(num_bins, by = "Participant")


# Function to calculate binned averages of blinks
calculate_bin_counts <- function(blink_data, num_bins) {
  # Determine if each data is "Blink"
  is_blink <- blink_data == "1"
  
  # Calculate the bin index for each data point
  bin_indices <- cut(seq_along(is_blink), breaks = num_bins, labels = FALSE)
  
  # Calculate the average (proportion of "Outside") for each bin
  bin_counts <- tapply(is_blink, bin_indices, sum)
  
  return(bin_counts)
}

# Apply binning and averaging for each participant
binned_blinks <- blinks %>%
  group_by(Participant) %>%
  summarise(
    Epoch = 1:unique(Epochs),  # Create bin index
    Blink = calculate_bin_counts(LEFT_IN_BLINK, unique(Epochs))  # Calculate the bin averages
  ) %>%
  ungroup()

# Calculate average value across time to include in the figure.
average_data <- binned_blinks %>%
  group_by(Epoch) %>%
  summarise(average_change = mean(Blink, na.rm = TRUE))

# Plot of 'Outside' fixation proportion over time for each participant
Blink_TOT <- ggplot(binned_blinks, aes(x = Epoch, y = Blink, )) +
    geom_line(aes(color = Participant), alpha = 0.5) +  # Slightly transparent individual lines
    geom_line(data = average_data, aes(x = Epoch, y = average_change), color = "blue", size = 1.5) +  # Thicker average line
  labs(title = "'Blink rate Over Time",
       x = "Epoch",
       y = "Blink rate") +
  sm_minimal()+
    geom_vline(data = num_bins, aes(xintercept = Epochs, y = 0), color = "black", size = 0.7, linetype="dotted")

Blink_TOT

invisible({capture.output({


url <- 'https://osf.io/k7xjv//?action=download'
filename <- "Blinks.txt"
GET(url, write_disk(filename, overwrite = TRUE))
blink <- read.delim(filename)


})})

Cognitive task

invisible({capture.output({

url <- 'https://osf.io/gkra7//?action=download'
filename <- "Performance_task_exp.txt"
GET(url, write_disk(filename, overwrite = TRUE))
Task_perf <- read.delim(filename)

})})

Task_perf <- Task_perf[, c(1,5)]

Task_perf$Participant <- as.numeric(Task_perf$Participant)

#Divide in Bins

# Merge the bin information into the main data frame
Task_perf <- Task_perf %>%
  left_join(num_bins, by = "Participant")


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

# Apply binning and averaging for each participant
binned_Task <- Task_perf %>%
  group_by(Participant) %>%
  summarise(
    Epoch = 1:unique(Epochs),  # Create bin index
    Performance = calculate_bin_averages(PERFORMANCE, unique(Epochs))  # Calculate the bin averages
  ) %>%
  ungroup()

binned_Task <- filter(binned_Task, Participant !="9")
binned_Task$Participant <- paste0("P", binned_Task$Participant)


TET_filter <- filter(TET_filter, Participant != "P9")


TET_perf <- merge(TET_filter, binned_Task, by = c("Participant", "Epoch"))

TET_perf <- TET_perf[, c(1:2, 10:11)]

Visualization of performance by cluster. This is to compare the mean performance of both cluster

# ggplot(TET_perf, aes(x = Performance, as.factor(Cluster))) +
#   geom_density(alpha = 0.5) +
#   labs(title = "Density Plot of Performance by cluster",
#        x = "Accuracy",
#        y = "Density")+
#   scale_colour_brewer(palette = "Set2")+
#   scale_fill_brewer(palette = "Set2")+ 
#   theme(plot.background = element_rect(colour = NA)) +
#   size


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

PerfTET <- ggplot(data = TET_perf, mapping = aes(x = Cluster, y = Performance)) +
  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 = Cluster, y = Performance, fill = Cluster),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  # geom_point(aes(fill = Cluster, group=Participant), size=2,shape=21)+
  # geom_line(aes(group=Participant),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs(y = "Performance (%)", x = element_blank(),
        
        )+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  scale_x_discrete(labels = custom_x_labels) +   # Custom x-axis labels
  size

PerfTET

TET_perf_long <- TET_perf %>% 
  pivot_wider(
    names_from  = c(Cluster), 
    values_from = c(Performance)
  )


ttestBF(formula = Performance ~ Cluster, data = TET_perf)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 1.177398 ±0.02%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS
Cluster1 <- na.omit(data.frame(Cluster1 = TET_perf_long$`1`))
Cluster2 <- na.omit(data.frame(Cluster2 = TET_perf_long$`2`))

bf_Cluster1 <- ttestBF(x = Cluster1$Cluster1, nullInterval = NULL)  # Against 0
bf_Cluster2 <- ttestBF(x = Cluster2$Cluster2, nullInterval = NULL)  # Against 0


# Extract posterior samples for each mean
posterior_cluster1 <- posterior(bf_Cluster1, iterations = 10000)
posterior_cluster2 <- posterior(bf_Cluster2, iterations = 10000)

# Convert to data frames
posterior_cluster1_df <- as.data.frame(posterior_cluster1)
posterior_cluster2_df <- as.data.frame(posterior_cluster2)

# Calculate the 95% credible intervals for each mean
credible_interval_cluster1 <- quantile(posterior_cluster1_df$mu, probs = c(0.025, 0.975))
credible_interval_cluster2 <- quantile(posterior_cluster2_df$mu, probs = c(0.025, 0.975))


# Print the credible intervals
print("Credible Interval for Cluster 1 for Performance:")
## [1] "Credible Interval for Cluster 1 for Performance:"
print(credible_interval_cluster1)
##      2.5%     97.5% 
## 0.7965451 0.8109904
print("Credible Interval for Cluster 2 for Performance:")
## [1] "Credible Interval for Cluster 2 for Performance:"
print(credible_interval_cluster2)
##      2.5%     97.5% 
## 0.7871323 0.7971450
summary_stats <- TET_perf %>%
  group_by(Cluster) %>%
  summarise(
    mean = mean(Performance),
    median= median(Performance),
    sd = sd(Performance),
    n = n()
  )

print(summary_stats)
## # A tibble: 2 × 5
##   Cluster  mean median     sd     n
##   <fct>   <dbl>  <dbl>  <dbl> <int>
## 1 1       0.804  0.823 0.0943   658
## 2 2       0.792  0.799 0.0994  1455
# model <- lmer(Performance ~ Cluster + (1 | Epoch), data = TET_perf)
# 
# # Summary of the model
# summary(model)
# 
# anova_results <- anova(model)
# anova_table <- as.data.frame(anova_results)
# print(anova_table)
# 
 

t_test_result <- t.test(Performance ~ Cluster, data = TET_perf, paired =FALSE)
# print(t_test_result)

The mean performance is lower in the Cluster 2 than in Cluster 1, and it is statistically significant, although the mean performance is almost the same.

ExploratoryTET <- ggarrange(PupilTET, FixationTET,PerfTET,
                     labels = c("A", "B", "C"),
                     ncol = 2, nrow = 2, 
                     # common.legend = T,
                     # legend = "right",
                     align = "none"
)

ExploratoryTET

Physical Performance in the running exercise

We are going first to visualize the performance time of each participant for the two conditions

perf_long <- pivot_longer(perf, cols =2:3, names_to = c("Condition"),
                         values_to = "Performance", values_drop_na = TRUE)

PerformancePlot <- ggplot(data = perf_long, mapping = aes(x = Condition, y = Performance)) +
  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 = Performance, 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 = "Time (s)", x = element_blank(),
        title ="Performance in the time to exhaustion test",
        caption = "n = 21")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

PerformancePlot

Now we perform the main confirmatory analysis and also the analysis we used as stopping rule. Our criteria was:

1)The sample size will be determined by using sequential tests with two-sided Bayes factor with a minimum of 20 participants and controlling the Bayes factor for the main index of physical performance (i.e., reduction of average time in running exercise after the mental effort task) until it reaches strong evidence in favor of the alternative hypothesis (BF10 > 6) or the null hypothesis (BF10 < 1/6).

2)If the Bayes factor does not reach the criteria, we will collect participants in batches of two (to keep the randomization and counterbalancing) until it fulfills the criteria.

3)If not, we also plan to stop the experiment when we reach the maximum number of participants we would be able to compensate economically (40 participants).

4)In any case, even if we were not able to reach the 40 participants sample, we would stop the experiment in July 2024

Thus, the results show:

# Perform the Bayesian paired sample t-test
bf_result_per <- ttestBF(x = perf$EXP, y = perf$CON, paired = TRUE)

bf_result_per
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.2666684 ±0.02%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# # Extract posterior samples (default is 10000 iterations)
# posterior_samples_perf <- posterior(bf_result_per, iterations = 10000)
# 
# # The default parameter extracted in a paired t-test is the effect size delta
# # Convert to a data frame for easier manipulation
# posterior_df_perf <- as.data.frame(posterior_samples_perf)
# 
# # Plot the posterior distribution of the effect size (delta)
# ggplot(posterior_df_perf, aes(x = delta)) +
#   geom_density(fill = "blue", alpha = 0.5) +  # Density plot
#   labs(title = "Posterior Distribution of the Effect Size (delta)",
#        x = "Effect Size (delta)",
#        y = "Density") +
#   theme_minimal()

# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = perf$EXP, nullInterval = NULL)  # Against 0
bf_con <- ttestBF(x = perf$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% 
## 476.4993 658.1403
print("Credible Interval for Perfomance CON:")
## [1] "Credible Interval for Perfomance CON:"
print(credible_interval_con)
##     2.5%    97.5% 
## 495.4362 670.7810
t.test(x = perf$EXP, y = perf$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  perf$EXP and perf$CON
## t = -0.59402, df = 20, p-value = 0.5592
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -66.59973  37.07592
## sample estimates:
## mean difference 
##        -14.7619
summary_stats <- perf_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Performance),
    sd = sd(Performance),
    n = n()
  )

summary_stats
## # A tibble: 2 × 4
##   Condition  mean    sd     n
##   <chr>     <dbl> <dbl> <int>
## 1 CON        589.  187.    21
## 2 EXP        575.  195.    21

The Bayes factor is 0.27 (with a default Cauchy of 0.7), which shows a moderate evidence in favor of the null hypothesis.

However, let’s perform a robustness check to understand how stable the results are as more data is added and provide insights into the robustness of the Bayesian evidence.

# Function to calculate Bayes Factor for given data
calculate_bf <- function(perf) {
  bf <- ttestBF(x = perf$EXP, y = perf$CON, paired = TRUE)
  return(as.numeric(extractBF(bf)$bf))
}

# Initialize an empty data frame to store Bayes Factor values
bf_values <- data.frame(step = integer(), BF = numeric())

# Loop through the data, adding one participant at a time
for (i in 2:nrow(perf)) {
  # Subset the data to include the first i participants
  subset_data <- perf[1:i, ]

  # Calculate the Bayes Factor for the current subset
  current_bf <- calculate_bf(subset_data)

  # Store the results
  bf_values <- rbind(bf_values, data.frame(step = i, BF = current_bf))
}


ggplot(bf_values, aes(x = step, y = BF)) +
  geom_line(color = "blue") +
  geom_point(color = "red") +
  geom_text(aes(label = round(BF, 2)), vjust = -0.5, color = "black", size = 3) +  # Annotate BF values
  labs(title = "Bayes Factor Evolution with Sequential Addition of Participants",
       x = "Number of Participants Included",
       y = "Bayes Factor") +
  theme_minimal()

It seems that the results are robust and the evidence seems to support moderately the null hypothesis with the increment of sample size.

RPE

Given that there were no difference in physical performance and to simplify the analysis, the RPE over time was average across participants

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


summary_stats_rpe <- summarySE(RPE_long, measurevar = "Value",
                       groupvars=c("Participant", "Condition"))
summary_stats_rpe <- summary_stats_rpe[,c(1,2,4)]


summary_stats_rpe <-summary_stats_rpe %>% 
  pivot_wider(
    names_from = "Condition",
    values_from = "Value"
  )



# Perform the Bayesian paired sample t-test
bf_result_rpe <- ttestBF(x = summary_stats_rpe$EXP, y = summary_stats_rpe$CON, paired = TRUE)

bf_result_rpe
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.2717819 ±0.02%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# # Extract posterior samples (default is 10000 iterations)
# posterior_samples_perf <- posterior(bf_result_per, iterations = 10000)
# 
# # The default parameter extracted in a paired t-test is the effect size delta
# # Convert to a data frame for easier manipulation
# posterior_df_perf <- as.data.frame(posterior_samples_perf)
# 
# # Plot the posterior distribution of the effect size (delta)
# ggplot(posterior_df_perf, aes(x = delta)) +
#   geom_density(fill = "blue", alpha = 0.5) +  # Density plot
#   labs(title = "Posterior Distribution of the Effect Size (delta)",
#        x = "Effect Size (delta)",
#        y = "Density") +
#   theme_minimal()

# Perform Bayesian one-sample t-tests for each group
bf_exp <- ttestBF(x = summary_stats_rpe$EXP, nullInterval = NULL)  # Against 0
## t is large; approximation invoked.
bf_con <- ttestBF(x = summary_stats_rpe$CON, nullInterval = NULL)  # Against 0
## t is large; approximation invoked.
# 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.556001 8.189900
print("Credible Interval for RPE CON:")
## [1] "Credible Interval for RPE CON:"
print(credible_interval_con)
##     2.5%    97.5% 
## 7.660380 8.342721
t.test(x = summary_stats_rpe$EXP, y = summary_stats_rpe$CON, paired = TRUE)
## 
##  Paired t-test
## 
## data:  summary_stats_rpe$EXP and summary_stats_rpe$CON
## t = -0.62884, df = 20, p-value = 0.5366
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.5529853  0.2968041
## sample estimates:
## mean difference 
##      -0.1280906
summary_stats <- RPE_long %>%
  group_by(Condition) %>%
  summarise(
    mean = mean(Value)
  )

summary_stats
## # A tibble: 2 × 2
##   Condition  mean
##   <chr>     <dbl>
## 1 CON        7.93
## 2 EXP        7.82