---
title: "DataProcessing"
author: "Alex van den Berg"
output:
  html_document:
    toc: true
    theme: united
    highlight: tango
date: "`r Sys.Date()`"
runtime: shiny
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(data.table)
library(RColorBrewer)
library(ggplot2)
library(ggpubr)
library(plyr)
library(readxl)
library(rstatix)
library(plotrix)
library(car)
library(Rmisc)

#LMER stuff
library(lme4)
library(emmeans)
library(sjstats)
library(lmerTest)
library(MuMIn)

# saving figures
library(svglite)
```

## Loading in the Data

### Obtaining the paths

First, we get list of all participants, which are saved into folders in data directory.

```{r}
dataFolder = file.path(".","data")
participants = list.dirs(path = dataFolder, full.names = FALSE, recursive = FALSE)
participants

get_pdir = function(pnum){ 
  return(file.path(dataFolder, pnum))
}
```

The experimental conditions are saved in terms of the [UXF sessions](https://github.com/immersivecognition/unity-experiment-framework/wiki/Concepts), which are the same for each participant. We create a data table to easily translate between the session ID and the condition name

```{r}
participantFolder = get_pdir(participants[1]) # we pick some participant
sessions = list.dirs(path = participantFolder, full.names = FALSE, recursive = FALSE)
conditions = data.table(name=c("A","F","N"),dir=c(sessions[1],sessions[2],sessions[3]))
setkey(conditions,name)
```

We'll need the paths to the [trackers](https://github.com/immersivecognition/unity-experiment-framework/wiki/Data-collection). Which contain all the data which is tracked for each of the trials. The same tracker names are used for each participant for each session (condition in this case), so we'll just get a list we can refer back to later. Additionally, we'll define another function to get a certain tracker dataset for a selected participant.

```{r}
trackerPath = file.path(get_pdir(participants[1]),conditions["A"],"trackers")
trackers = list.files(path = trackerPath, full.names = FALSE)
trackers

# get any type of data
get_tdata = function(pnum,condition,tnum) {
  return(read.csv(file.path(get_pdir(pnum),conditions[condition]$dir,"trackers",trackers[tnum])))
}
# We use these to get the right and left foot data
left = function(pnum,c) {
  return(get_tdata(pnum,c,5))
}
right = function(pnum,c) {
  return(get_tdata(pnum,c,6))
}
```

### Get the target positions

We need the target positions actually defined in Unity to compare to the detected step locations to calculate the distance of the foot placement to the targets.

```{r}
path = file.path(dataFolder,"targetlocations.csv")
targets = read.csv(path)
```

### Get the Distances

These distances need to be extracted from a separate excel sheet for all participants and tucked into a data frame so that we can plot things later.

We add 0's to the single digit numbers in stepVec to make sure they make sense when alphabetically ordered. We need them as characters otherwise some plot methods don't recognize the column as categorical.

```{r, include=FALSE}
# this has a bunch of output I can't really suppress in another way...
distancesPath = file.path(dataFolder,"distances_new.xlsx")
distances <- read_excel(distancesPath,"Distances_corrected") 
step_indices <- read_excel(distancesPath,"Indices") # the index of the step detection in the data
```

```{r, cache=TRUE}
colNames_ = colnames(distances)
distancesTable = data.table(participant=c(),condition=c(),stepNr=c(),distance=c(), index=c(), leg=c())

for (i in 1:length(distances)) {
  split = strsplit(colNames_[i],"_")[[1]]
  currParticipant = split[1]
  cnum = as.numeric(split[[2]])
  
  d = distances[[i]]
  stepVec = as.character(1:16) # vector with step numbers
  stepVec[which(nchar(stepVec)==1)] = paste0("0",stepVec[which(nchar(stepVec)==1)]) # add 0's to single digit numbers to make sure they are ordered correctly when plotting
  c_ = conditions[[1]][cnum]
  t = data.table(participant=currParticipant,condition=c_,stepNr=stepVec,distance=d,index=step_indices[[i]], leg=c("right","left"))
  # Note: Legs are alternating, starting with right as per experiment instructions
  
  distancesTable=rbind(distancesTable,t)
}
setkey(distancesTable,participant)

distancesTable$dist_calc = 0.0
distancesTable$xrel = 0.0
distancesTable$zrel = 0.0
tx = targets$pos_x
tz = targets$pos_z
for (p in participants) {
  for (c in conditions$name) {
    xright = right(p,c)$pos_x
    zright = right(p,c)$pos_z
    right_rows = distancesTable[p,condition==c & leg=="right"]
    right_indices = distancesTable[p]$index[right_rows]
    right_stepnrs = as.numeric(distancesTable[p]$stepNr[right_rows])
    dx = xright[right_indices] - tx[right_stepnrs]
    dz = zright[right_indices] - tz[right_stepnrs]
    d_right = sqrt(dx^2+dz^2)
    distancesTable[p][condition==c][leg=="right"]$xrel = dx
    distancesTable[p][condition==c][leg=="right"]$zrel = dz
    distancesTable[p][condition==c][leg=="right"]$dist_calc = d_right

    xleft = left(p,c)$pos_x
    zleft = left(p,c)$pos_z
    left_rows = distancesTable[p,condition==c & leg=="left"]
    left_indices = distancesTable[p]$index[left_rows]
    left_stepnrs = as.numeric(distancesTable[p]$stepNr[left_rows])
    dx = xleft[left_indices] - tx[left_stepnrs]
    dz = zleft[left_indices] - tz[left_stepnrs]
    d_left = sqrt(dx^2+dz^2)
    distancesTable[p][condition==c][leg=="left"]$xrel = dx
    distancesTable[p][condition==c][leg=="left"]$zrel = dz
    distancesTable[p][condition==c][leg=="left"]$dist_calc = d_left
  }      
}

```

### Get the Condition order and trial number
We also have the **order in which the conditions** are presented, and want to add this to our data table.

```{r}
cOrderPath = file.path(dataFolder,"condition_orders.xlsx")
conditionOrders = data.table(read_excel(cOrderPath))
conditionOrders$first = conditions[conditionOrders$first]$name
conditionOrders$second = conditions[conditionOrders$second]$name
conditionOrders$third = conditions[conditionOrders$third]$name

