# This is supplementary code for the manuscript titled
# "Effectiveness of a Time to Fixate for Fitness to Drive
# Evaluation in Neurological Patients" authored by Nadica
# Miljkovic and Jaka Sodnik shared as preprint on arXiv

# R version
# platform       x86_64-w64-mingw32          
# arch           x86_64                      
# os             mingw32                     
# system         x86_64, mingw32             
# status                                     
# major          4                           
# minor          1.2                         
# year           2021                        
# month          11                          
# day            01                          
# svn rev        81115                       
# language       R                           
# version.string R version 4.1.2 (2021-11-01)
# nickname       Bird Hippie        

# setting the directory to the same path where the table with parameters
# and R code are saved on a computer ("tableParametersAll.csv")
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

# libraries used for the analysis/visualization
library(ggplot2)
library(dplyr)
library(ggpubr)

# reading data
dat <- read.csv("tableParametersAll.csv")

dat$ttf <- as.double(dat$ttf) # time to fixate in ms
dat$fitness <- factor(dat$fitness, levels = 
                           c("fit", "cond fit", "unfit"))

# excluding outliers of reaction times
s1 <- summary(dat$R1)
iq1 <- IQR(dat$R1, na.rm = T)
dat$R1[dat$R1 > (s1[5]+3*iq1)] <- NA
dat$R1[dat$R1 < (s1[2]-3*iq1)] <- NA

s2 <- summary(dat$R2)
iq2 <- IQR(dat$R2, na.rm = T)
dat$R2[dat$R2 > (s2[5]+3*iq2)] <- NA
dat$R2[dat$R2 < (s2[2]-3*iq2)] <- NA

s3 <- summary(dat$R3)
iq3 <- IQR(dat$R3, na.rm = T)
dat$R3[dat$R3 > (s3[5]+3*iq3)] <- NA
dat$R3[dat$R3 < (s3[2]-3*iq3)] <- NA

dat$R1[dat$R1 < 500] <- NA
dat$R1[dat$R1 > 4000] <- NA
dat$R2[dat$R2 < 500] <- NA
dat$R2[dat$R2 > 4000] <- NA
dat$R3[dat$R3 < 500] <- NA
dat$R3[dat$R3 > 4000] <- NA

# averaging repeated measurements
dat$R <- dat$R1 # just to initialize
for (ix in 1:length(dat$R1)){
  dat$R[ix] <- mean(c(dat$R1[ix], dat$R2[ix], dat$R3[ix]),
                    na.rm = T)
}

# excluding outliers in time to fixate (ttf)
# note: please, skip this part if you want get results with outliers
s <- summary(dat$ttf)
s
iq <- IQR(dat$ttf, na.rm = T)
dat$ttf[dat$ttf > (s[5]+3*iq)] <- NA
dat$ttf[dat$ttf < (s[2]-3*iq)] <- NA

dat$ttf[dat$ttf > 500] <- NA # additional rule

# boxplots of ttf and reaction times (PR in the paper)
stplt <- ggplot(data = dat, aes(x = fitness, y = ttf, col = fitness)) + 
  geom_boxplot() + 
  theme_minimal() +
  scale_color_manual(breaks = c("fit", "cond fit", "unfit"),
                     values = c("chartreuse2", "grey50", "red")) +
  ylab("TTF [ms]") +
  xlab("fitness")

rplt <- ggplot(data = dat, aes(x = fitness, y = R, col = fitness)) + 
  geom_boxplot() +
  theme_minimal() +
  scale_color_manual(breaks = c("fit", "cond fit", "unfit"),
                     values = c("chartreuse2", "grey50", "red")) +
  ylab("PR [ms]") +
  xlab("fitness")

sctplt <- ggplot(data = dat,
       aes(x = ttf, y = R, col = fitness)) +
  geom_point() +
  scale_color_manual(breaks = c("fit", "cond fit", "unfit"),
                     values = c("chartreuse2", "grey50", "red"))+
  theme_minimal() +
  ylab("PR [ms]") +
  xlab("TTF [ms]")

ggarrange(sctplt,                                                 
          ggarrange(stplt, rplt, nrow = 2, labels = c("b)", "c)"),
                    label.x = 0.1,
                    legend = FALSE),
          ncol = 2, 
          labels = "a)", common.legend = TRUE,
          label.x = 0.1)

