---
title: "CHI 2025 - Facial Inputs Study 1"
author: "Mark McGill"
output:
  html_document:
    theme: flatly
    toc: true
    toc_float: 
      collapsed: true
        
date: "`r Sys.Date()`"
---
Analysis for Study 1 - participants performed 3 trials * 53 target facial action units sensed by a Meta Quest Pro, each trial lasting around 3 seconds in the "do/hold this as comfortably and accurately as possible" or some such. In each trial, we logged not just the FAU they were asked to perform, but all the other 52 FAU "scores" (a magic activation number between 0..1 reported by the Meta SDK). 

```{r, results='hide',  echo=FALSE}
library(plyr) 
library(dplyr) 
library(ggplot2)
library(ARTool)
library(scales)
library(ggtext)
library(purrr)
library(tidyverse)
library(readr)
library(patchwork)
```

```{r, results='hide',  echo=FALSE}
ExportPlot <- function(gplot, filename, width=2, height=1.5, doPNG=TRUE, doPDF=TRUE, doEPS=FALSE) {
  # Export plot in PDF and EPS.
  # Notice that A4: width=11.69, height=8.27
  if (doPDF){
    ggsave(paste(filename, '.pdf', sep=""), gplot, width = width, height = height)
  }
  
  if (doEPS){
    postscript(file = paste(filename, '.eps', sep=""), width = width, height = height, horizontal = FALSE, onefile = FALSE, paper = "special")
  }
  #print(gplot)
  #dev.off()
  if (doPNG){
    #png(file = paste(filename, '_.png', sep=""), width = width * 100, height = height * 100)
    ggsave(paste(filename, '_.png', sep=""), gplot, width = width, height = height)
  }
  #print(gplot)
  #dev.off()
}
```


# Parsing the FAU data from Study 1
Lots of csv files split in folders by participant. Common headers we'll use are *PID* for participant ID; and *AUofFocus* for the AU they were performing during the trial/task.

```{r, results='hide',  echo=FALSE}
options(readr.show_col_types = FALSE)

read_resting_face_file <-function(x){
    t <- read.csv(x, header=FALSE) # load file
    names(t) <- c('AUofFocus','Calibration Threshold')
    PID <- parse_number(x)
    t <- t %>% mutate(PID = PID)
    return(t)
}

# This is the participants resting face state captured at the start
resting_face_files <- list.files(path = "Raw Data/",pattern="Calibration_Participant_.*\\.csv$", recursive=TRUE, include.dirs = TRUE, full.names=TRUE)

# This is the answers to performance, fatigue, comfort given after each set of trials
questionnaire_files <- list.files(path = "Raw Data/",pattern="Questionnaire_Participant_.*\\.csv$", recursive=TRUE, include.dirs = TRUE, full.names=TRUE)

# This is the real-time performance data of the FAUs
recording_files <- list.files(path = "Raw Data/",pattern="Recording_Participant_.*\\.csv$", recursive=TRUE, include.dirs = TRUE, full.names=TRUE)

data.questionnaire_raw <- questionnaire_files %>%
  map_dfr( read_csv ) %>% rename(PID = ParticipantID)

 data.questionnaire_raw$PID <- factor(data.questionnaire_raw$PID)

data.resting_face_raw <- resting_face_files %>%
  map_dfr( read_resting_face_file ) 

data.recording_raw <- recording_files %>%
  map_dfr( read_csv ) 

# Should be 1060 rows - 53 AUs * 20 people
data.questionnaire_raw
#nrow(data.questionnaire_raw)
```
Grab the released data as "resting face" data for the model, written to "data.recording_released.csv".

```{r}
data.recording_released <- data.recording_raw %>% filter(TimeInTrial > 3500 & TrialPhase == "released") # %>% select (-c(TrialNo, TrialPhase, TimeInTrial, EndOfTrialFrame, TrialFrameNo))

# get the various IDs we need for later  
FAUs <- levels(factor(data.recording_released$AUofFocus))
PIDs <- levels(factor(data.recording_released$PID))

for( i in rownames(data.recording_released) ){
   data.recording_released[i, "AUofFocusValue"] = data.recording_released[i, data.recording_released[i, "AUofFocus", drop=TRUE], drop=TRUE]
}
write.csv(data.recording_released, "r_outputs/data.recording_released.csv", row.names=FALSE)
data.recording_released
```
# Filter FAU Data by Trial Time
Filter the recording data based on a sensible window within the trial, to account for noise at the start/end of the FAU, and extract the values for the AUofFocus from the massive table of FAU values for simplicity