distancesTable = merge(distancesTable, conditionOrders, by="participant")
```

We can now also derive the trial number for each data point
```{r}
distancesTable$trialNo = "0"

for (i in 1:length(conditions$name)) {
  distancesTable[condition==conditions$name[i]][first==conditions$name[i]]$trialNo = "1"
  distancesTable[condition==conditions$name[i]][second==conditions$name[i]]$trialNo = "2"
  distancesTable[condition==conditions$name[i]][third==conditions$name[i]]$trialNo = "3"
}
```

### Get the Questionnaire Data

The questionnaire data is in a separate file. We extract this data and calculate some resulting scores for the embodiment. We then add the questionnaire data to our distances data table.

```{r}
questionnairePath = file.path(dataFolder,"questionnaires.xlsx")
questionnaires <- data.table(read_excel(questionnairePath))

setnames(questionnaires, "1. I felt as if the virtual feet were my own feet", "Q1_Feet") # ownership
setnames(questionnaires, "2. It felt like I could control the virtual feet as if they were my own", "Q2_Feet") # ownership
setnames(questionnaires, "3. I felt as if my own feet were located where I saw the virtual feet", "Q3_Feet") # ownership
setnames(questionnaires, "4. I felt like I was wearing different shoes from when I came to the laboratory", "Q4_Feet") # Appearance & Response
setnames(questionnaires, "5. I felt as if the virtual  feet I saw when looking in the mirror were my own.", "Q5_Feet") # low relevance question, but related to ownership in original publication
setnames(questionnaires, "6. It felt as if my feet were turning into an avatar (so becoming artificial)", "Q6_Feet") # Appearance
setnames(questionnaires, "7. The movements of the virtual feet were caused by my movements", "Q7_Feet") # Appearance, Agency in original publication.

setnames(questionnaires, "8. I felt as if the virtual body was my own body", "Q1_Body")
setnames(questionnaires, "9. It felt like I could control the virtual body as if it was my own", "Q2_Body")
setnames(questionnaires, "10. I felt as if my own body was located where I saw the virtual body", "Q3_Body")
setnames(questionnaires, "11. I felt like I was wearing different clothes from when I came to the laboratory", "Q4_Body")
setnames(questionnaires, "12. I felt as if the virtual  body I saw when looking in the mirror was my own", "Q5_Body")
setnames(questionnaires, "13. It felt as if my body were turning into an avatar", "Q6_Body")
setnames(questionnaires, "14. The movements of the virtual body were caused by my movements", "Q7_Body")

setnames(questionnaires, "15. Which version did you prefer?", "S1")
setnames(questionnaires, "16. Which version did you feel like you scored best? Regardless of the score?", "S2")

setkey(questionnaires,participant)

questionnaires$embodiment_Feet = (questionnaires$Q1_Feet+questionnaires$Q2_Feet+questionnaires$Q3_Feet+questionnaires$Q4_Feet+questionnaires$Q5_Feet+questionnaires$Q6_Feet+questionnaires$Q7_Feet)/7
questionnaires$embodiment_Body = (questionnaires$Q1_Body+questionnaires$Q2_Body+questionnaires$Q3_Body+questionnaires$Q4_Body+questionnaires$Q5_Body+questionnaires$Q6_Body+questionnaires$Q7_Body)/7

distancesTable = merge(distancesTable,questionnaires,by="participant")
```


## Pre-processing
### Remove NAs
Some steps had to be removed because of drift. We remove those from the dataset first:

```{r}
na.omit(distancesTable, cols=c("distance"), invert = TRUE)[,c("participant","stepNr","condition","distance")]

distancesTable = na.omit(distancesTable, cols=c("distance"))
```
### Get subset of data

With these inputs we may filter the data to only include certain conditions / participants as we wish.

```{r filters, echo=FALSE}
participants = unique(distancesTable$participant)

inputPanel(
  selectInput("cnum", label = "condition",
              choices = conditions$name,  selected = conditions$name,  multiple = TRUE),
  selectInput("pnum", label = "Participant",
              choices = participants, selected = participants, multiple = TRUE)
)

```

```{r, echo=FALSE}
# Filter the DT based on selected inputs
dtFiltered_condition <- reactive({
  dt = distancesTable[condition %in% input$cnum]
  return (dt)
})

dtFiltered_participant <- reactive({
  dt = dtFiltered_condition()
  if (length(input$pnum)>0) {dt = dt[participant %in% input$pnum]}
  return (dt)
})

# Also for the questionnaires
qFiltered <- reactive({
  dt = questionnaires
  if (length(input$pnum)>0) {dt = dt[participant %in% input$pnum]}
  return (dt)
})

```

### Mark outliers

We then mark the steps with extreme distance values, grouping them by condition.

First, we define some methods to find our outliers:

```{r}
filter_outliers = function(dt_,coeff){ 
  quartiles <- quantile(dt_$distance, probs=c(.25, .75), na.rm = TRUE)
  IQR <- IQR(dt_$distance)
  Lower <- quartiles[1] - coeff*IQR
  Upper <- quartiles[2] + coeff*IQR 
   
  return (subset(dt_, dt_$distance > Lower & dt_$distance < Upper))
}

get_outliers = function(dt_,coeff){ 
  quartiles <- quantile(dt_$distance, probs=c(.25, .75), na.rm = TRUE)
  IQR <- IQR(dt_$distance)
  Lower <- quartiles[1] - coeff*IQR
  Upper <- quartiles[2] + coeff*IQR 
  A = subset(dt_, dt_$distance <= Lower)
  B = subset(dt_, dt_$distance >= Upper)
  return (rbind(A,B))
}
```

We identify the outliers & mark them in our data table.

```{r, echo=FALSE}
inputPanel(
  sliderInput("coeff", "Coefficient cutoff for outliers",
              min = 0.5, max = 3.0, value = 3.0, step = 0.1),
  checkboxInput("outliers", "Filter outliers", 
                value=FALSE),
)