ggsave("fig-ttf-rt.tiff", dpi = 300, units = "in",
       width = 6.5, height = 3.58)

# creating other data frames
fit <- filter(dat, fitness == "fit")
# cond fit stands for conditionally fit
condFit <- filter(dat, fitness == "cond fit")
unFit <- filter(dat, fitness == "unfit")

# calculating summary statistics of all parameters
# summary() function could be used with sd() instead
mean(fit$ttf, na.rm = T) # time to fixate (ttf) in ms
mean(unFit$ttf, na.rm = T)
mean(condFit$ttf, na.rm = T)

sd(fit$ttf, na.rm = T)
sd(unFit$ttf, na.rm = T)
sd(condFit$ttf, na.rm = T)

median(fit$ttf, na.rm = T)
median(unFit$ttf, na.rm = T)
median(condFit$ttf, na.rm = T)

min(fit$ttf, na.rm = T)
min(unFit$ttf, na.rm = T)
min(condFit$ttf, na.rm = T)

max(fit$ttf, na.rm = T)
max(unFit$ttf, na.rm = T)
max(condFit$ttf, na.rm = T)

mean(fit$R, na.rm = T) # average reaction time (R)
mean(unFit$R, na.rm = T)
mean(condFit$R, na.rm = T)

sd(fit$R, na.rm = T)
sd(unFit$R, na.rm = T)
sd(condFit$R, na.rm = T)

median(fit$R, na.rm = T)
median(unFit$R, na.rm = T)
median(condFit$R, na.rm = T)

min(fit$R, na.rm = T)
min(unFit$R, na.rm = T)
min(condFit$R, na.rm = T)

max(fit$R, na.rm = T)
max(unFit$R, na.rm = T)
max(condFit$R, na.rm = T)

mean(fit$ttc, na.rm = T) # time-to-collision (ttc) in s
mean(unFit$ttc, na.rm = T)
mean(condFit$ttc, na.rm = T)

sd(fit$ttc, na.rm = T)
sd(unFit$ttc, na.rm = T)
sd(condFit$ttc, na.rm = T)

min(fit$ttc, na.rm = T)
min(unFit$ttc, na.rm = T)
min(condFit$ttc, na.rm = T)

max(fit$ttc, na.rm = T)
max(unFit$ttc, na.rm = T)
max(condFit$ttc, na.rm = T)

median(fit$ttc, na.rm = T)
median(unFit$ttc, na.rm = T)
median(condFit$ttc, na.rm = T)

mean(fit$speed, na.rm = T) # speed at collision onset in km/h
mean(unFit$speed, na.rm = T)
mean(condFit$speed, na.rm = T)

min(fit$speed, na.rm = T)
min(unFit$speed, na.rm = T)
min(condFit$speed, na.rm = T)

max(fit$speed, na.rm = T)
max(unFit$speed, na.rm = T)
max(condFit$speed, na.rm = T)

sd(fit$speed, na.rm = T)
sd(unFit$speed, na.rm = T)
sd(condFit$speed, na.rm = T)

median(fit$speed, na.rm = T)
median(unFit$speed, na.rm = T)
median(condFit$speed, na.rm = T)

# initial gaze distance (igd) in pixels
# transformed into angles
k <- (82/960 + 52/540)/2 # common scaling
fit$igd <- k*fit$igd
unFit$igd <- k*unFit$igd
condFit$igd <- k*condFit$igd

mean(fit$igd, na.rm = T) 
mean(unFit$igd, na.rm = T)
mean(condFit$igd, na.rm = T)

min(fit$igd, na.rm = T)
min(unFit$igd, na.rm = T)
min(condFit$igd, na.rm = T)

max(fit$igd, na.rm = T)
max(unFit$igd, na.rm = T)
max(condFit$igd, na.rm = T)

sd(fit$igd, na.rm = T)
sd(unFit$igd, na.rm = T)
sd(condFit$igd, na.rm = T)

median(fit$igd, na.rm = T)
median(unFit$igd, na.rm = T)
median(condFit$igd, na.rm = T)