```{r}
data.recording_filtered <- data.recording_raw %>% filter(TimeInTrial > 500 & TimeInTrial <2500 & TrialPhase == "holding") # %>% select (-c(TrialNo, TrialPhase, TimeInTrial, EndOfTrialFrame, TrialFrameNo))

# get the various IDs we need for later  
FAUs <- levels(factor(data.recording_filtered$AUofFocus))
PIDs <- levels(factor(data.recording_filtered$PID))

for( i in rownames(data.recording_filtered) ){
   data.recording_filtered[i, "AUofFocusValue"] = data.recording_filtered[i, data.recording_filtered[i, "AUofFocus", drop=TRUE], drop=TRUE]
}

head(data.recording_filtered)
write.csv(data.recording_filtered, "r_outputs/data.recording_filtered.csv", row.names=FALSE)

```

# FAU Activation - Resting State
For each AU captured alongside the AUofFocus, we want to calculate the value - resting face value captured at the start of the study
```{r}

data.recording_filtered_resting_face = NULL

# Overwrite the FAU values to be trial value - resting face value
for(i in PIDs ){
  recorded_subset = subset(data.recording_filtered, PID == i)

  for (FAU in FAUs){
    resting_face_subset <- data.resting_face_raw %>% filter(PID == i & AUofFocus == FAU)
    resting_face_default <- resting_face_subset$`Calibration Threshold`
    recorded_subset[paste0(FAU)] <- recorded_subset[paste0(FAU)] - resting_face_default
  }
  
  data.recording_filtered_resting_face <- rbind(data.recording_filtered_resting_face, recorded_subset)
}

# zero any negative values as this ranges from 0..1
#FAUs_RF <- paste(FAUs, "_Minus_RF", sep="")
data.recording_filtered_resting_face[FAUs] <- lapply(data.recording_filtered_resting_face[FAUs], function(x) ifelse(x < 0.0, 0.0, x))

# then overwrite the AUofFocusValue data with the new values
for( i in rownames(data.recording_filtered_resting_face) ){
   data.recording_filtered_resting_face[i, "AUofFocusValue"] = data.recording_filtered_resting_face[i, data.recording_filtered_resting_face[i, "AUofFocus", drop=TRUE], drop=TRUE]
}

# order AUofFocus alphabetically
data.recording_filtered_resting_face$AUofFocus <- factor(data.recording_filtered_resting_face$AUofFocus)

head(data.recording_filtered_resting_face)
write.csv(data.recording_filtered_resting_face, "r_outputs/data.recording_filtered_resting_face.csv", row.names=FALSE)
```

# FAU Activation
Calculate a per-participant per-FAU activation threshold based on the participant's performance data

Reasonable approach is to calculate a confidence interval for the values and pick a threshold based on the lower bound.

```{r}
# standard error
se <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))
} 

# confidence interval
conf <- 0.01 
quantile <- qnorm(1-conf/2)

data.thresholds <- data.recording_filtered_resting_face %>% 
  group_by(PID, AUofFocus) %>%
  select(c(AUofFocusValue)) %>% 
  summarise(across(everything(), list(mean=mean, sd=sd, sum=sum, count=~ n(), max=max, se=se)))

# what's a reasonable way of calculating everyones activation thresholds based on their performance???
# could take a percentage of peak?
data.thresholds$AUofFocus_CalculatedThreshold_Peak <- 0.8 * data.thresholds$AUofFocusValue_max
# or mean - std dev?
data.thresholds$AUofFocus_CalculatedThreshold_Mean <- data.thresholds$AUofFocusValue_mean - data.thresholds$AUofFocusValue_sd
# or base it on a lower CI?
data.thresholds$AUofFocus_CIUpper <- data.thresholds$AUofFocusValue_mean + (quantile* data.thresholds$AUofFocusValue_se)
data.thresholds$AUofFocus_CILower <- data.thresholds$AUofFocusValue_mean - (quantile* data.thresholds$AUofFocusValue_se)

# whatever it is we pick, stick it into this column
data.thresholds$AUofFocus_Threshold <- data.thresholds$AUofFocus_CILower
#data.thresholds$AUofFocus_Threshold <- data.thresholds$AUofFocus_CalculatedThreshold_Mean
data.thresholds <- data.thresholds %>% select(c(PID, AUofFocus, AUofFocus_Threshold))

# can't have thresholds lower than 0, and if a threshold is 0 then the AU is effectively unactivatable...
data.thresholds$AUofFocus_Threshold <- ifelse(data.thresholds$AUofFocus_Threshold < 0.0, 0.0, data.thresholds$AUofFocus_Threshold)

head(data.thresholds)
write.csv(data.thresholds, "r_outputs/data.thresholds.csv", row.names=FALSE)

```