renderPrint({
  dt = dtFiltered_participant()
  data_no_outlier <- NULL
  data_outliers   <- NULL
  coeff = input$coeff
  for (c in conditions$name){
    dt_ = dt[condition==c]
    remainder = filter_outliers(dt_,coeff)
    outliers = get_outliers(dt_,coeff)
    print(outliers[,c("participant","stepNr","condition","distance")])
    print( paste("Number of outliers:", length(dt_$participant) - length(remainder$participant)) )
    data_no_outlier = rbind(data_no_outlier, remainder)
    data_outliers   = rbind(data_outliers, outliers)
  }
})
```

We also filter them out if the check box above is ticked.

```{r,echo=FALSE} 
dtFiltered <- reactive({ # THE FINAL DATASET, WITH OUTLIERS REMOVED IF NECESSARY (No outliers were removed in the end)
  dt = dtFiltered_participant()
  dt$outlier = FALSE
  coeff = input$coeff
  for (c in conditions$name){
    outliers = get_outliers(dt[condition==c],coeff)
    nrows = length(outliers$participant)
    for (o in 1:nrows){
      row = outliers[o,]
      dt[row$participant][stepNr==row$stepNr][condition==row$condition]$outlier = TRUE
    }
  }
  if (input$outliers) {return (dt[outlier==FALSE])}
  return (dt)
})
```

### Calculate means
We also calculate the means for each of the columns. This is done in a reactive container so that our filters from the pre-processing gets applied.

```{r, echo=FALSE}
get_summ <- reactive(function(column_name){
  return ( ddply(dtFiltered(), column_name, summarise, grp.mean=mean(distance)) )
})

get_mu <- reactive({
  column_names = colnames(dtFiltered())
  mu <- lapply(column_names, get_summ())
  return (setNames(mu, column_names))
})
```

We additionally want to get the means for each combination of trial number and condition and print them.

```{r, echo=FALSE}
get_summ_combined <- reactive({
  return ( ddply(dtFiltered(), c("condition","trialNo"), summarise, 
                 grp.mean=mean(distance), 
                 grp.sd=sd(distance),
                 grp.count = length(distance)) )
})

renderPrint({ 
  # Calculate the mean distance & sd for each group based on condition and trialNo
  summ = get_summ_combined()
  #print(summ)
  # Iterate through the unique combinations of condition and trialNo and print the mean distance  & sd for each group
  for (i in 1:nrow(summ)) {
    cat("Condition:", summ$condition[i], ", TrialNo:", summ$trialNo[i], ", MeanDistance:", summ$grp.mean[i],", SD:", summ$grp.sd[i],", N:", summ$grp.count[i], "\n")
  }
})
```

## Plot data

To get a feel for our data, we start by plotting out the our data.

To disable certain plots, use the check boxes in their input window (plotting takes a while because data is loaded every time).

Click here to save the figures
```{r, echo=FALSE}
inputPanel(
  actionButton("save", "Save figures")
)
```

### Feet trajectory

```{r feetplot, echo=FALSE}
inputPanel(
  checkboxInput("plotTraj", "Plot trajectories", 
                value=FALSE),
  checkboxInput("targets", "Show targets", 
                value=TRUE),
  checkboxInput("legend", "Show legend", 
                value=TRUE)
)

xLims = 2.0
yLims = 2.0
ratio = xLims / yLims
rTargets = 0.25 # m
figHeight = 600 # px

renderPlot({
  if (!input$plotTraj) {return()}
  
  plot(1,type="n",xlim = c(-xLims , xLims), ylim = c( -yLims, yLims),xlab="x pos (m)", ylab="z pos (m)")
  
  if (input$targets){
    #points(targets$pos_x,targets$pos_z,pch = 10,cex=5,lwd = 2)
    draw.ellipse(targets$pos_x, targets$pos_z, a = rTargets*ratio, b = rTargets)
  }
  
  n = length(input$cnum)
  currCol=n
  legendNames=c()
  legendColors=c()
  
  for (p in input$pnum) {
    currCol=n
    for (c in input$cnum) {
      rightCol = brewer.pal(n = n, name = 'YlOrRd')[currCol]
      leftCol = brewer.pal(n = n, name = 'GnBu')[currCol]
      
      # get all vectors
      xright = right(p,c)$pos_x
      zright = right(p,c)$pos_z
      xleft = left(p,c)$pos_x
      zleft = left(p,c)$pos_z
      
      # plot the detected steps
      right_rows = dtFiltered()[p,condition==c & leg=="right"]
      right_indices = dtFiltered()[p]$index[right_rows]
      left_rows = dtFiltered()[p,condition==c & leg=="left"]
      left_indices = dtFiltered()[p]$index[left_rows]
      
      points(xright[right_indices],zright[right_indices],col=rightCol,cex=3,lwd = 5)
      points(xleft [left_indices], zleft [left_indices], col=leftCol, cex=3,lwd = 5)
      
      # Get trajectories
      lines(xright,zright,col=rightCol)
      lines(xleft,zleft,col=leftCol)
      
      legendNames = append(legendNames,paste("Right, ",c))
      legendNames = append(legendNames,paste("Left, ",c))
      
      legendColors = append(legendColors,rightCol)
      legendColors = append(legendColors,leftCol)
      
      currCol = currCol-1
    }
  }
  if (input$legend) { legend(0, 1, legend=legendNames,col=legendColors,lty=1) }
},width = figHeight * ratio, height = figHeight)

```

and for the height data

```{r, echo=FALSE}
inputPanel(
  numericInput("minHeight", "Min height",
              min = 0, max = 1.5, value = 1, step = 0.1),
  numericInput("maxHeight", "Max height",
              min = 1, max = 2, value = 1.3, step = 0.1),
  numericInput("timeLim", "Max time",
              min = 1, max = 60, value = 25, step = 1)
)

