Individualized mental fatigue does not impact neuromuscular function and exercise performance

knitr::kable(head(IMF[, 1:10]), "pipe", caption = "Data for both cognitive task in wide format")

# Change to long format. We merge col 2 to 17 (separated by_) and we keep only 4 cols. Participant, Time, Condition and the variable with the performance. 
IMF_long<- pivot_longer(IMF, cols =2:17, names_to = c("Time", "Condition"),
                       names_sep = "_", values_to = "Performance", values_drop_na = TRUE) 
# Rename the name of the conditions for illustrative purpose
IMF_long<-IMF_long %>% mutate(Condition= recode(Condition, LCL ="Control", HCL= "Mental fatigue"))

knitr::kable(head(IMF_long[, 1:4]), "pipe", caption = "Data for both cognitive task in long format")

#First step is to calculate the average.
# Calculate the average across participants for each time and condition for the global performance data.
IMF_longavg <- IMF_long %>% 
  group_by(Condition, Time) %>% 
  summarise_at(vars(Performance),           
               list(average = mean))              

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

# Create the plot with global performance for each task
IMF.plot <- ggplot(IMF_longavg) +
  aes(x = Time, y = average, colour = Condition, group = Condition) +
  geom_smooth(method = "lm")+
  scale_color_brewer(palette = "Set2", direction = 1) +
  sm_minimal() +  
  guides(fill = "none") +
  labs( y = "Correct answers (%)",
        caption = "n = 22",
        title ="Performance of the individualized and control task") +
  size

# Filter the data set for getting mental fatigue condition
HCL <-  filter(IMF_long, Condition == "Mental fatigue")

HCL.plot <- ggplot(HCL, aes(x = Time, y = Performance, color = Participant, group = Participant)) +
  geom_line(linetype = 1,
            linewidth = 1.2) +
  guides(color = guide_legend(title = "ID"))+
  scale_color_identity(labels = paste("ID", 1:22))+
  sm_minimal() +
  guides(fill = "yes") +
  labs( y = "Performance (%)",
        # caption = "n = 22",
        title ="Individual performance across the task") +
  geom_dl(aes(label = Participant), method = list(dl.combine("first.points", "last.points")), cex = 0.3, alpha = 0.5) +
  size




HCL.plot

# Pass the ggplot to a ggplotly
ggplotly(HCL.plot)

FatigueRain <- ggplot(data = data, mapping = aes(x = Condition, y = VAS_FAT_NOR)) +
  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 = VAS_FAT_NOR, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  geom_point(aes(fill = Condition, group=ID), size=2,shape=21)+
  geom_line(aes(group=ID),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Normalize score (AU)", x = element_blank(),
        title ="Subjective level of mental fatigue",
        caption = "n = 21")+
  size


ActivationRain <- ggplot(data = data, mapping = aes(x = Condition, y = VAS_ARO_NOR)) +
  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 = VAS_ARO_NOR, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  geom_point(aes(fill = Condition, group=ID), size=2,shape=21)+
  geom_line(aes(group=ID),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Normalize score (AU)", x = element_blank(),
        title ="Subjective level of arousal",
        # subtitle = "BF10 = 0.487; Anecdotal evidence in favour of the null hypothesis 
        # (i.e., activation is similar for bot tasks) ",
        caption = "n = 21" ) +
  size

# Arrange the 2 figures in 2 cols and 1 row. Put a common legend on the right side (if desired)
VASplots <- ggarrange(FatigueRain, ActivationRain,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1, 
                     # common.legend = T,
                     # legend = "right",
                     align = "none"
)

VASplots

# Alternatively, we could this to place 2 figures in 2 columns (only for a Rm document)
# `{r out.width=c('50%', '50%'), fig.show='hold'}
# FatigueRain
# ActivationRain

# interactive graph
ggplotly(FatigueRain, width = 600, height= 900)

# Rain cloud  plot
PerformanceRPlot <- ggplot(data = data, 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=ID), size=2,shape=21)+
  geom_line(aes(group=ID),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Time to exhaustion (s)", x = element_blank(),
        title ="Performance",
        caption = "n = 22")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