For each AU captured alongside the AUofFocus, calculate whether it met the calibrated activation threshold 
```{r}
data.recording_filtered_with_thresholds = NULL

# hacky, but doing this in a loop as brain is fried
for(i in PIDs ){
  recorded_subset = subset(data.recording_filtered_resting_face, PID == i)

  for (FAU in FAUs){
    threshold_subset <- data.thresholds %>% filter(PID == i & AUofFocus == FAU)
    threshold_value <- threshold_subset$`AUofFocus_Threshold`
    recorded_subset[paste0(FAU, "_Activation")] <- ifelse(recorded_subset[FAU] > threshold_value , 1, 0)
  }
  
  data.recording_filtered_with_thresholds <- rbind(data.recording_filtered_with_thresholds, recorded_subset)
}

head(data.recording_filtered_with_thresholds)
```

And do the same calculation for the AUofFocus (wrote this first, but no harm doing it this way as well - but this should end up just reproducing same results as if you looked up the AUofFocus's Activation column from the previous block...). We'll also join the "AU_categorisation.csv" data at this stage. 

```{r}
categories <- read.csv("Raw Data/AU_categorisation.csv", stringsAsFactors = FALSE)
```


```{r}
# Join the threshold data with the results
data.recording_filtered_with_thresholds <- data.recording_filtered_with_thresholds %>% 
  left_join(., data.thresholds, by = join_by(PID, AUofFocus)) 

# And join the AUofFocus data with our categorisation for plotting based on regions (if we want that? who knows!)
data.recording_filtered_with_thresholds <- data.recording_filtered_with_thresholds %>% 
  left_join(., categories, by = join_by(AUofFocus))

# Calculate AUofFocus activations based on threshold
data.recording_filtered_with_thresholds$AUofFocus_PositiveActivations = ifelse(data.recording_filtered_with_thresholds$AUofFocusValue > data.recording_filtered_with_thresholds$`AUofFocus_Threshold`, 1, 0)
#recording_data_filtered$PositiveActivationDuration = recording_data_filtered$PositiveActivations * 38 # rough average of milliseconds between each frame
head(data.recording_filtered_with_thresholds)

write.csv(data.recording_filtered_with_thresholds, "r_outputs/data.recording_filtered_with_thresholds.csv", row.names=FALSE)
```

# Responsiveness
To examine responsiveness, find the first timepoint for each PID/trial where they hit the threshold for the AUofFocus
```{r}
library(dplyr)
library(data.table)

data.responsiveness = data.recording_filtered_with_thresholds %>%
  group_by(PID, AUofFocus, TrialNo) %>%
  mutate(want = ifelse(row_number()==min(which((AUofFocus_PositiveActivations >= 1) == TRUE)),
                       "yes", NA_character_)) %>% filter(want=="yes")


# this data is missing rows at this point because in some cases they never hit the threshold
# so we need to complete the data 
head(data.responsiveness)

df1 <- expand.grid(PID = PIDs, AUofFocus = FAUs, TrialNo = c(1,2,3), TimeInTrial = 2000)
df1$PID <- factor(df1$PID)
data.responsiveness$PID <-factor(data.responsiveness$PID)

data.responsiveness = data.responsiveness %>% select(c(PID, AUofFocus, TrialNo, TimeInTrial))

data.responsiveness =  data.responsiveness  %>% full_join(df1, by=join_by(PID, AUofFocus, TrialNo))

data.responsiveness = data.responsiveness %>%
  rowwise() %>%
  mutate(TimeToActivation = min(TimeInTrial.x, TimeInTrial.y)) %>%
  mutate(across(everything(), ~ ifelse(is.na(.), 2000, .))) %>%
  select(-c(TimeInTrial.x, TimeInTrial.y))

head(data.responsiveness)

data.recording_filtered_with_thresholds$PID <- factor(data.recording_filtered_with_thresholds$PID)
data.recording_filtered_with_thresholds = data.recording_filtered_with_thresholds %>% left_join(data.responsiveness, by=join_by(PID, AUofFocus, TrialNo))
head(data.recording_filtered_with_thresholds)
```

# % of Positive Activations Over Trial
Now we want to average by participant and AUOfFocus and calculate some percentages for the amount of activations of each FAU over the trials

```{r}
# generate summary stats
data.recording_means <- data.recording_filtered_with_thresholds %>% 
  group_by(PID, AUofFocus) %>%
  select(-c(BodyRegion, Side, TrialPhase,   "AUofFocusSymmetrical")) %>% 
  relocate (c(AUofFocusValue, AUofFocus_PositiveActivations), .after = AUofFocus) %>% # move the new calculated columns to the front before all the interference stuff
  summarise(across(everything(), list(mean=mean, sd=sd, sum=sum, count=~ n())))

# calculate the percentage of time across the three trials the AUofFocus was active
data.recording_means$PositiveActivationPercentage <- (data.recording_means$AUofFocus_PositiveActivations_sum / data.recording_means$AUofFocus_PositiveActivations_count) * 100

# and do the same calculation for all the other FAUs to see what accidental activations occurred
for (FAU in FAUs){
  data.recording_means[paste0(FAU, "_Activation_Percentage")] <- (data.recording_means[paste0(FAU, "_Activation_sum")] / data.recording_means[paste0(FAU, "_Activation_count")]) * 100
}

# figure out the columns needed to sum/count over all FAUs activation summary stats
SumFAUs <- paste(FAUs, "_Activation_sum", sep="")
CountFAUs <- paste(FAUs, "_Activation_count", sep="")

# and then calculate summaries for total accidental activations (excluding the AUofFocus)
data.recording_means <- data.recording_means %>%
  mutate(TotalAccidentalActivations_Sum = rowSums(pick(SumFAUs)), na.rm = TRUE) %>%
  mutate(TotalAccidentalActivations_Count = rowSums(pick(CountFAUs)), na.rm = TRUE)

data.recording_means$AccidentalActivationsPercentage <- ( (data.recording_means$TotalAccidentalActivations_Sum - data.recording_means$AUofFocus_PositiveActivations_sum) / (data.recording_means$TotalAccidentalActivations_Count-data.recording_means$AUofFocus_PositiveActivations_count)  ) * 100

write.csv(data.recording_means, "r_outputs/data.recording_means.csv", row.names=FALSE)

head(data.recording_means)
```

# Combined Dataset and Overall Score
Now we want to create the combined dataset - recorded data (with activation thresholds/responsiveness calculated), questionnaire data (after each set of trials), and calibration data - all grouped by PID and AUofFocus.

Values captured: Average Calibration Value	Calibration Value STD	%Calibration Value STD	Avg Discomfort (lower is better)	Avg Discomfort STD	Avg Borg RPE (lower is better)	Avg Borg RPE STD	Avg Performance	Avg Performance STD	Usability Score (lower is better)

We also calculate our *overall score* here based on fatigue, performance, discomfort, and responsiveness. 

```{r}
data.all_results <- data.recording_means %>%
  # full_join(calibration_data, by = join_by(PID, AUofFocus)) %>%
  full_join(data.questionnaire_raw, by = join_by(PID, AUofFocus)) %>%
  #relocate (c("Calibration Threshold", "Discomfort", "Borg RPE", "Performance", "PositiveActivationPercentage"), .after = PositiveActivationDuration_count)  %>% 
  rename(BorgRPE = "Borg RPE")

# calculate a score to assess the "best" FAUs based on our measures
data.all_results$Score <- data.all_results$BorgRPE + data.all_results$Performance + data.all_results$Discomfort + ((data.all_results$TimeToActivation_mean/2000))*10 + ((1.0-data.all_results$AUofFocusValue_mean)*10)

data.all_results$TimeToActivation_mean = data.all_results$TimeToActivation_mean/1000

data.all_results <- data.all_results %>% 
  left_join(., categories, by = join_by(AUofFocus))

write.csv(data.all_results, "r_outputs/data.all_results.csv", row.names=FALSE)

head(data.all_results)
```

# Plots
## Plots for Responsiveness, Discomfort, Fatigue, Performance, Score
Generate the boxplots for each measure.

```{r, fig.height=8, fig.width=8}
data.plots = data.all_results
data.plots$AUofFocus = factor(data.plots$AUofFocus)
data.plots$AUofFocusSymmetrical = factor(data.plots$AUofFocusSymmetrical)
data.plots = data.plots %>%  mutate(AUofFocusSymmetrical = fct_rev(AUofFocusSymmetrical))



plot_AUofFocusVar <- function(plot_data, measure, reorder=FALSE){
  #plot_data <- data
  #plot_data$AUofFocus = plot_data$AUofFocusSymmetrical
  #plot_data$AUofFocus = as.factor(plot_data$AUofFocus)
  #plot_data <- plot_data %>% mutate(AUofFocus = factor(plot_data$AUofFocus))
  #plot_data$AUofFocus <- factor(plot_data$AUofFocus, levels=rev(levels(plot_data$AUofFocus)))
  
  if (reorder){
    p <- ggplot(plot_data, aes(y=reorder(AUofFocusSymmetrical, .data[[measure]], mean), x=.data[[measure]], colour=AUofFocusSymmetrical)) 
  } else {
      p <- ggplot(plot_data, aes(y=AUofFocusSymmetrical, x=.data[[measure]], colour=AUofFocusSymmetrical)) 
  }
  
  p<- p + geom_boxplot(outlier.colour=NULL, 
                outlier.size=1) + theme(legend.position = "none") + ggtitle(measure) +   stat_summary(fun = mean, geom = "point", size = 2) + geom_point(pch = ".", size = 3)  + theme(plot.margin = margin(10, 10, 10, 10, "pt"))

  print(p)
  return(p)
}


plot.AUofFocusValue_mean = plot_AUofFocusVar(data.plots, "AUofFocusValue_mean") +xlab("Mean Activation") +ggtitle("Mean Activation") +labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

#plot.Thresholds = plot_AUofFocusVar(data.thresholds, "AUofFocus_Threshold") + xlab("Threshold")+ ggtitle("Threshold")

plot.TimeToActivation = plot_AUofFocusVar(data.plots, "TimeToActivation_mean") + xlab("Seconds")+ ggtitle("Responsiveness")+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

plot.PositiveActivationPercentage = plot_AUofFocusVar(data.plots, "PositiveActivationPercentage") + xlim(0,100) + xlab("%")+ ggtitle("Activation %")+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

plot.FalseActivationPercentage = plot_AUofFocusVar(data.plots, "AccidentalActivationsPercentage") + xlim(0,100)+ xlab("%")+ ggtitle("Other FAU\n Activations %")+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

plot.Discomfort = plot_AUofFocusVar(data.plots, "Discomfort")  + scale_x_continuous(breaks=c(0,2.5,5.0,7.5,10), limits=c(0,10), labels=c("Low", "", "", "", "High")) + ggtitle("Perceived\nDiscomfort")+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")
plot.BorgRPE = plot_AUofFocusVar(data.plots, "BorgRPE") + ggtitle("Perceived\nEffort") + xlab("Effort") + scale_x_continuous(breaks=c(0,2.5,5.0,7.5,10), limits=c(0,10), labels=c("Low", "", "", "", "High"))+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

plot.Performance = plot_AUofFocusVar(data.plots, "Performance") + ggtitle("Perceived\nPerformance") + scale_x_continuous(breaks=c(0,2.5,5.0,7.5,10), limits=c(0,10), labels=c("Best", "", "", "", "Worst "))+labs(y="AUs of Focus - Alphabetical Order", colour=" - Alphabetical Order")

plot.Score = plot_AUofFocusVar(data.plots, "Score", reorder=TRUE) + scale_x_continuous(breaks=c(0,10,20,30,40,50), limits=c(0,50), labels=c("Best", "", "", "", "", "Worst   ")) + theme(text = element_text(size = 14)) +labs(y="AUs of Focus - Descending Order", colour="AUs of Focus - Descending Order")


```
## Combi-plot for Paper - Figure 2
Layout 1.
```{r, fig.width=18, fig.height=10}
#combi_plot <- plot.AUofFocusValue_mean  + plot.TimeToActivation + plot.PositiveActivationPercentage  + plot.FalseActivationPercentage + plot.Discomfort +plot.BorgRPE + #plot.Performance +plot.Score + plot_layout(ncol = 9, guides = "collect", axes= "collect")
#ExportPlot(combi_plot, "r_outputs/combi_plot", 14, 9)


layout <- "
ABCD
AEF#
"
combi_plot <- plot.Score + plot.AUofFocusValue_mean  + plot.TimeToActivation  + plot.Discomfort +plot.BorgRPE + plot.Performance + plot_layout(ncol = 3, guides = "collect", axes= "collect", design = layout) 
ExportPlot(combi_plot, "r_outputs/combi_plot", 11, 11)

print(combi_plot)
```

Layout 2.
```{r, fig.width=13, fig.height=7}
combi_plot = plot.AUofFocusValue_mean + plot.TimeToActivation  + plot.Discomfort +plot.BorgRPE + plot.Performance +plot.Score + plot_layout(ncol = 6, guides = "collect", axes= "collect")

ExportPlot(combi_plot, "r_outputs/paper_combi_plot", 13, 7)
print(combi_plot)

```

## Plot - Activation Matrix - Figure 3 Top
Plot by activation percentage.

```{r, fig.width=12, fig.height=10}
#cols_to_select <- paste0(FAUs, "_Activation_Percentage", sep="")
cols_to_select = paste0(FAUs, "_mean", sep="")
cols_to_select <- prepend(cols_to_select, "AUofFocus")
cols_to_select <- prepend(cols_to_select, "PID")

# average everything together
data.corr_plot <- data.all_results %>% 
  select(all_of(c(cols_to_select))) %>% 
  rename_with(~str_remove(., '_mean')) %>%
  group_by(AUofFocus) %>%
  summarise(across(FAUs, list(mean=mean))) %>%
  rename_with(~str_remove(., '_mean'))# %>%
  #select(-c(PID)) 

library(reshape2)

melted_cormat <- melt(data.corr_plot, id=c("AUofFocus")) 
melted_cormat <- rename(melted_cormat, ActivatedAUs = variable)

#head(data.corr_plot)
#head(melted_cormat)

tile_plot_1 <- ggplot(data = melted_cormat, aes(x=AUofFocus, y=ActivatedAUs, fill=value)) + 
  geom_tile() + theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) + scale_fill_viridis_c(direction = -1, option = "plasma", limits=c(0,0.5))+ labs(fill="Mean\nActivation") + xlab("FAU of Focus") + ylab("Activated FAUs")

print(tile_plot_1)
ExportPlot(tile_plot_1, "r_outputs/tile_plot_activation", 10, 8)
# x / along the way tells you if the FAU was robust to false activation
# y / up and down tells you the FAUs that were likely falsely activated when trying that FAU
```
## Plot - Correlation Matrix - Figure 3 Bottom
Calculate correlations between the AUofFOcus and the other FAUs at the time the FAU gesture was being performed. 

```{r}
subset_FAUs <- FAUs # c("MouthLeft", "MouthRight", "CheekPuffL", "CheekPuffR", "JawDrop", "LidTightenerL", "LidTightenerR")

cols_to_select = paste0(subset_FAUs, "_mean", sep="")
cols_to_select <- prepend(cols_to_select, "AUofFocus")
cols_to_select <- prepend(cols_to_select, "PID")


data.corr_plot = data.all_results %>% 
  select(all_of(c(cols_to_select))) %>% 
  rename_with(~str_remove(., '_mean')) %>%
  group_by(AUofFocus, PID) %>%
  summarise(across(subset_FAUs, list(mean=mean))) %>%
  rename_with(~str_remove(., '_mean')) %>% filter(AUofFocus %in% subset_FAUs)

# get the AU of focus value specifically for that row
for( i in rownames(data.corr_plot) ){
   data.corr_plot[i, "AUofFocusValue"] = data.corr_plot[i, data.corr_plot[i, "AUofFocus", drop=TRUE], drop=TRUE]
}



#library(smplot2)
#ggplot(data = data.corr_plot, mapping = aes(x = MouthLeft, y = AUofFocusValue)) +
#  geom_point(shape = 21, fill = "#0f993d", color = "white", size = 3) +
#  sm_statCorr(
#    color = "#0f993d", corr_method = "spearman",
#    linetype = "dashed"
#  ) +facet_wrap(.~AUofFocus)

data.corr_plot
```
```{r}
# convert data to a matrix format for plotting
library(Hmisc)

mat <- outer(1:53, 1:53)
colnames(mat)<- FAUs
rownames(mat)<- FAUs

for( FAUofFocus in FAUs ){
  for (FAUActivated in FAUs){
   data.subset <- subset(data.corr_plot, AUofFocus== FAUofFocus)
   x <- data.subset$AUofFocusValue
   y <- data.subset[FAUActivated]
   
   mycor<-cor(x=x, y=y, method = c("pearson"))
   mat[FAUofFocus, FAUActivated] = mycor
  }
}
# mat
```
```{r, fig.width=12, fig.height=10}
## convert to tibble, add row identifier, and shape "long"
dat2 <-
  mat %>%
  as_tibble(rownames = "Var1") %>%
  pivot_longer(-Var1, names_to = "Var2", values_to = "value") 

tile_plot_2 <- ggplot(data = dat2, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile() + theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) + scale_fill_viridis_c(limits=c(-1, 1), direction = -1, option = "H") + xlab("FAU of Focus") + ylab("Activated FAUs") + labs(fill="Pearson\nCorrelation")

tile_plot_2

ExportPlot(tile_plot_2, "r_outputs/tile_plot_correlation", 10, 8)

```

## Combi-plot for paper - Figure 3
```{r}
tile_combi <- tile_plot_1 + tile_plot_2 + plot_layout(ncol = 1, guides = "collect", axes= "collect")
tile_combi
ExportPlot(tile_combi, "r_outputs/tile_combi", 10, 13)

```


# Statistical Testing
ART Anovas to show differences here by AUofFocus (discomfort, exertion, performance, score) 
```{r} 
data_symmetrical <- data.plots %>% 
  group_by(PID, AUofFocusSymmetrical) %>%
  select(-c(AUofFocus, BodyRegion)) %>% 
  summarise(across(everything(), list(mean=mean, sd=sd)))

fau_art <- function(data, variable){
  
  m <- art(as.formula(paste(variable, "~ AUofFocusSymmetrical + Error(PID)")), data=data)
  m.anova <- anova(m)
  print(m.anova)
  
  contrasts = art.con(m, "AUofFocusSymmetrical") %>%
    summary() %>%
    mutate(sig = ifelse(p.value < 0.05, "*", "")) #%>% 
    #filter(p.value <0.05 & grepl('JawDrop|CheekPuff|MouthLR|LidTightener', contrast) )
  
  # print(contrasts)
  return(contrasts)
} 

```
```{r} 
discomfort_contrasts = fau_art(data_symmetrical, "Discomfort_mean")
```


```{r}
rpe_contrasts = fau_art(data_symmetrical, "BorgRPE_mean")
```

```{r}
score_contrasts = fau_art(data_symmetrical, "Score_mean")
write.csv(score_contrasts, "r_outputs/score_contrasts", row.names=FALSE)
```

# Ranked Optimal FAUs - Table 2
```{r}
library(rPref)
pareto_data  = data_symmetrical %>% group_by(AUofFocusSymmetrical) %>% select(c(AUofFocusSymmetrical, AUofFocusValue_mean_mean, TimeToActivation_mean_mean, Discomfort_mean, BorgRPE_mean, Performance_mean)) %>%   summarise(across(everything(), list(mean=mean)))

result = psel(pareto_data, high(AUofFocusValue_mean_mean_mean) * low(TimeToActivation_mean_mean_mean) * low(Discomfort_mean_mean) * low(BorgRPE_mean_mean) * low(Performance_mean_mean), top=100, show_level=TRUE)

write.csv(result, "r_outputs/rpref_result.csv", row.names=FALSE)

print(result)
```