renderPlot({
  if (!input$plotTraj) {return()}
  
  timeLim = 1
  minHeight = 0
  maxHeight = 0
  my_plot = plot(1,type="n",xlim = c(0 , input$timeLim), ylim = c( input$minHeight, input$maxHeight),xlab="time (s)", ylab="y pos (m)")
  
  n = length(input$cnum)
  currCol=n
  legendNames=c()
  legendColors=c()
  
  for (p in input$pnum) {
    currCol=n
    for (c in input$cnum) {
      rightCol = brewer.pal(n = n, name = 'YlOrRd')[currCol]
      leftCol = brewer.pal(n = n, name = 'GnBu')[currCol]
      
      # get all vectors
      yright = right(p,c)$pos_y
      yleft = left(p,c)$pos_y
      time = (right(p,c)$time - right(p,c)$time[1])
      
      # plot the detected steps
      right_rows = dtFiltered()[p,condition==c & leg=="right"]
      right_indices = dtFiltered()[p]$index[right_rows]
      left_rows = dtFiltered()[p,condition==c & leg=="left"]
      left_indices = dtFiltered()[p]$index[left_rows]
      
      points(time[right_indices],yright[right_indices],col=rightCol,cex=3,lwd = 5)
      points(time[left_indices], yleft [left_indices], col=leftCol, cex=3,lwd = 5)
      
      # Get trajectories
      lines(time,yright,col=rightCol)
      lines(time,yleft,col=leftCol)
      
      legendNames = append(legendNames,paste("Right, ",c))
      legendNames = append(legendNames,paste("Left, ",c))
      
      legendColors = append(legendColors,rightCol)
      legendColors = append(legendColors,leftCol)
      
      currCol = currCol-1
    }
  }
  legend(0, 0, legend=legendNames,col=legendColors,lty=1)
})
```

### Heatmaps

In order to get an idea of the distribution of our data and check for anomalies, we plot some heatmaps of our detected steps and their relation to the target.

```{r, echo=FALSE}
inputPanel(
  checkboxInput("plotHeatmaps", "Plot Heatmaps", 
                value=TRUE),
  selectInput("stepNr", label = "Step Number(s)",
              choices = stepVec, selected=stepVec, multiple = TRUE),
  selectInput("leg", label = "Leg",
              choices = c("right","left"), selected=c("right","left"), multiple = TRUE),
)

circleFun <- function(center = c(0,0),r = 1, npoints = 100){
    tt <- seq(0,2*pi,length.out = npoints)
    xx <- center[1] + r * cos(tt)
    yy <- center[2] + r * sin(tt)
    return(data.frame(x = xx, y = yy))
}

renderPlot({
  if (!input$plotHeatmaps) {return()}
  
  circle = circleFun(c(0,0),0.25,100)
  dt = dtFiltered()
  dt = dt[condition %in% input$cnum][leg %in% input$leg][stepNr  %in% input$stepNr]
  
  ggplot() +
    geom_path(data = circle, aes(x = x, y = y), color = "black") +
    geom_point(data = dt[outlier==FALSE], aes(x = xrel, y = zrel,col=condition), fill=rgb(0,0,0,0.2),shape = 21,size=10) + 
    geom_point(data = dt[outlier==TRUE], aes(x = xrel, y = zrel,col=condition), color=rgb(1,0,0,0.2),fill=rgb(1,0,0,0.5), shape = 21,size=10) + 
    xlim(-0.25,0.25) +
    ylim(-0.25,0.25)

}, height = 500, width = 560)


observeEvent(input$save, {
  if (!input$plotHeatmaps) {return()}
    
  circle = circleFun(c(0,0),0.25,100)
  dt = dtFiltered()
  dt = dt[condition %in% input$cnum][leg %in% input$leg][stepNr  %in% input$stepNr]
  
  p = ggplot() +
    geom_path(data = circle, aes(x = x, y = y), color = "black") +
    geom_point(data = dt[outlier==FALSE], aes(x = xrel, y = zrel,col=condition), fill=rgb(0,0,0,0.2),shape = 21,size=5) + 
    geom_point(data = dt[outlier==TRUE], aes(x = xrel, y = zrel,col=condition), color=rgb(1,0,0,0.2),fill=rgb(1,0,0,0.5), shape = 21,size=10) + 
    xlim(-0.25,0.25) +
    ylim(-0.25,0.25) + theme(legend.position = "none")
  ggsave("heatmap.svg", plot = p, device = "svg", height = 3, width = 3, unit = 'in')   
})

```


### Embodiment scores

```{r, echo=FALSE}
embodimentTable = reactive({ # Dynamically retrieve the data for the embodiment scores from the questionnaires table.
  dt = data.table()
  for (i in 1:length(qFiltered()[[1]])) {
    
    tA = data.table(
      participant = qFiltered()$participant[i], 
      condition   = "A",
      embodiment  = qFiltered()$embodiment_Body[i])
    tF = data.table(
      participant = qFiltered()$participant[i], 
      condition   = "F",
      embodiment  = qFiltered()$embodiment_Feet[i])
    
    dt=rbindlist(list(dt,tA))
    dt=rbindlist(list(dt,tF))
  }
  setkey(dt,participant)
  return (dt)
})

renderPlot({
  if (input$group != "None") {
    ggplot(embodimentTable(), aes_string("condition","embodiment")) +
      geom_boxplot() + geom_jitter()
  }
})

observeEvent(input$save, {
  p = ggplot(embodimentTable(), aes_string("condition","embodiment")) +
      geom_boxplot() + geom_jitter()
  ggsave("embodiment_bp.svg", plot = p, device = "svg") #install.packages("svglite")
})