RPERPlot <- ggplot(data = data, mapping = aes(x = Condition, y = RPE)) +
  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 = RPE, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  sm_minimal() +
  geom_point(aes(fill = Condition, group=ID), size=2,shape=21)+
  geom_line(aes(group=ID),  linewidth=0.3, color='gray', alpha=0.6) +
  guides(fill = "none") +
  labs( y = "Borg scale 1-10 (AU)", x = element_blank(),
        title ="RPE",
        caption = "n = 22")+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+ 
  theme(plot.background = element_rect(colour = NA)) +
  size

# Arrange the plots together
PerRPE <- ggarrange(PerformanceRPlot, RPERPlot,
                            labels = c("A", "B"),
                            ncol = 2,
                            common.legend = T,
                            legend = "right",
                            align = "h"
)
PerRPE

# Cleveland plot
# adapt y axis 
brks <- c(seq(-900, 900, by = 60))
# Configure labels for the graph
lbls = c(seq(900, 0, -60), seq(60, 900, 60))

# Generate Plot
PerfoIND <- data %>%  # Cast the users table as a number
  mutate(Performance = as.numeric(Performance)) %>%
  # Pipe into ggplot
  ggplot(aes(x = reorder(ID,(Performance)), y = Performance, fill = Condition)) +
  # Plot the point. Alpha to increase transparency of the points and avoid overlapping
  geom_point(aes(color = Condition, alpha =0.5))+
  # join the point (they disappear in the interactive plot)
  geom_line(aes(group = ID), linewidth = 0.1) +
  # Shift the y axis
  scale_y_continuous(breaks = brks, labels = lbls) +
  # Flip the coordinates
  coord_flip() +
  # Add a theme
  sm_classic() +
  # Add a title
  labs(title="Individual Performance", 
       y = "Time to exhaustion (s)", x = "ID",
       caption = "n = 22") +
  # Add more theming options
  theme(plot.title = element_text(hjust = .5),
        axis.ticks = element_blank())+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  size

ggplotly(PerfoIND)

knitr::kable(head(force[, 1:5]), "pipe", caption = "Knee extensor evaluation variables")

# We specify the order that should appear in the graph the time.
# By default they are ordered alphabetically 
# If we do  not specify Post would come first than Pre (alphabetically)
force$Time <- factor(force$Time, levels=c("Pre", "Post Task", "Post TTE"))
force$Condition <- factor(force$Condition, levels=c("Control", "Mental Fatigue"))

# filter a data.frame with only the variable of interest
MVC <-  filter(force, Variable == "MVC")

#We calculate some data to add to the graph. 
sumrepdatMVC <- summarySE(MVC, measurevar = "Value",
                       groupvars=c("Condition", "Time"))