# one-way ANOVA with Bonferroni correction for ttf
pairwise.t.test(dat$ttf, dat$fitness, p.adj = "bonferroni")
# one-way ANOVA with Bonferroni correction for reaction time
pairwise.t.test(dat$R, dat$fitness, p.adj = "bonferroni")

# Additional features tested (optional)
# pairwise.t.test(dat$ttc, dat$fitness, p.adj = "bonferroni")
# pairwise.t.test(dat$speed, dat$fitness, p.adj = "bonferroni")
# pairwise.t.test(dat$igd, dat$fitness, p.adj = "bonferroni")

# Bartlett test to check the variance
bartlett.test(ttf ~ fitness, data = dat)
bartlett.test(R ~ fitness, data = dat)

# p value from the Bartlett test suggests to use Welch's ANOVA for ttf
oneway.test(ttf ~ fitness, data = dat, var.equal = FALSE)
# p value from the Bartlett test suggests to use F test for R
oneway.test(R ~ fitness, data = dat, var.equal = TRUE)

# Additional tests with Welch's ANOVA (optional)
# bartlett.test(igd ~ fitness, data = dat)
# oneway.test(igd ~ fitness, data = dat, var.equal = FALSE)

# Tukey's post hoc test for st
model <- aov(ttf ~ fitness, data = dat)
summary(model)
TukeyHSD(model, conf.level=.95) 
plot(TukeyHSD(model, conf.level=.95), las = 2)
# testing normality of model residuals
shapiro.test(model$residuals)

# Tukey's post hoc test for R
model2 <- aov(R ~ fitness, data = dat)
summary(model2)
TukeyHSD(model2, conf.level=.95) 
plot(TukeyHSD(model2, conf.level=.95), las = 2)
# testing normality of model residuals
shapiro.test(model2$residuals)

# caculating cross-correlation coefficients
cor.test(dat$ttf, dat$speed)
cor.test(dat$ttf, dat$ttc)
cor.test(dat$ttf, dat$igd)
cor.test(dat$ttf, dat$R)

# checking (non)homogeneity (optional)
# library(car)
# leveneTest(dat$st~dat$fitness)
# leveneTest(dat$speed~dat$fitness)
# leveneTest(dat$ttc~dat$fitness)
# leveneTest(dat$igd~dat$fitness)

# multi factor ANOVA
# with interactions, rigor p < 0.001 was considered due to the (non)homogeneity
modelAll <- aov(ttf ~ fitness * speed * ttc * igd, data = dat)
summary(modelAll)
shapiro.test(modelAll$residuals)

# now comparing variables independently (additive general linear model, optional)
modelAll <- aov(ttf ~ fitness + ttc + speed + igd, data = dat)
summary(modelAll)
shapiro.test(modelAll$residuals)

# data visualization (box plots)
ggplot(dat, aes(x = fitness, y = ttf)) + geom_boxplot()
ggplot(dat, aes(x = fitness, y = speed)) + geom_boxplot()
ggplot(dat, aes(x = fitness, y = ttc)) + geom_boxplot()
ggplot(dat, aes(x = fitness, y = igd)) + geom_boxplot()
ggplot(dat, aes(x = fitness, y = R)) + geom_boxplot()

# percents of manual corrections for st calculation
sum(dat$manual_correction, na.rm = T)
# in %
100 * sum(dat$manual_correction, na.rm = T)/length(dat$manual_correction) 

# additional analysis on manual correction

100*sum(dat$manual_correction == 0 & dat$fitness == "fit", na.rm = T)/sum(dat$fitness == "fit")
100*sum(dat$manual_correction == 0 & dat$fitness == "cond fit", na.rm = T)/sum(dat$fitness == "cond fit")
100*sum(dat$manual_correction == 0 & dat$fitness == "unfit", na.rm = T)/sum(dat$fitness == "unfit")

100*sum(dat$manual_correction == 0 & dat$fitness == "fit", na.rm = T)/length(dat$fitness == "fit")
100*sum(dat$manual_correction == 0 & dat$fitness == "cond fit", na.rm = T)/length(dat$fitness == "cond fit")
100*sum(dat$manual_correction == 0 & dat$fitness == "unfit", na.rm = T)/length(dat$fitness == "unfit")