```
It's quite clear there is no significant difference between the two conditions, but just to be complete we quickly calculate the difference using a some statistical tests:

```{r, echo=FALSE}
renderPrint({
  dt = embodimentTable()
  dtA=dt[condition=="A"]
  dtF=dt[condition=="F"]
  # print the means etc
  print("Avatar")
  print(summary(dtA))
  sd1=sd(dtA$embodiment)
  print(sd1)
  print("Feet")
  print(summary(dtF))
  sd2=sd(dtF$embodiment)
  print(sd2)
  # do the t test & print it
  x = dtA$embodiment
  y = dtF$embodiment
  tt_res = t.test(y, x, paired=TRUE)
  print(tt_res)
  
  # nonparametric version
  print(wilcox.test(x, y, paired=TRUE))
})
```

## First look


### Distances

Let's start by plotting the distribution of the distances for each of the conditions. Use the input fields to split the histograms in different ways.

```{r, echo=FALSE}
categories = append(colnames(distancesTable),"None")
categories = categories[-which(categories=="distance")]

inputPanel(
  selectizeInput("group", "Group by", 
                choices = categories, selected = "condition",multiple = FALSE),
  selectizeInput("position", "Bar positions", 
                choices = c("stack","identity","dodge"), selected = "identity",multiple = FALSE),
  selectizeInput("split", "Split by", 
                choices = categories, selected = "None",multiple = FALSE),
  checkboxInput("means", "Show means", 
                value=TRUE),
  numericInput("height", "Height per plot",
              min = 50, max = 500, value = 200, step = 50),
  numericInput("bin", "Bin width",
              min = 0.001, max = 0.5, value = 0.01, step = 0.01)
)

plotCount <- reactive({
  colName = as.character(input$split)
  return ( uniqueN(dtFiltered()[,..colName]) )
})

plotHeight <- reactive( 
  if (input$split=="None") {
    return (input$height)
  }
  else { return (input$height*plotCount()) } 
  )

renderPlot({
  clength = length(input$cnum)
  
  aes = aes_string(x="distance")
  a = 1
  fill="grey"
  if (input$group!="None") { 
    fill="white"
    aes = modifyList(aes,aes_string(col=input$group)) 
    if (input$position == "identity") {a = 1/clength}
  }

  p = ggplot(dtFiltered(), aes) + 
    geom_histogram(binwidth=input$bin,fill=fill,alpha=a, position=input$position)
    
  mu = get_mu()
  if (input$means && input$split != "None") {p = p + geom_vline(mu[input$split][[1]],mapping=aes_string(xintercept="grp.mean",col=input$split),linetype="dashed") }
  if (input$means && input$split == "None" && input$group != "None") {p = p + geom_vline(mu[input$group][[1]],mapping=aes_string(xintercept="grp.mean",col=input$group),linetype="dashed") }

  
  if (input$split != "None") {p = p+facet_grid(sym(input$split))}
  p
  },height=reactive(plotHeight())
)

```

Here we show a boxplot of the distances divided based on condition. We additionally split the datapoints based on trial number.
```{r, echo=FALSE}
renderPlot({
  ggplot(dtFiltered_participant(), aes(x = interaction(trialNo, condition), y = distance, fill = condition)) +
  geom_boxplot() +
  geom_jitter(aes(), width = 0.2)
})

renderPlot({
ggplot(dtFiltered_participant(), aes(x = condition, y = distance, fill = trialNo)) +
  geom_boxplot(aes(group = interaction(condition, trialNo)), outlier.shape = NA) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.13, dodge.width = 0.75)) +
  scale_fill_grey(name = "Trial No", start = 0.2, end = 0.8)  # scale_color_discrete for colors
})

observeEvent(input$save, {
  p = ggplot(dtFiltered_participant(), aes(x = condition, y = distance, fill = trialNo)) +
    geom_boxplot(aes(group = interaction(condition, trialNo)), outlier.shape = NA) +
    geom_jitter(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.7)) +
    scale_fill_grey(name = "Trial No", start = 0.2, end = 0.8)  # scale_color_discrete for colors
  ggsave("distances_bp.svg", plot = p, device = "svg") #install.packages("svglite")
})

```

### Questionnaire Results

Some initial histograms and scatter plots of the questionnaire data:

```{r, echo=FALSE}
inputPanel(
  selectizeInput("showQ", "Show question", 
                choices = colnames(questionnaires), selected = "embodiment_Feet",multiple = FALSE),
  numericInput("binwidth", "Bin Width",
              min = 0, max = 500, value = 1, step = 0.1),
  selectizeInput("positionQ", "Bar positions", 
                choices = c("stack","identity","dodge"), selected = "identity",multiple = FALSE),
  selectizeInput("groupQ", "Group by", 
                choices = append(colnames(questionnaires),"None"), selected = "None",multiple = FALSE),
)


renderPlot({
  aes = aes_string(x=input$showQ)
  a = 1
  fill="grey"
  if (input$groupQ!="None") { 
    fill="white"
    aes = modifyList(aes,aes_string(col=input$groupQ)) 
    if (input$positionQ == "identity") {a = 1/3}
  }
  
  #nb = geom_histogram(fill=fill,alpha=a,position=input$positionQ)
  b = geom_histogram(fill=fill,alpha=a,position=input$positionQ,binwidth=input$binwidth)
  
  #print(questionnaires[,input$showQ])
  p = ggplot(qFiltered(), aes) + b
  p
  
})

# Scatter plot
renderPlot({
  muQ_ = data.table(merge(get_mu()$participant,qFiltered(),by="participant"))
  
  aes = aes_string(x=input$showQ,y="grp.mean",size=2)
  if (input$groupQ!="None") { 
    aes = modifyList(aes,aes_string(col=input$groupQ, shape=input$groupQ))
    muQ_[[input$groupQ]] = as.factor(muQ_[[input$groupQ]])
  }
  
  p = ggplot(muQ_, aes) + labs(y = "Mean distance [m]") +
    geom_point() + 
    geom_smooth(method='lm')
  p
})

```

The fit summary:

```{r, echo=FALSE}
inputPanel(
  limitTo <- selectizeInput("limitTo", "Limit to group", choices = c(), multiple = TRUE)
)