MVCrain <- ggplot(MVC, aes(x = Time, y = Value, 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_point(aes(x = as.numeric(Time)-.15, y = Value, colour = Condition),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Time, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  geom_line(data = sumrepdatMVC, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), linetype = 3)+
  geom_point(data = sumrepdatMVC, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "MVC Force (N)", x = element_blank(),
        title ="Maximal Voluntary Contraction",
        caption = "n = 21")+
size

MVCrain

# VAL rain
VAL <- filter(force, Variable == "VAL")

sumrepdatVAL <- summarySE(VAL, measurevar = "Value",
                       groupvars=c("Condition", "Time"))

VALrain <- ggplot(VAL, aes(x = Time, y = Value, 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_point(aes(x = as.numeric(Time)-.15, y = Value, colour = Condition),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Time, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  geom_line(data = sumrepdatVAL, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), linetype = 3)+
  geom_point(data = sumrepdatVAL, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "Voluntary activation (%)", x = element_blank(),
        title ="Voluntary Activation Level",
        caption = "n = 21")+
  size
VALrain

# Doublet rain
Doublet <-  filter(force, Variable == "Doublet100Hz")
# Drop a missing value
Doublet = drop_na(Doublet)
sumrepdatDoub <- summarySE(Doublet, measurevar = "Value",
                          groupvars=c("Condition", "Time"))

Doubletrain <- ggplot(Doublet, aes(x = Time, y = Value, 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_point(aes(x = as.numeric(Time)-.15, y = Value, colour = Condition),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Time, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  geom_line(data = sumrepdatDoub, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), linetype = 3)+
  geom_point(data = sumrepdatDoub, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "100 Hz paired stimulation (N)", x = element_blank(),
        title ="Potentiated doublet amplitude (100Hz)",
        caption = "n = 20")+
  size


Doubletrain

# Ratio rain
Ratio <- filter(force, Variable == "Ratio 10Hz/100Hz")
# Drop a missing value
Ratio = drop_na(Ratio)
sumrepdatRatio <- summarySE(Ratio, measurevar = "Value",
                           groupvars=c("Condition", "Time"))


Ratiorain <- ggplot(Ratio, aes(x = Time, y = Value, 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_point(aes(x = as.numeric(Time)-.15, y = Value, colour = Condition),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Time, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  geom_line(data = sumrepdatRatio, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), linetype = 3)+
  geom_point(data = sumrepdatRatio, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "Ratio 10Hz/100Hz (%)", x = element_blank(),
        title ="Ratio 10 Hz / 100 Hz",
        caption = "n = 21")+
  size

Ratiorain

# Individual Plot MVC

MVCIND<-ggplot(MVC, aes(x = Time, y = Value))
MVCIND <- MVCIND + geom_line(aes(color = Condition, group = ID)) + geom_point()


# Divide panels for conditions (rows) and participants (columns)
MVCIND <- MVCIND + facet_grid(Condition ~ ID) 


MVCIND <- MVCIND +
  labs (y= "MVC Force (N)",
    x= element_blank(),
    title = "MVC for each participant",
    caption = "n = 21") +
  theme(
    strip.background = element_blank(),
    legend.title = element_blank(),
    strip.text.x = element_blank())+
  sm_minimal()+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 7))+
  easy_remove_legend()


ggplotly(MVCIND, width = 800, height = 800)

# Individual Plot VAL

VALIND<-ggplot(VAL, aes(x = Time, y = Value))
VALIND <- VALIND + geom_line(aes(color = Condition, group = ID)) + geom_point()

# Divide panels for conditions (rows) and participants (columns)
VALIND <- VALIND + facet_grid(Condition ~ ID) 


VALIND <- VALIND +
  labs (y= "Voluntary Activation (%)",
    x= element_blank(),
    title = "VAL for each participant"
    ) +
  theme(
    strip.background = element_blank(),
    legend.title = element_blank(),
    strip.text.x = element_blank())+
  sm_minimal()+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 7))+
  easy_remove_legend()


ggplotly(VALIND, width = 800, height = 800)

# Individual Plot Doublet

DoubletIND<-ggplot(Doublet, aes(x = Time, y = Value))
DoubletIND <- DoubletIND + geom_line(aes(color = Condition, group = ID)) + geom_point()

# Divide panels for conditions (rows) and participants (columns)
DoubletIND <- DoubletIND + facet_grid(Condition ~ ID) 


DoubletIND <- DoubletIND +
  labs (y= "100 Hz paired stimulation (N)", 
    x= element_blank(),
    title = "Potentiated Doublet amplitude (100 Hz) for each participant"
    ) +
  theme(
    strip.background = element_blank(),
    legend.title = element_blank(),
    strip.text.x = element_blank())+
  sm_minimal()+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 7))+
  easy_remove_legend()


ggplotly(DoubletIND, width = 800, height = 800)

# Individual Plot Ratio

RatioIND<-ggplot(Ratio, aes(x = Time, y = Value))
RatioIND <- RatioIND + geom_line(aes(color = Condition, group = ID)) + geom_point()

# Divide panels for conditions (rows) and participants (columns)
RatioIND <- RatioIND + facet_grid(Condition ~ ID) 


RatioIND <- RatioIND +
  labs (y= "Ratio 10 Hz / 100 Hz (%)",
    x= element_blank(),
    title = "Ratio 10 Hz / 100 Hz for each participant"
    ) +
  theme(
    strip.background = element_blank(),
    legend.title = element_blank(),
    strip.text.x = element_blank())+
  sm_minimal()+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 7))+
  easy_remove_legend()


