##Set working directory

setwd("C://Users//Russell Chan//OneDrive - Universiteit Twente//2020_MSCA_IF//2_Bachelor_thesis//Step_Emma//Data Analysis//Dataframes//")

##Packages that are required for the analysis

library(readxl)
library(haven)
library(tidyverse)
library(lme4)
library(effects)
library(lattice)
library(car)
library(ggplot2)
library(knitr)
library(reshape2)
library(dplyr)
library(forcats)
library(DHARMa)
library(Hmisc)
library(phia)
library(lsmeans)
library(emmeans)
library(multcomp)
library(nlme)
library(gplots)
library(scales)
library(Rmisc)

##Import dataset from the directory

df.FW <- readxl::read_excel("C://Users//Russell Chan//OneDrive - Universiteit Twente//2020_MSCA_IF//2_Bachelor_thesis//Step_Emma//Data Analysis//Dataframes//df_triallevel_23678_F_W1.xlsx")

##Create factors for analysis

#Subject is a factor to account for intraclass variance in linear models
df.FW$subject <- factor(df.FW$subject)

#Block is our function of time
#df.FW$Block <- factor(df.FW$session)

#Accuracy as a factor
df.FW$accuracy <- factor(df.FW$accuracy)

#Fatigue is a factor to know if it determines RT
#df.FW$fatigue <- factor(df.FW$fatigue)

#Workload is a factor to know if it determines RT
#df.FW$workload <- factor(df.FW$workload)

##Simple linear Models

#1 Fatigue and RT
lmfatigue = lm(RT~fatigue, data = df.FW) #Create the linear regression
summary(lmfatigue) #Review the results
## 
## Call:
## lm(formula = RT ~ fatigue, data = df.FW)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3521 -0.1560 -0.0235  0.1043  8.6861 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.70222    0.02220   31.62   <2e-16 ***
## fatigue     -0.12768    0.01115  -11.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3131 on 1822 degrees of freedom
##   (96 observations deleted due to missingness)
## Multiple R-squared:  0.06713,    Adjusted R-squared:  0.06662 
## F-statistic: 131.1 on 1 and 1822 DF,  p-value: < 2.2e-16
plot(lmfatigue)

#Effects lmfatigue
ae.m.fatigue <- allEffects(lmfatigue)
ae.m.df.fatigue <- as.data.frame(ae.m.fatigue[[1]])

#Fatigue Plot
fatigue <- ggplot(ae.m.df.fatigue, aes(x=fatigue, y=fit)) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
  geom_line() +
  geom_point() +
  ylab("RT (ms)") +
  xlab("fatigue") +
  ggtitle("RT~fatigue") +
  theme_classic()

print(fatigue)

#2 Workload and RT

lmworkload = lm(RT~workload, data = df.FW) #Create the linear regression
summary(lmworkload) #Review the results
## 
## Call:
## lm(formula = RT ~ workload, data = df.FW)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3587 -0.1442 -0.0374  0.0865  8.6499 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.652862   0.018244   35.79   <2e-16 ***
## workload    -0.031588   0.002768  -11.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3132 on 1822 degrees of freedom
##   (96 observations deleted due to missingness)
## Multiple R-squared:  0.06671,    Adjusted R-squared:  0.0662 
## F-statistic: 130.2 on 1 and 1822 DF,  p-value: < 2.2e-16
plot(lmworkload)

#Effects lmworkload
ae.m.workload <- allEffects(lmworkload)
ae.m.df.workload <- as.data.frame(ae.m.workload[[1]])

#Workload Plot
workload <- ggplot(ae.m.df.workload, aes(x=workload, y=fit)) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.1) +
  geom_line() +
  geom_point() +
  ylab("RT (ms)") +
  xlab("workload") +
  ggtitle("RT~worload") +
  theme_classic()

print(workload)

###Building more complex interaction models

#TBC

#Posthocs (Tukey tests for sanity checks. Not really needed.)

#Summary of posthocs
#summary(glht(m.MSL.RT, linfct = mcp(Block = "Tukey")), test = adjusted("holm"))