renderPrint({
  if (is.numeric(qFiltered()[[input$showQ]])) {
    print( paste("Mean of",input$showQ,":",mean(qFiltered()[[input$showQ]])) )
    print( paste("SD of",input$showQ,":",sd(qFiltered()[[input$showQ]])) )
    print( paste("Median ",input$showQ,":",median(qFiltered()[[input$showQ]])) )
  }
  else{
    print(qFiltered()[[input$showQ]])
  }
  
  muQ_ = data.table(merge(get_mu()$participant,qFiltered(),by="participant"))
  
  if (input$groupQ != "None") {
    included = which(muQ_[[input$groupQ]] %in% input$limitTo)
    print(included)
    muQ_ = muQ_[included]
    print("Included participants: ")
    print(muQ_$participant)
    print(".............................")
  }
  
  # limit
  summary(lm(formula = as.formula(paste0("grp.mean ~ ",input$showQ)), data = muQ_))
})

reactive({
  options = unique(questionnaires[[input$groupQ]])
  updateSelectInput(session = getDefaultReactiveDomain(),"limitTo",choices = options,selected=options)
  return()
})

```

## Analysis

We first need to check our assumptions

### Normality

```{r, echo=FALSE}
# qq plot
renderPlot({
  dt = dtFiltered()
  p = ggqqplot(dt, x="distance" )
  p
})

```

### Homogeneity of variance
Unfortunately, we do violate the assumption of homogeneity of variance.
```{r}

dtResiduals = reactive({
  dt = dtFiltered()
  mu = get_mu()
  dt$residual = 999 # just fill with some placeholder value
  dt$mean = 999

  for (i in 1:length(mu$condition$grp.mean)){
    row = mu$condition[i,]
    dt[condition==row$condition]$mean = row$grp.mean
    dt[condition==row$condition]$residual = 
      dt[condition==row$condition]$distance - row$grp.mean
  }
  return (dt)
})

renderPrint({
  dt = dtResiduals()

  print(leveneTest(distance ~ condition, dt))
  print(bartlett.test(distance ~ condition, dt))
})

renderPlot({
  dt = dtResiduals()
  ggplot(dt,
         aes(residual)) +
           geom_density(fill = "gray")
})


```

### Linear Mixed-Effects
Because we have unequal sample sizes (due to drift occurring for some steps) we cannot directly use an ANOVA. Since our data (& its residuals) are roughly normally distributed, we can use a linear mixed model regression.

```{r, echo=FALSE}
get_lmer = reactive({
  dt = dtFiltered()
  
  dt$condition = factor(dt$condition, levels = c("A","F","N"))
  summary(dt)
  
  m <- lmer(distance ~ condition + (1 | participant), data=dt, REML=F)
  return (m)
})

renderPrint({
  m = get_lmer()
  print(m)
  print("..........")
  print(anova(m))
  print(eta_sq(m, partial=TRUE))
  #print(summary(partial_eta_squared(m)))
  print(r.squaredGLMM(m))
  m.emm.s = emmeans(m, list(pairwise ~ condition), adjust = "bonferroni")
  print(m.emm.s)
  m.eff = eff_size(m.emm.s, sigma = sigma(m), edf = 561)
  print(m.eff)
})

renderPlot({
  m = get_lmer()
  qqnorm(resid(m))
})
```

We can now update our boxplot with the p-values:
```{r, echo=FALSE}
makeBoxPlots = reactive({
    stat.test <- tibble::tribble(
    ~group1, ~group2, ~p.adj,     ~y.position, ~p.adj.signif,
    
    #"A",    "F",    1.000000e+00, .24, "",
    "A",    "N",    9.164185e-99, .265, "****",
    "F",    "N",    2.917203e-95, .25, "****",
    
  )
  p = ggplot(dtResiduals(), aes_string("condition","distance")) +
    geom_boxplot() + geom_jitter() + stat_pvalue_manual(stat.test,label= "p.adj.signif")#, y.position = 35, step.increase = 0.1,hide.ns = TRUE)
  return (p)
})

renderPlot({
  p = makeBoxPlots()
  p
})

# Filter the DT based on selected inputs
observeEvent(input$save, {
  ggsave("boxplot.svg", plot = makeBoxPlots(), device = "svg")
})
```


### Fixed effects
We evaluate different models to understand the fixed effects of condition and trial number individually.

```{r, echo=FALSE}
renderPlot({
  ggplot(dtFiltered(), aes(x=trialNo,y=distance)) +
    geom_boxplot()# + geom_jitter()
})

renderPlot({
  ggplot(dtFiltered(), aes(x=trialNo,y=distance,fill=condition)) +
    geom_boxplot()# + geom_jitter()
})


observeEvent(input$save, { # to save the figure
  p = ggplot(dtFiltered(), aes(x=trialNo,y=distance,fill=condition)) +
    geom_boxplot()
  ggsave("exposureorder_bp.svg", plot = p, device = "svg", height = 4, width = 5, unit = 'in') #install.packages("svglite")
})

```

```{r, echo=FALSE}
get_lmer_exposure = reactive({ # Dynamically retrieve the lme results depending on the selected fixed effects
  
  if (input$exposureFixed=="trialNo"){
    dt = dtFiltered()[condition  %in% input$exposureCondition]
    dt$trialNo = factor(dt$trialNo, levels = c(1,2,3), labels("1","2","3"))
    m <- lmer(distance ~ trialNo + (1 | participant), data=dt, REML=F)
    return (m)
  }
  else if (input$exposureFixed=="condition")
  {
    dt = dtFiltered()[trialNo %in% input$exposureTrial]
    dt$condition = factor(dt$condition, levels = c("A","F","N"), labels("A","F","N"))
    m <- lmer(distance ~ condition + (1 | participant), data=dt, REML=F)
    return (m)
  }
  
})

inputPanel(
  selectInput("exposureFixed", label = "Fixed Effect",
              choices = c("trialNo","condition"), selected="trialNo"),
  selectInput("exposureCondition", label = "Condition (for trialNo effect)",
              choices = c("A","F","N"), selected=c("N"),multiple = TRUE),
  selectInput("exposureTrial", label = "Trial number (for condition effect)",
              choices = c(1,2,3), selected=c(1), multiple = TRUE)
)