ggplotly(RatioIND, width = 800, height = 800)

# take only the variable avg 20 of the data set.
AVG20<- filter(TMS_long, Variable =="AVG20") %>% convert_as_factor(ID, Condition,Time)
# Exclude participants with excessive noise in the signal
AVG20 <- AVG20 %>%  filter(!ID %in% c("8", "18"))

sumrepdatAVG20 <- summarySE(AVG20, measurevar = "Value",
                          groupvars=c("Condition", "Time"))

TMSrain <- ggplot(AVG20, aes(x = Time, y = Value, 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_point(aes(x = as.numeric(Time)-.15, y = Value, colour = Condition),position = position_jitter(width = .05), size = .25, shape = 20)+
  geom_boxplot(aes(x = Time, y = Value, fill = Condition),outlier.shape = NA, alpha = .5, width = .1, colour = "black")+
  geom_line(data = sumrepdatAVG20, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), linetype = 3)+
  geom_point(data = sumrepdatAVG20, aes(x = as.numeric(Time)+.1, y = Value, group = Condition, colour = Condition), shape = 18) +
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  labs( y = "MEP amplitude (mV)", x = element_blank(),
        title ="Average 20 MEPs",
        caption = "n = 20")

TMSrain

# Plot for individual data for Motor Evoked potentials
AVG20plot<-ggplot(AVG20, aes(x = Time, y = Value ))
AVG20plot <- AVG20plot + geom_line(aes(color = Condition, group = ID)) + geom_point()


# Divide panels for conditions (rows) and participants (columns)
AVG20plot <- AVG20plot + facet_grid(Condition ~ ID) 


AVG20plot <- AVG20plot +
  labs (y= "MEP amplitude (mV)",
    x= element_blank(),
    title = "Average 20 MEPs pre-post task for each participant",
    caption = "n = 20") +
  theme(
    strip.background = element_blank(),
    legend.title = element_blank(),
    strip.text.x = element_blank())+
  sm_minimal()+
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1, size = 8))+
  easy_remove_legend()

ggplotly(AVG20plot)

# Define the list of files with the desired pattern i.e., excel files
# you will need to download the data from OSF and put it in a folder and set that dir.
my_files <- list.files(pattern = "*.xlsx")
my_files

invisible({capture.output({

  # Loop for combining the excel files. 
# Choose the values for each participant from time 0 (row 0) to max time 1800s (the task lasted around 30min)
# I did some cleaning in the original file before, like e.g., removing the first 50 rows (descriptive data)...
# and change the time scale from 0 to 1800sec as original files did not begin at the same numerical time...probably could be done here (see later)...but i found easier it for the number of files i had to handle :)

nirs <- lapply(my_files, function(i){
  x=read_excel(i, sheet = 1)
  x=x [c(0:1800), c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)] # 2nd argument is to select the columns variables. There are 15 col with data
  x$file = i
  x
})

})})

# Bind multiple data frames by row and column
nirs <- nirs %>%
  dplyr::bind_rows(nirs)

# Separate column to create participant and condition
nirs=separate(nirs, file, sep="_", into = c("ID", "Condition"))
# Change label of participants
nirs <- nirs %>%  
  mutate(Condition= recode(Condition, A.xlsx ="Control", B.xlsx ="Mental Fatigue"))  

# change 0 values to NA
  nirs[nirs == 0] <- NA
# Drop all NA values
nirs <- na.omit(nirs)