renderPrint({ # print the analysis results depending on the selected effect
  m = get_lmer_exposure()
  print(m)
  print("..........")
  print(anova(m))
  print(eta_sq(m, partial=TRUE))
  #print(summary(partial_eta_squared(m)))
  print(r.squaredGLMM(m))
  if (input$exposureFixed=="trialNo"){ 
    m.emm.s = emmeans(m, list(pairwise ~ trialNo), adjust = "bonferroni")
    m.eff = eff_size(m.emm.s, sigma = sigma(m), edf = 16)
  }
  else if (input$exposureFixed=="condition"){
    m.emm.s = emmeans(m, list(pairwise ~ condition), adjust = "bonferroni")
    m.eff = eff_size(m.emm.s, sigma = sigma(m), edf = 15.9)
  }
  print(m.emm.s)
  print(m.eff)
})
```

### Combined model
To see if there is an (interaction) effect of the exposure order/conditions, we perform an LME analysis like above where we extend our previous model by incorporating both fixed effects: `Distance ∼ condition*trialNo + (1|participant)`

```{r, echo=FALSE}

renderPrint({
  dt = dtFiltered()
  #print(dt)
  dt$trialNo = factor(dt$trialNo, levels = c("1","2","3"))#, labels("1","2","3"))
  dt$condition = factor(dt$condition, levels = c("A","F","N"))#, labels("A","F","N"))
  m <- lmer(distance ~ condition*trialNo + (1 | participant), data=dt, REML=F)
  print(m)
  print("..........")
  print(anova(m))
  print("Partial Eta-squared:")
  print(eta_sq(m, partial=TRUE))
  print(r.squaredGLMM(m))
  
  # Calculate AIC and BIC
  print("AIC:")
  print(AIC(m))
  print("BIC:")
  print(BIC(m))
  # Calculate the random effects variance
  print("random effects variance: ")
  print(VarCorr(m)$participant)
  # Extract the fixed effects coefficients, standard errors, t-values, and p-values
  #print("fixed effects coefficients: ")
  #print(summary_m$coefficients)
  print("--------------")
  ### Post-hoc
  #trialno
  print("pairwise ~ trialNo:")
  m.trialNo.emm.s = emmeans(m, list(pairwise ~ trialNo), adjust = "tukey")
  m.trialNo.eff = eff_size(m.trialNo.emm.s, sigma = sigma(m), edf = 568)
  print(m.trialNo.emm.s)
  print(m.trialNo.eff)
  #condition
  print("pairwise ~ condition:")
  m.condition.emm.s = emmeans(m, list(pairwise ~ condition), adjust = "tukey")
  m.condition.eff = eff_size(m.condition.emm.s, sigma = sigma(m), edf = 568)
  print(m.condition.emm.s)
  print(m.condition.eff)
  #interaction Perform pairwise comparisons with Tukey's HSD adjustment
  print("pairwise ~ condition * trialNo:")
  m.interaction.emm.s = emmeans(m, list(pairwise ~ condition*trialNo), adjust = "tukey")
  m.interaction.eff = eff_size(m.interaction.emm.s, sigma = sigma(m), edf = 33.6)
  print(m.interaction.emm.s)
  print(m.interaction.eff)
  pairwise_comparisons <- pairs(m.interaction.emm.s, adjust = "tukey") # same results
  #print(pairwise_comparisons)
  jt_results <- joint_tests(m.interaction.emm.s, by = NULL)


})
```


```{r, echo=FALSE}
renderPrint({
  m = get_lmer_exposure_combined()
  em = emmeans(m, list(pairwise ~ condition * trialNo), adjust = "tukey")
  pvals <- summary(em)$`pairwise differences of condition, trialNo`
  pvals <- pvals[, "p.value", drop = FALSE]
  # Create a new data frame to store the p-values
  pvals_df <- data.frame(combo = rep(NA, 36), p.value = pvals)
  # Assign the condition and trialNo values to the "combo" column
  pvals_df$combo <- paste(combos$Var1, combos$Var2, sep = "_")
  # Merge the p-values with the original data frame
  dtFiltered_pvals <- merge(dtFiltered(), pvals_df, by = "combo")
  # Add an asterisk to pairs of groups with p-values less than 0.05
  #pvals_combinations$label <- ifelse(pvals_combinations$p.value < 0.05, "*", "")
  print(dtFiltered_pvals)
})

renderPlot({
ggplot(dtFiltered_participant(), aes(x = condition, y = distance, fill = trialNo)) +
  geom_boxplot(aes(group = interaction(condition, trialNo)), outlier.shape = NA) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.13, dodge.width = 0.75)) +
  scale_color_discrete(name = "Trial No")
})
```

Indeed, we see that there is an interaction effect. In the following subsections we explore the data a bit further, but with a more subjective basis. Results should be interpreted from the above model.

### Variability between participants
We can take a look at the data for each participant. Although these results should be interpreted with caution, as they belong to just one participant and therefore lack statistical power.
```{r, echo=FALSE}
### For each participant
renderPlot({
  p = ggplot(dtFiltered(), aes(x=participant,y=distance,fill=trialNo,colour=condition)) +
      geom_boxplot(outlier.alpha=0) + #geom_jitter(aes(colour=condition))+
      geom_point(position=position_jitterdodge(),alpha=0.7) + facet_wrap(~first,scale="free_x")
  
  p
})

observeEvent(input$save, {
  p = ggplot(dtFiltered(), aes(x=participant,y=distance,fill=trialNo,colour=condition)) +
      geom_boxplot(outlier.alpha=0) + #geom_jitter(aes(colour=condition))+
      geom_point(position=position_jitterdodge(),alpha=0.7) + facet_wrap(~first,scale="free_x")
  ggsave("participant_variability_bp.svg", plot = p, device = "svg", height = 4, width = 10, unit = 'in') #install.packages("svglite")
})
```
We can observe that the variability of a participant is quite dependent on which condition they have started with. So to make these results slightly more insightful, we group all participants based on which condition they started with. Note: this is a more subjective analysis and should be interpreted with caution. These figures just show my thought process and all results should be interpreted from the combined model above.

```{r, echo=FALSE}
### For all participants
renderPlot({
  ggplot(dtFiltered(), aes(x=first,y=distance,fill=condition)) +
      geom_boxplot()# + geom_jitter()
})