# Change column names()
nirs <- setNames(nirs, (c("Time", "Tx1,Tx2,Tx3 TSI%", "Tx1,Tx2,Tx3 TSI Fit",
                   "Tx1-O2Hb", "Tx1-HHb", "Tx1-tHb", "Tx1-HbDiff",
                   "Tx2-O2Hb", "Tx2-HHb", "Tx2-tHb", "Tx2-HbDiff",
                   "Tx3-O2Hb", "Tx3-HHb", "Tx3-tHb", "Tx3-HbDiff",
                   "ID", "Condition"
                   )))

# Reorder the columns. Just for aesthetics purpose
nirs <- nirs %>% 
select(ID, Time, Condition, everything()) 

# Calculate the average for the 3 Nirs "lengths" and add to the data.frame
nirs$O2Hbtotal<- ((nirs$`Tx1-O2Hb` + nirs$`Tx2-O2Hb` + nirs$`Tx3-O2Hb`)/3)
nirs$HHbtotal <- ((nirs$`Tx1-HHb` + nirs$`Tx2-HHb` + nirs$`Tx3-HHb`)/3)
nirs$HbDiff <- ((nirs$`Tx1-HbDiff` + nirs$`Tx2-HbDiff` + nirs$`Tx3-HbDiff`)/3)
nirs$tHb <- (nirs$O2Hbtotal + nirs$HHbtotal)

# Create bins of of time (8 blocks of 4 min) (with this, one could easily avoid to adapt the time of each individual file, but i realized late...it's learning by doing and it s-works)

# create the variable to cut the times
time <- nirs$Time
# create 8 bins of equal length 
nirs$Period <-  cut(time, 8, labels = c("t1", "t2", "t3", "t4", "t5", "t6", "t7", "t8"))

#change to long formata
nirs_long <- pivot_longer(nirs, cols = 4:21, names_to = "Variable", values_to = "Values")

# convert long format with the average of each period and variable for each participant
nirs_stat <- nirs_long %>% 
  group_by(ID, Condition, Period, Variable) %>% 
  summarise_at(vars(Values),              # Specify column
               list(average = mean))

# Select only the variable of (possible) interest
target <- c("O2Hbtotal", "HHbtotal", "HHbtotal", "HbDiff", "tHb", "Tx1,Tx2,Tx3 TSI%")


# include these variables in the data frame
nirs_stat <- filter(nirs_stat, Variable %in% target)

# Filter the variables of interest for the figures
O2 <- filter(nirs_stat, Variable == "O2Hbtotal") %>% 
   group_by(Period, Condition) %>% summarise_at(vars(average),           
                                                list(average = mean))

HHB <- filter(nirs_stat, Variable == "HHbtotal") %>% 
  group_by(Period, Condition) %>% summarise_at(vars(average),           
                                               list(average = mean))



O2Plot <- ggplot(na.omit(O2), aes(x = Period, y = average, colour = Condition, group = Condition) ) +
  geom_point() + 
  geom_line()+
  scale_y_continuous(limits = c(-5, 2))+ 
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  easy_remove_x_axis(what = "text")+
  geom_hline(yintercept = 0, linetype='dotted', col = 'gray')+
  labs( y = "Δ oxy-HB  (μm)", x = element_blank(),
        title =" Oxyhaemoglobin",
        subtitle = " ",
        caption = "")+
  size


HHBPlot <- ggplot(na.omit(HHB), aes(x = Period, y = average, colour = Condition, group = Condition) ) +
  geom_point() + 
  geom_line()+
  scale_y_continuous(limits = c(-6, 2))+ 
  scale_colour_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2")+
  sm_classic() +
  easy_remove_x_axis(what = "text")+
  geom_hline(yintercept = 0, linetype='dotted', col = 'gray')+
  labs( y = "Δ deoxy-HB  (μm)", x = element_blank(),
        title =" Deoxyhemoglobin",
        subtitle = " ",
        caption = "")+
  size




NIRSplot <- ggarrange(O2Plot, HHBPlot,
                        labels = c("A", "B"),
                        ncol = 1, nrow = 2, 
                        common.legend = T,
                        legend = "right",
                        align = "h"
)

NIRSplot