observeEvent(input$save, {
  p =   ggplot(dtFiltered(), aes(x=first,y=distance,fill=condition)) + geom_boxplot()
  ggsave("participant_variability_grouped_bp.svg", plot = p, device = "svg", height = 4, width = 5, unit = 'in') #install.packages("svglite")
})
```

And for the statistical analysis:

```{r, echo=FALSE}
inputPanel(
  selectInput("firstFixed", label = "Fixed Effect",
              choices = c("first","condition"), selected="first"),
  selectInput("firstFirst", label = "First condition (for condition effect)",
              choices = c("A","F","N"), selected="N"),
  selectInput("firstCondition", label = "Condition (for first condition effect)",
              choices = c("A","F","N"), selected="N"),
)


get_lmer_first = reactive({
  dtF = dtFiltered()[first==input$firstFirst]
  dtC = dtFiltered()[condition==input$firstCondition]
  
  if (input$firstFixed=="first")
  { # see effect between different starting conditions, given a certain condition
    dtC$first = factor(dtC$first, levels = c("A","F","N"), labels("A","F","N"))
    m <- lmer(distance ~ first + (1 | participant), data=dtC, REML=F)
    return (m)
  }
  else if (input$firstFixed=="condition")
  { # see effect between different conditions, given a certain starting condition
    dtF$condition = factor(dtF$condition, levels = c("A","F","N"), labels("A","F","N"))
    m <- lmer(distance ~ condition + (1 | participant), data=dtF, REML=F)
    return (m)
  }
  
})

renderPrint({
  m = get_lmer_first()
  print(m)
  print(anova(m))
  
  print(eta_sq(m, partial=TRUE))
  print(r.squaredGLMM(m))

  if (input$firstFixed=="first"){ 
    m.emm.s = emmeans(m, list(pairwise ~ first), adjust = "bonferroni")
    m.eff = eff_size(m.emm.s, sigma = sigma(m), edf = 16)
  }
  else if (input$firstFixed=="condition"){
    m.emm.s = emmeans(m, list(pairwise ~ condition), adjust = "bonferroni")
    m.eff = eff_size(m.emm.s, sigma = sigma(m), edf = 188)
  }
  print(m.emm.s)
  print(m.eff)
})
```


### Pearson Correlation
```{r}
library(ggstatsplot)
library(ggside)

inputPanel({
  selectInput("par", "Parametric data", 
                choices = c("parametric", "nonparametric","bayes","robust"),         selected = "parametric",     multiple = FALSE)
})

mu_A <- reactive({
  return ( ddply(dtFiltered()[condition=="A"], "participant", summarise, grp.mean=mean(distance)) )
})
mu_F <- reactive({
  return ( ddply(dtFiltered()[condition=="F"], "participant", summarise, grp.mean=mean(distance)) )
})

renderPlot({
  muAQ_ = data.table(merge(mu_F(),qFiltered(),by="participant"))
  
  ggplot(muAQ_,aes(embodiment_Body)) + geom_density(fill = "gray")
})
renderPlot({
  muFQ_ = data.table(merge(mu_F(),qFiltered(),by="participant"))
  ggplot(muFQ_,aes(embodiment_Feet)) + geom_density(fill = "gray")
})

renderPrint({
  muAQ_ = data.table(merge(mu_A(),qFiltered(),by="participant"))
  muFQ_ = data.table(merge(mu_F(),qFiltered(),by="participant"))
  
  # Check normality of averages
  print(shapiro.test(muAQ_$grp.mean))
  print(shapiro.test(muAQ_$embodiment_Body))
  print(shapiro.test(muFQ_$grp.mean))
  print(shapiro.test(muFQ_$embodiment_Feet))
  
  correlationA <- cor(muAQ_$embodiment_Body, muAQ_$grp.mean, method = 'pearson')
  correlationF <- cor(muFQ_$embodiment_Feet, muFQ_$grp.mean, method = 'pearson')
  print(paste("Embodiment Avatar vs Distance:", correlationA))
  print(paste("Embodiment Feet vs Distance:", correlationF))
})

renderPlot({
  muQ_ = data.table(merge(mu_A(),qFiltered(),by="participant"))
  ggscatterstats(data=muQ_,x=embodiment_Body, y = grp.mean, type=input$par)
})
renderPlot({
  muQ_ = data.table(merge(mu_F(),qFiltered(),by="participant"))
  ggscatterstats(data=muQ_,x=embodiment_Feet, y = grp.mean, type=input$par)
})

```

```{r, echo=FALSE}
observeEvent(input$save, {
  muQ_F = data.table(merge(mu_F(),qFiltered(),by="participant"))
  p = ggscatterstats(data=muQ_F,x=embodiment_Feet, y = grp.mean, type=input$par,  title="Relationship between feet embodiment and average distance", conf.level=.1 ,
                     xsidehistogram.args = list(binwidth=0.5, fill="#009E73",color = "black"),
                     ysidehistogram.args = list(binwidth=0.0025,fill="#D55E00",color = "black") )
  ggsave("correlationFeet.svg", plot = p, device = "svg", height = 6, width = 6, unit = 'in') #install.packages("svglite")
  
  muQ_A = data.table(merge(mu_A(),qFiltered(),by="participant"))
  p = ggscatterstats(data=muQ_A,x=embodiment_Body, y = grp.mean, type=input$par, title="Relationship between full-body avatar embodiment and average distance", conf.level=.1 ,
                     xsidehistogram.args = list(binwidth=0.5, fill="#009E73",color = "black"),
                     ysidehistogram.args = list(binwidth=0.0025,fill="#D55E00",color = "black") )
  ggsave("correlationAvatar.svg", plot = p, device = "svg", height = 6, width = 6, unit = 'in')
})
```